#include <iostream>
#include <cstdlib>
#include <fstream>
#include <vector>
#include <cmath>
#include <string>
#include <sstream>

#include "fiTQun.h"
#include "fiTQun_shared.h"
#include "spliTChan.h"
#include "PDK_MuGamma.h"

#include "fQROOTOut.h"

#include <TFile.h>
#include <TTree.h>
#include <TVector3.h>

#ifdef NOSKLIBRARIES
#include "WCSimWrap.h"
#else
#include "skheadC.h"
#include "apmringC.h"
#include "apmueC.h"
#include "fillntC.h"
#include "vcworkC.h"
#include "mcguts.h"
#endif

#include "spliTChanOutC.h"
#include "fitqunoutC.h"
#include "TRuntimeParameters.hxx"

#include <time.h>

using namespace fiTQun_parameters;
using namespace std;

extern"C" {
  
#ifndef NOSKLIBRARIES
  void fortinit_(char*, char*, int, int);
  int fortread_();
  void skclosef_(int*);
  
  void trginfo_(float *);
  void vcrdvccm_();
  void mcguts(int *ntracks);
  
  void aprstbnk_(int *);
  void filldstvars_();
  void readfqzbsbank_(int*);
  void fillfqzbsbank_(int*);
#endif
}

const double pi=3.1415926535898;

int main(int argc, char* argv[]){
  
  bool printStuff = false;
  
  int idum;
  
  int ilist[500];
  int nlist=0;
  bool flglst=false;
  
  float trgofst;
  
  time_t timer;
  time_t timeevt;
  time_t tcurrevt;
  
  double tmpnglnL;
  int PCdumflg;
  
  TRuntimeParameters::Get().ParseOptions(argc,argv);
  
  // fiTQun parameter file
  std::string fqParOverrideFileName = TRuntimeParameters::Get().GetFQParOverrideFileName();
  if (fqParOverrideFileName.size() >= 1) {
    TRuntimeParameters::Get().ReadParamOverrideFile(fqParOverrideFileName);
  }
  
  // spliTChan parameter file
  std::string scParOverrideFileName = TRuntimeParameters::Get().GetSCParOverrideFileName();
  if (scParOverrideFileName.size() >= 1) {
    TRuntimeParameters::Get().ReadParamOverrideFile(scParOverrideFileName);
  }
  
  std::string inFileName = TRuntimeParameters::Get().GetInFileName();
  std::string outZBSFileName = TRuntimeParameters::Get().GetOutZBSFileName();
  std::string outROOTFileName = TRuntimeParameters::Get().GetOutROOTFileName();
  std::string eventListFileName = TRuntimeParameters::Get().GetEventListFileName();
  
  if (eventListFileName.size() >= 1) {
    std::cout << "Event list file: " << eventListFileName << std::endl;
    ifstream flist(eventListFileName.c_str());
    while( !flist.eof() ) {
      flist >> ilist[nlist];
      if ((nlist!=0) && !(ilist[nlist]>ilist[nlist-1])) break;
      std::cout << ilist[nlist] << std::endl;
      nlist++;
    }
    flist.close();
    flglst=true;
    std::cout << "Number of events on list: " << nlist << std::endl;
    
  }
  std::cout << "Use event list = " << flglst << std::endl;
  
  int nevents = TRuntimeParameters::Get().GetReadCount();
  std::cout << "Read count = " << nevents << std::endl;
  
  int SkipToEvent = TRuntimeParameters::Get().GetSkipCount();
  int StopatEvent;
  if (SkipToEvent>0) {
    StopatEvent=SkipToEvent+nevents;
  }
  else {
    StopatEvent=nevents;
  }
    
  // Get file name lengths for input into fortinit.
  char* cinFileName = (char*)inFileName.c_str();
  int inFileLength;
  for (int i=0; i<256; i++) {
    if (cinFileName[i]=='\0') {
      inFileLength=i;
      break;
    }
  }
  
  char* coutZBSFileName = (char*)outZBSFileName.c_str();
  int outZBSFileLength;
  for (int i=0; i<256; i++) {
    if (coutZBSFileName[i]=='\0') {
      outZBSFileLength=i;
      break;
    }
  }
  
  int writeHistos = TRuntimeParameters::Get().GetWriteHistos();

  std::cout << "Input file: " << inFileName.c_str() << ";" << inFileLength << std::endl;
  
#ifdef NOSKLIBRARIES
  // setup geometry info, and read first event
  WCSimWrap::Get( (char *)inFileName.c_str() );
  
  // initialize the fit with number of PMTs to use
  fiTQun * thefit = new fiTQun( WCSimWrap::Get()->NPMT() );
  
  spliTChan * theSubEvents = new spliTChan(writeHistos, WCSimWrap::Get()->NPMT());
  
#else
  
  fortinit_((char*)inFileName.c_str(),(char*)outZBSFileName.c_str(),inFileLength,outZBSFileLength);
  
  std::cout << "Output file: " << outZBSFileName.c_str() << ";" << outZBSFileLength << std::endl;
  
  fiTQun *thefit = new fiTQun();
  
  // If memory is not dynamically allocated, program crashes for some reason
  spliTChan *theSubEvents = new spliTChan(writeHistos);
  // Runtime parameters loaded within constructor
  
#endif
  
  fiTQun_shared::Get()->SetWriteHistos(writeHistos);
  int useSubEvents = TRuntimeParameters::Get().GetUseSubEvents();
  int copyTrees = TRuntimeParameters::Get().GetCopyTrees();
  std::string copyTreesList = TRuntimeParameters::Get().GetCopyTreesList();
  int doTOFind    =  TRuntimeParameters::Get().GetParameterI("spliTChan.UseTOFMethod");
  std::cout<<"Use TOF Method? "<<doTOFind<<std::endl;
 
  int fQuiet = TRuntimeParameters::Get().GetParameterI("fiTQun.fQuiet");
  fiTQun_shared::Get()->Cantyoubequiet(fQuiet);

  // Which fits should be run?
  
  bool flgSubEvt = TRuntimeParameters::Get().GetParameterI("fiTQun.DoSubEvent");
  
  bool flg1Rfit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoSingleRingFits");
  
  int flgRingType[nPID];
  flgRingType[ig] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoGamma1RFit");
  flgRingType[ie] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoElectron1RFit");
  flgRingType[imu] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoMuon1RFit");
  flgRingType[ipip] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoPionPlus1RFit");
  flgRingType[ikp] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoKaonPlus1RFit");
  flgRingType[ip] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoProton1RFit");
  flgRingType[icg] = TRuntimeParameters::Get().GetParameterI("fiTQun.DoConeGenerator1RFit");
  
  int flgPDK_MuGamma = TRuntimeParameters::Get().GetParameterI("fiTQun.DoPDK_MuGammaFit");
  PDK_MuGamma *PDK_MuGammafit = NULL;
  if (flgPDK_MuGamma) PDK_MuGammafit = new PDK_MuGamma();
  int PDK_MuGammaSeedType = TRuntimeParameters::Get().GetParameterI("fiTQun.PDKSeedType");
  int PDK_MuGammaFitType = TRuntimeParameters::Get().GetParameterI("fiTQun.PDKFitType");

  int flgpi0fit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoPi0Fit");
  int flgpi0fit2 = TRuntimeParameters::Get().GetParameterI("fiTQun.DoPi0Fit2");
  
  bool flggamma2ElecFit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoGamma2ElecFit");
  bool flgtrueGammaOutput = TRuntimeParameters::Get().GetParameterI("fiTQun.OutputTrueGammaInfo");
  bool flgpipfit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoPiPlusFit");
  
  int doMSFit     =  TRuntimeParameters::Get().GetParameterI("fiTQun.DoMSFit");
  int MSFitMethod =  TRuntimeParameters::Get().GetParameterI("fiTQun.MSFitMethod");
  fiTQun_shared::Get()->SetMSFitMethod(MSFitMethod);
  std::cout<< "Multi-segment fit method: "<<MSFitMethod<<std::endl;
  fiTQun_shared::Get()->SetMSElossMin(TRuntimeParameters::Get().GetParameterD("fiTQun.MSElossMin"));
  fiTQun_shared::Get()->SetMSThetRes(TRuntimeParameters::Get().GetParameterD("fiTQun.MSThetRes"));

  int flgMRfit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoMRFit");
  int flgMoreMRfit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoMoreMRFit");
  
  if (flgRingType[icg]) {
    for (int i=0; i<nPID; i++) flgRingType[i]=0;
    flgRingType[icg]=1;//only run CG hypothesis
    flgpi0fit=0;
    flgpi0fit2=0;
    flggamma2ElecFit=false;
    flgpipfit=false;
    flgMRfit=0;
    flgMoreMRfit=0;
    std::cout << "Cone generator reconstruction!!!" << std::endl;
  }
  std::cout << "flgSubEvt = " << flgSubEvt << std::endl;
  std::cout << "flg1Rfit = " << flg1Rfit << std::endl;
  std::cout << "flgpi0fit = " << flgpi0fit << std::endl;
  std::cout << "flgpi0fit2 = " << flgpi0fit2 << std::endl;
  std::cout << "flgpipfit = " << flgpipfit << std::endl;
  std::cout << "flgtrueGammaOutput: " << flgtrueGammaOutput << std::endl;
  std::cout << "flgMSFit = " << doMSFit << std::endl;
  std::cout << "flgMRfit = " << flgMRfit << std::endl;
  std::cout << "flgMoreMRfit = " << flgMoreMRfit << std::endl;
  
  // Options for all fits
  bool fFitCprof = TRuntimeParameters::Get().GetParameterI("fiTQun.UseFitCProfile");
  if (fFitCprof) std::cout << "Using fit Cherenkov profile!" << std::endl;
  bool fFittprof = TRuntimeParameters::Get().GetParameterI("fiTQun.UseFittProfile");
  if (fFittprof) std::cout << "Using fit time profile!" << std::endl;
  
  int iDetGeom=TRuntimeParameters::Get().GetParameterI("fiTQun.DetGeomType");
  
  thefit->ReadSharedParams(fiTQun_shared::Get()->GetScatflg(),flgRingType,fFitCprof,fFittprof,iDetGeom);
#ifdef NOSKLIBRARIES
  TString strDetName = "WCSim";
#else
  TString strDetName = Form("SK%d",fiTQun_shared::Get()->GetSKVer());
#endif
  
  int iForceFitDefWnd = TRuntimeParameters::Get().GetParameterI("fiTQun.ForceFitDefaultWindow");
  
  bool flgRePrefitTwnd = TRuntimeParameters::Get().GetParameterI("fiTQun.RePrefitTimeWindow");
  int flgInGatefit = TRuntimeParameters::Get().GetParameterI("fiTQun.DoInGateDcyeFit");
    
  if( TRuntimeParameters::Get().HasParameter("fiTQun.CorrTQ") )
    thefit->SetTQCorrflg(TRuntimeParameters::Get().GetParameterI("fiTQun.CorrTQ"));
  
  int flgVtx4Pk = TRuntimeParameters::Get().GetParameterI("fiTQun.PeakSearchVtx");
  
  thefit->SetRingDirTfmFit(TRuntimeParameters::Get().GetParameterI("fiTQun.DoRingDirTfmFit"));
  
  // Read in vertex pre-fit parameters:
  if( TRuntimeParameters::Get().HasParameter("fiTQun.vtxPreFitNfitr") ){
    fiTQun_shared::npfitr = TRuntimeParameters::Get().GetParameterI("fiTQun.vtxPreFitNfitr");
    char sgmarbuff[50];
    if (fiTQun_shared::npfitr > 10) {
      std::cerr << "Error fiTQun.vtxPreFitNfitr must be <= 10." << std::endl;
      exit(-1);
    }
    for (int iNpfitr = 0; iNpfitr < fiTQun_shared::npfitr; iNpfitr++){
      sprintf(sgmarbuff, "fiTQun.vtxPreFitSigmar%i", iNpfitr);
      if( TRuntimeParameters::Get().HasParameter(sgmarbuff) ){
	fiTQun_shared::sgmar[iNpfitr] = TRuntimeParameters::Get().GetParameterD(sgmarbuff);
      } else {
	std::cerr << "Error parameter " << sgmarbuff << " not found. Need " << fiTQun_shared::npfitr << " values of sgmar" << std::endl;
	exit(-1);
      }
    }
  }

  if( TRuntimeParameters::Get().HasParameter("fiTQun.vtxPreFitGrdszVtx") )
    fiTQun_shared::grdszVtx = TRuntimeParameters::Get().GetParameterD("fiTQun.vtxPreFitGrdszVtx");

  if( TRuntimeParameters::Get().HasParameter("fiTQun.vtxPreFitGrdsztime") )
    fiTQun_shared::grdsztime = TRuntimeParameters::Get().GetParameterD("fiTQun.vtxPreFitGrdsztime");
    
  // Options for single ring fits
  int flgtotq[nPID];
  flgtotq[ig] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintGamma1RFit");
  flgtotq[ie] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintElectron1RFit");
  flgtotq[imu] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintMuon1RFit");
  flgtotq[ipip] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintPionPlus1RFit");
  flgtotq[ikp] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintKaonPlus1RFit");
  flgtotq[ip] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintProton1RFit");
  flgtotq[icg] = TRuntimeParameters::Get().GetParameterI("fiTQun.ChargeConstraintConeGenerator1RFit");
  

  thefit->SetQEEffCorr(TRuntimeParameters::Get().GetParameterD(Form("fiTQun.QEEffCorr%s",strDetName.Data())));
  
  int flgTune = TRuntimeParameters::Get().GetParameterI("fiTQun.TuningMode");
  if (flgTune) {
    std::cout << "Using true vector information as initial track parameters!" << std::endl;
  }
  
  // Pi0 fit options
  int gammaSeedType = TRuntimeParameters::Get().GetParameterI("fiTQun.GammaSeedType");
  int pi0SeedType = TRuntimeParameters::Get().GetParameterI("fiTQun.Pi0SeedType");
  int Pi0SeedType2 = TRuntimeParameters::Get().GetParameterI("fiTQun.Pi0SeedType2");
  
  int fPi0fitDcye = 0;
  if (TRuntimeParameters::Get().GetParameterI("fiTQun.Pi0fitOnDcye")) {
    fPi0fitDcye=1;
    std::cout << "Running pi0 fit on the second time window!" << std::endl;
  }
  
  // MR fit options
  int flgMRConvfit = TRuntimeParameters::Get().GetParameterI("fiTQun.FitMRConvL");

  int flgModeSeqFit = TRuntimeParameters::Get().GetParameterI("fiTQun.ModeSeqMRFit");
  
  int flgQTdist = TRuntimeParameters::Get().GetParameterI("fiTQun.OutputQTdist");
  double d_vtx_for_QT = TRuntimeParameters::Get().GetParameterD("fiTQun.dVtxQT");
 
  //Decay E search related
  int doDecayESearch = TRuntimeParameters::Get().GetParameterI("fiTQun.DecayESearch");
  Int_t tmptwl = TRuntimeParameters::Get().GetParameterI("fiTQun.DecayESearchWindowStart");
  Int_t tmptwr = TRuntimeParameters::Get().GetParameterI("fiTQun.DecayESearchWindowEnd");

  TRuntimeParameters::Get().PrintListOfParameters();

  double cqtot = fiTQun_shared::Get()->GetTotqCoef();
  
  // Check if it's an atmospheric MC APFIT file and extract the run/subrun
  string  sinFileName = inFileName.c_str();
  stringstream ssinFileName;
  ssinFileName << inFileName.c_str();
  bool isATMPD = 0;
  bool isMC = 0;
  int mcrun=-1, mcsubrun=-1;
  if (sinFileName.find("atmpd")!=string::npos) isATMPD = 1;
  if (sinFileName.find("mc")!=string::npos) isMC = 1;
#ifndef NOSKLIBRARIES
  if (isATMPD && isMC && sinFileName.find("apfit")!=string::npos) {
    
    vector<string> parsed;
    while(ssinFileName.good()) {
      string substr;
      getline(ssinFileName, substr, '.' );
      parsed.push_back(substr);
    }
    mcrun = atoi(parsed[2].c_str());
    mcsubrun = atoi(parsed[3].c_str());
    
    cout << "Loading atmospheric neutrino MC run " << mcrun << ", subrun " << mcsubrun << endl;
  }
#endif
  
  fiTQun_shared::Get()->Setofilenm(coutZBSFileName);
  
  const int ntestiter=1000;
  
  // ---Variable declarations
  const int nStages = 2;
  const int nPi0Stages = 3;
  
  const int ntwoElecPhotonParams = 8;
  const int ntGammaParams = 13;
  
  // True parameters
  double tGammaParams[ntGammaParams];
  
  // Single track fiTQun parameters
  double snglTrkParams[fiTQun_shared::nSnglTrkParams];
  
  // pi0 fiTQun parameters
  double pi0Params[nPi0Stages][fiTQun_shared::nPi0Params];
  
  //2 Electron (Gamma) parameters
  double twoElecPhotonParams[nStages][ntwoElecPhotonParams];
  double twoElecPhotonnglnL [nStages];
  double twoElecPhotonTotMu [nStages];
  
  enum peakFindingMethods {fqpeak_fqpftof, fqpeak_fq1retof, muechk_notof, muechk_fqpftof, muechk_fq1retof, muechk_aptof};
  
  int nevt=0;
  
  const int nEvisRang = 6;
  Double_t EvisRang[nEvisRang+1] = {0,300,700,1330,2500,5000,15000};
  TH2D *hTQdist[nEvisRang][3] = {};
  
  fQROOTOut *theROOTOut = NULL;
  if (outROOTFileName.size() >= 1) {
    theROOTOut = new fQROOTOut((char*)outROOTFileName.c_str());
      if(copyTrees) theROOTOut->CopyInput( (char *)inFileName.c_str(), copyTreesList );
    theROOTOut->GetTree()->Branch("nevt",&nevt,"nevt/I");
    
#ifndef NOSKLIBRARIES
    
    theROOTOut->GetTree()->Branch("nrun",&apdstnt_.nrun,"nrun/I");
    theROOTOut->GetTree()->Branch("nsub",&apdstnt_.nsub,"nsub/I");
    theROOTOut->GetTree()->Branch("nev",&apdstnt_.nev,"nev/I");
    theROOTOut->GetTree()->Branch("date",apdstnt_.date,"date[3]/I");
    theROOTOut->GetTree()->Branch("time",apdstnt_.time,"time[4]/I");
    
    theROOTOut->GetTree()->Branch("nhit",&apdstnt_.nhit,"nhit/I");
    theROOTOut->GetTree()->Branch("potot",&apdstnt_.potot,"potot/F");
    
    theROOTOut->GetTree()->Branch("nhitac",&apdstnt_.nhitac,"nhitac/I");
    theROOTOut->GetTree()->Branch("wall",&apdstnt_.wall,"wall/F");
    theROOTOut->GetTree()->Branch("evis",&apdstnt_.evis,"evis/F");
    theROOTOut->GetTree()->Branch("nring",&apdstnt_.nring,"nring/I");
    theROOTOut->GetTree()->Branch("ip",apdstnt_.ip,"ip[nring]/I");
    theROOTOut->GetTree()->Branch("pos",apdstnt_.pos,"pos[3]/F");
    theROOTOut->GetTree()->Branch("dir",apdstnt_.dir,"dir[nring][3]/F");
    theROOTOut->GetTree()->Branch("amome",apdstnt_.amome,"amome[nring]/F");
    theROOTOut->GetTree()->Branch("amomm",apdstnt_.amomm,"amomm[nring]/F");
    
    theROOTOut->GetTree()->Branch("nmue",&apntmue_.nmue,"nmue/I");
    theROOTOut->GetTree()->Branch("etype",apntmue_.etype,"etype[nmue]/I");
    theROOTOut->GetTree()->Branch("etime",apntmue_.etime,"etime[nmue]/F");
    theROOTOut->GetTree()->Branch("egood",apntmue_.egood,"egood[nmue]/F");
    theROOTOut->GetTree()->Branch("ehit",apntmue_.ehit,"ehit[nmue]/F");
    
#endif
    
    if (flgQTdist) {
      if (flgQTdist>0) thefit->SetTreeQTdist(theROOTOut->GetTree());
      for (int iErang=0; iErang<nEvisRang; iErang++) {
        for (int iFV=0; iFV<3; iFV++) {
          hTQdist[iErang][iFV] = new TH2D(Form("hTQdist_Erang%d_FV%d",iErang,iFV),"",180,0.,180.,80,-20.,20.);
        }
      }
    }
    
  }
  
  timer = time(NULL);
  
  int nlstctr=0;

  while (1) {
    
    for (int idx=0; idx<20; idx++) {
      fitquninfo_.fqproctime[idx]=0.;
    }
    
    tcurrevt = time(NULL);
#ifndef NOSKLIBRARIES
    
    // Read next ZBS event
    if (nevt>0) {// first event is read in fortinit
      int iret = fortread_();
      if (iret>0) break;
    }
    
    std::cout << std::endl;
    
    std::cout << "Processing ZBS event # " << nevt << std::endl;
    
#else
    
    WCSimWrap * wc = WCSimWrap::Get();
    int iret = wc->LoadEntry( nevt );
    int curevnum = wc->SubEvt(0)->GetHeader()->GetEvtNum();
    
    std::cout <<"================================================="<< std::endl;
    std::cout <<" Load Entry "<<nevt<<" of "<<wc->NEvt()<<" EvtNum "<<curevnum<<" retval="<<iret<<std::endl;
    
    if (iret<0) break;
#endif
    
    if (flglst) {
      if (nlstctr>=nlist) break;
#ifndef NOSKLIBRARIES
      if ( (isATMPD && ilist[nlstctr]==skhead_.nevsk) || (ilist[nlstctr]==nevt) ) {
        nlstctr++;
      } else {
        nevt++;
        continue;
      }
#else
      if ( (isATMPD && ilist[nlstctr]==curevnum) ||
          (ilist[nlstctr]==nevt) ) {
        nlstctr++;
      } else {
        nevt++;
        continue;
      }
#endif
    }
    
    if (SkipToEvent>0) {
      if (nevt<SkipToEvent) {
        nevt++;
        continue;
      }
    }
    if (StopatEvent>0) {
      if (nevt>=StopatEvent) break;
    }
    
#ifndef NOSKLIBRARIES
    /**************Read information from zbs banks************************/
    int inputbnk=0;
    aprstbnk_(&inputbnk);
    filldstvars_();
    
    readfqzbsbank_(&idum);
    
    /* load mc_particle structure */
    int nmctrks;
    mcguts(&nmctrks);
    
    vcrdvccm_();
    trginfo_(&trgofst);
    
    /**************End reading zbs banks**********************************/
#endif
    
    thefit->InitEvent(nevt);
    
    /**************Subevent Section***************************************/
    while(1) {
      
      if (!flgSubEvt) break;
      
      thefit->ClrInfo(0);
      
      if (flgTune>=2) break;
      
      thefit->SetTimeWindow(0);

      // fiTQun pre-fit vertex on default time window
      double vtx_prefit[4] = {0};
      thefit->DoVtxPrefit(vtx_prefit);
      
      // Do fiTQun 1R e-hypothesis fit for better vertex
      if (flgRingType[ie] && flgVtx4Pk!=0) {
        if (flgtotq[ie]) fiTQun_shared::Get()->SetTotqCoef(cqtot);
        else fiTQun_shared::Get()->SetTotqCoef(0.);
        tmpnglnL=thefit->Do1RFit(ie,snglTrkParams,PCdumflg,1);
        
        // TOF corrected fiTQun/spliTChan peak search with FQ 1Re vertex
        thefit->DoPeakSearch(snglTrkParams);
        theSubEvents->CopyPeakInfo(fqpeak_fq1retof,fitquntwnd_.fqtwnd_npeak[0],fitquntwnd_.fqtwnd_peakt0[0], fitquntwnd_.fqtwnd_peakiness[0]);
      }
      
      // TOF corrected fiTQun/spliTChan peak search with FQ prefit vertex
      thefit->DoPeakSearch(vtx_prefit);
      theSubEvents->CopyPeakInfo(fqpeak_fqpftof,fitquntwnd_.fqtwnd_npeak[0],fitquntwnd_.fqtwnd_peakt0[0], fitquntwnd_.fqtwnd_peakiness[0]);
      
      // Search for clusters with simple scanning algorithm.
      // This fills the Cluster_* variables, which is used for shifting the time
      // window below with "FillTQbanks".
      theSubEvents->Cluster_Search(0);
      
      // Search for decay-e peaks based on original APFIT::Muechk algorithm.
      theSubEvents->MuechkWrap(0,muechk_notof); // No TOF correction
      
      // Use default SW trigger window for prefit and 1Re fit below
      // to mimic behaviour of original Muechk
      
      // TOF corrected Muechk with FQ prefit vertex
      theSubEvents->MuechkWrap(vtx_prefit,muechk_fqpftof);
      
      // TOF corrected Muechk with FQ 1Re vertex
      if (flgRingType[ie] && flgVtx4Pk!=0)
        theSubEvents->MuechkWrap(snglTrkParams, muechk_fq1retof);
      
      // Count the number of peaks found by Muechk in each cluster
      theSubEvents->Cluster_AssociateMuechk();
      cout<<"find clusters using peaks..?";
      if (doTOFind) cout<<"YES!"<<endl;
      if (doTOFind){
        theSubEvents->MakeClusters(fitquntwnd_.fqtwnd_npeak[0],
                                   fitquntwnd_.fqtwnd_peakt0[0],vtx_prefit,nevt);
      }
      
#ifndef NOSKLIBRARIES
      if (fiTQun_shared::Get()->GetSKVer()<=3) theSubEvents->SetATMClusters();
#endif
      
      for (int icluster=0; icluster<theSubEvents->Cluster_nCand; icluster++) {
        fitqunse_.cluster_npeaks[icluster][fqpeak_fqpftof]=0;
      }
      
      // Associate peaks to clusters
      for (int ipeak=0; ipeak<fitquntwnd_.fqtwnd_npeak[0]; ipeak++) {
        vtx_prefit[3]=fitquntwnd_.fqtwnd_peakt0[0][ipeak];
        for (int icluster=0; icluster<theSubEvents->Cluster_nCand; icluster++) {
          
          if (thefit->IsPeakInCluster(vtx_prefit,theSubEvents->Cluster_tStart[icluster],
                                      theSubEvents->Cluster_tEnd[icluster])){//is the peak caused by hits in this cluster?
            
            fitqunse_.muechk_icluster[fqpeak_fqpftof][ipeak] = icluster;
            
            fitqunse_.cluster_ipeak[icluster][fqpeak_fqpftof][fitqunse_.cluster_npeaks[icluster][fqpeak_fqpftof]] = ipeak;
            fitqunse_.cluster_timeofpeak[icluster][fqpeak_fqpftof][fitqunse_.cluster_npeaks[icluster][fqpeak_fqpftof]]=fitquntwnd_.fqtwnd_peakt0[0][ipeak];
            fitqunse_.cluster_npeaks[icluster][fqpeak_fqpftof]++;
            
            break;
          }
        }
      }
      
      // Tag bad clusters as those without a corresponding FQ peak
      if (flgRingType[ie] && flgVtx4Pk==2) theSubEvents->RemoveNoPeakClusters(fqpeak_fq1retof);
      else theSubEvents->RemoveNoPeakClusters(fqpeak_fqpftof);
      
      fiTQun_shared::Get()->ClearPeaks(0);
      
      break;
    }
    /**********End Subevent Section***************************************/
    
    fitquninfo_.fqproctime[1] = time(NULL)-tcurrevt;
    
    fitquntwnd_.trgoff=trgofst;
    
    // Sometimes cluster is not found, or does not exist,
    // but run fiTQun on original TQREAL window anyways
    int ncluster_min = 0;
    if (iForceFitDefWnd) ncluster_min = 1;
    int ncluster = max(ncluster_min,fitqunse_.cluster_ncand);
    
    int icluster_good=0;
    for (int icluster=0; icluster<ncluster; icluster++) {
      
      bool noGoodClusters = 0;
      if (!fitqunse_.cluster_goodflag[icluster]) {
        if (!theSubEvents->fQuiet) cout << "Skipping bad cluster " << icluster << endl;
        
        if (icluster==ncluster-1 && icluster_good==0)
          noGoodClusters = 1;
        else
          continue;
      }
      //      if (!(fitqunse_.cluster_npeaks[icluster][0]>1)) continue;//comment this out!!!!!!!!!!!!!!!!!!!!!!!!!
      bool flgCluster=true;
      // Modify time window and offset
      if (useSubEvents && fitqunse_.cluster_ncand && noGoodClusters==0) {
      	
      	// Shift fiTQun time window and offset
      	thefit->SetTimeWindow(icluster_good,fitqunse_.cluster_tstart[icluster],
                              fitqunse_.cluster_tend[icluster]);
      } 
      else {
        flgCluster=false;
      	thefit->SetTimeWindow(icluster_good);
      }
      
      if (flgSubEvt) {
        if (flgCluster) {
          fitquntwnd_.fqtwnd_iclstr[icluster_good]=icluster;
          thefit->DoIngateFit(flgInGatefit,flgRePrefitTwnd,fitqunse_.cluster_npeaks[icluster][fqpeak_fqpftof],fitqunse_.cluster_timeofpeak[icluster][fqpeak_fqpftof]);
        }
        else {
          fitquntwnd_.fqtwnd_iclstr[icluster_good]=-1;
          thefit->DoIngateFit(flgInGatefit,flgRePrefitTwnd,0,NULL);
        }
      }
      else {
        thefit->MaskIngates();
      }
      
      if (flgTune==5) fitqun1r_.fq0rnll[icluster_good]=thefit->FitMu();
      
      int pidtmp = 11;
      int iVctTmp = 0;
      
      bool flgcallPID=false;
      for (int iPID=0; iPID<nPID; iPID++) {
        if (!flg1Rfit) continue;
        
        if (!flgRingType[iPID] || (!(iPID==ie || iPID==imu || iPID==ipip) && icluster_good>0)) {
          continue;
        }
        
        timeevt = time(NULL);
        
        int flg1Rseed;
        
        if (flgTune!=0) {//use true information for initial track parameters
          
#ifndef NOSKLIBRARIES
          pidtmp = abs(vcwork_.ipvc[iVctTmp]);
          if (pidtmp==12 || pidtmp==14 || pidtmp==16) {// neutrino event!
            iVctTmp = 2;// outgoing lepton!
            pidtmp = abs(vcwork_.ipvc[iVctTmp]);
            std::cout << "Neutrino event! Outgoing lepton is: " << pidtmp << std::endl;
          }
#else
          pidtmp = 13;
#endif
          
          snglTrkParams[3]=trgofst;//time
          double momtmp[3],trkdir[3];
          for (int j=0; j<3; j++) {
#ifndef NOSKLIBRARIES
            snglTrkParams[j]=vcwork_.posivc[iVctTmp][j];//vertex
            momtmp[j]=vcwork_.pvc[iVctTmp][j];//momentum (MeV/c)
#else
            snglTrkParams[j]=0.0;//vcwork_.posivc[0][j];//vertex
            momtmp[j]=350.0;//vcwork_.pvc[0][j];//momentum (MeV/c)
#endif
          }
          snglTrkParams[6]=sqrt(momtmp[0]*momtmp[0]+momtmp[1]*momtmp[1]+momtmp[2]*momtmp[2]);
          for (int j=0; j<3; j++) trkdir[j]=momtmp[j]/snglTrkParams[6];
          snglTrkParams[4]=acos(trkdir[2]);//theta
          snglTrkParams[5]=atan2(trkdir[1],trkdir[0]);
          if (fiTQun_shared::PIDarr[iPID]==48) snglTrkParams[6]/=10.;
          snglTrkParams[7]=0.;//should be different for pi+
          
          if (flgTune>=2) {
            flg1Rseed=-1;//do not fit!
          }
          else {
            flg1Rseed=0;
          }
        }
        else {
          if (!flgcallPID || iPID==ie) {
            flgcallPID=true;
            flg1Rseed=1;
          }
          else {
            flg1Rseed=2;
          }
          
        }
        
        if (flgtotq[iPID]) fiTQun_shared::Get()->SetTotqCoef(cqtot);
        else fiTQun_shared::Get()->SetTotqCoef(0.);
        
        if (flgTune==3) {//time testing mode
          thefit->Get1Rmudist(iPID,snglTrkParams);
          int ncall=4000;
          std::cout << ncall << " calls..." << std::endl;
          for (int iLoop=0; iLoop<ncall; iLoop++) {
            tmpnglnL=thefit->GetOneRngnglogL(iPID,snglTrkParams,PCdumflg);
          }
        }
        else if (flgTune==4) {// Fit constants
          thefit->FitQEEff(pidtmp,snglTrkParams);// PID based on MC truth
          break;
        }
        else {
          tmpnglnL=thefit->Do1RFit(iPID,snglTrkParams,PCdumflg,flg1Rseed);
        }
        std::cout << "Time=" << time(NULL)-timeevt << std::endl << std::endl;
        
      }
      
      if (icluster_good==0) fitquninfo_.fqproctime[2] = time(NULL)-tcurrevt;
      
      if (icluster_good==0){
        if (flg1Rfit && doMSFit) {
          
          if (flgtotq[imu]) fiTQun_shared::Get()->SetTotqCoef(cqtot);
          else fiTQun_shared::Get()->SetTotqCoef(0.);
          
          timeevt = time(NULL);
          std::cout<<"Doing MS Stand-alone Fit"<<std::endl;
          thefit->DoMSFit(0,1);
          std::cout << "MS Stand-alone fit total: " << time(NULL)-timeevt << "s" << std::endl << std::endl;
        }
        fitquninfo_.fqproctime[3] = time(NULL)-tcurrevt;
      }
      
      fiTQun_shared::Get()->SetTotqCoef(0.);
      
      if (flg1Rfit && flgRingType[imu]) {
        double dum1RParams[fiTQun_shared::nSnglTrkParams]={0.,0.,0.,900.,1e-6,0.,0.,0.};
        fitqun1r_.fq0rnll[icluster_good]=thefit->GetOneRngnglogL(imu,dum1RParams,PCdumflg);
        fitqun1r_.fq0rtotmu[icluster_good]=thefit->GetTotmu();
      }
      
      
#ifndef NOSKLIBRARIES
      if (flgtrueGammaOutput) {
        
        //For Gamma fit
        float gammaMom [50];
        float gammaPhi [50];
        float gammaTheta [50];
        float gammaX [50];
        float gammaY [50];
        float gammaZ [50];
        //        float gammaT [50];
        for (int i = 0; i < 50; i++) {
          gammaMom[i] = 0;
          gammaPhi[i] = 0;
          gammaTheta[i] = 0;
          //          gammaT[i] = 0;
          gammaX[i]=0;
          gammaY[i]=0;
          gammaZ[i]=0;
        }
        int gammaList[50];
        int nGamma = 0;
        
        for (int itrk=0; itrk<nmctrks; itrk++) {
          int parentid =mc_particle.parent[itrk];
          int pdgCode = mc_particle.type[itrk];
          float momentum = mc_particle.momentum[itrk];
          float phi = atan2(mc_particle.dir[itrk][1],mc_particle.dir[itrk][0]);
          float theta = acos(mc_particle.dir[itrk][2]);
          float inivtx[3] = {0,0,0};
          //          float time = trgofst;
          for (int ix=0; ix<3; ix++) {
            inivtx[ix] = mc_particle.initial_vertex[itrk][ix];
          }
          if (pdgCode == 22) {
            gammaMom[nGamma] = momentum;
            gammaList[nGamma] = itrk;
            gammaPhi[nGamma] = phi;
            gammaTheta[nGamma] = theta;
            gammaX[nGamma] = inivtx[0];
            gammaY[nGamma] = inivtx[1];
            gammaZ[nGamma] = inivtx[2];
            //            gammaT[nGamma] = time;
            for (int iGamma=nGamma-1;iGamma>=0;iGamma--) {
              if (momentum > gammaMom[iGamma]) {
                gammaMom[iGamma + 1] = gammaMom[iGamma];
                gammaList[iGamma + 1] = gammaList[iGamma];
                gammaTheta[iGamma + 1] = gammaTheta[iGamma];
                gammaPhi[iGamma + 1] = gammaPhi[iGamma];
                gammaX[iGamma + 1] = gammaX[iGamma];
                gammaY[iGamma + 1] = gammaY[iGamma];
                gammaZ[iGamma + 1] = gammaZ[iGamma];
                //                gammaT[iGamma + 1] = gammaT[iGamma];
                //                gammaT[iGamma] = time;
                gammaX[iGamma] = inivtx[0];
                gammaY[iGamma] = inivtx[1];
                gammaZ[iGamma] = inivtx[2];
                gammaMom[iGamma] = momentum;
                gammaList[iGamma] = iGamma;
                gammaTheta[iGamma] = theta;
                gammaPhi[iGamma] = phi;
              } else {
                break;
              }
            }
            nGamma ++;
            if (nGamma >= 50) break;
          }
        }
        
        int electronCounter = 0;
        std::cout << "This event has " << nmctrks << " true tracks" << std::endl;
        
        for (int itrk=0; itrk<nmctrks; itrk++) {
          int parentid = mc_particle.parent[itrk];
          int inOutOrSecondary = mc_particle.origin[itrk];
          int pdgCode = mc_particle.type[itrk];
          float mass = mc_particle.mass[itrk];
          float momentum = mc_particle.momentum[itrk];
          float energy = mc_particle.energy[itrk];
          float beta = mc_particle.beta[itrk];
          float dir[3];
          float inivtx[3];
          float finvtx[3];
          for (int ix=0; ix<3; ix++) {
            dir[ix] = mc_particle.dir[itrk][ix];
            /* NUANCE outgoing particles get mislabelled as incoming, so if the
             initial_vertex is not set, use the the_vertex instead */
            inivtx[ix] = mc_particle.initial_vertex[itrk][ix]; // so this is broken for Nuance input files
            finvtx[ix] = mc_particle.final_vertex[itrk][ix];
          }
          float time = mc_particle.time[itrk];
          std::cout << "track number " << itrk << " has PDG code, momentum, parent = " << pdgCode << ", " << momentum << ", " << parentid << std::endl;
          std::cout << "direction = (" << dir[0] << "," << dir[1] << "," << dir[2] << ")" << std::endl;
          std::cout << "inivtx = (" << inivtx[0] << "," << inivtx[1] << "," << inivtx[2] << ")" << std::endl;
          std::cout << "finvtx = (" << finvtx[0] << "," << finvtx[1] << "," << finvtx[2] << ")" << std::endl;
          std::cout << "time, trgofst, mass, beta = " << time << ", " << trgofst << ", " << mass << ", " << beta << std::endl;
          
        }
        
        tGammaParams[0] = gammaX[0];
        tGammaParams[1] = gammaY[0];
        tGammaParams[2] = gammaZ[0];
        tGammaParams[3] = trgofst;
        tGammaParams[4] = gammaTheta[0];
        tGammaParams[5] = gammaPhi[0];
        tGammaParams[6] = gammaMom[0];
        std::cout << "###########True parameters: ";
        for (int i=0; i<ntGammaParams; i++) {
          std::cout << tGammaParams[i] << ", ";
        }
        std::cout << "###############" << std::endl;
      }
#endif
      
      if (flggamma2ElecFit) {
        double gammafitpars[fiTQun_shared::nPi0Params];
        
        if (gammaSeedType == 1) { // seed with truth
          // True Gamma
          // {X, Y, Z, T, theta, phi, p, Etheta, Ephi, Ep, Ptheta, Pphi, Pp}
          // {0, 1, 2, 3, 4,     5,   6, 7,      8,    9,  10,     11,   12}
          // 2 Electron Gamma Fit
          //      0  1  2  3  4       5     6   7
          // X = {X, Y, Z, T, theta, phi, p1, p2}
          if (flgtrueGammaOutput) {
            for (int i=0; i<7; i++) {
              gammafitpars[i] = tGammaParams[i];
            }
            gammafitpars[6] /= 2;
            gammafitpars[7] = gammafitpars[6];
          } else {
            std::cerr << "ERROR: Seed gamma fit with truth requested,"
            << "but truth determination turned off." << std:: endl;
          }
        } else if (gammaSeedType == 2) { // scan with 1-ring e-like fit fixed
          // start with 1-ring e-like fit result
          for (int ipar=0; ipar<fiTQun_shared::nSnglTrkParams; ipar++) {
            gammafitpars[ipar] = snglTrkParams[ipar];
          }
          twoElecPhotonnglnL[0] = thefit->SetTwoColinearRingSeed(gammafitpars);
          twoElecPhotonTotMu[0] = thefit->GetTotmu();
        } else {
          std::cerr << "ERROR, unknown gammaSeedType: " << gammaSeedType << std::endl;
        }
        
        std::cout << "###########Initial parameters: ";
        for (int i=0; i<ntwoElecPhotonParams; i++) {
          std::cout << gammafitpars[i] << ", ";
          twoElecPhotonParams[0][i] = gammafitpars[i];
        }
        std::cout << "###############" << std::endl;
        
        if (fiTQun_shared::Get()->GetDrawflg()) thefit->DrawPreddist();
        
        /*
         if (flgpfit) {
         int fixAllButP1P2[nPi0Params];
         for (int i=0; i<nPi0Params; i++) {
         if ((i==6)||(i==10)) {
         fixAllButP1P2[i] = 0;
         } else {
         fixAllButP1P2[i] = 1;
         }
         }
         
         twoElecPhotonnglnL [1] = thefit->TwoColinearRingFit(1,1,fitpars,PCdumflg,fixAllButP1P2);
         
         std::cout << "###########After PFit parameters: ";
         for (int i=0; i<nPi0Params; i++) {
         std::cout << fitpars[i] << ", ";
         }
         std::cout << "###############" << std::endl;
         }
         */
        
        twoElecPhotonnglnL[1] = thefit->TwoColinearRingFit(1,1,gammafitpars,PCdumflg);
        twoElecPhotonTotMu[1] = thefit->GetTotmu();
        if (fiTQun_shared::Get()->GetDrawflg()) thefit->DrawPreddist();
        
        std::cout << "###########        Fit result: ";
        for (int i=0; i<ntwoElecPhotonParams; i++) {
          std::cout << gammafitpars[i] << ", ";
          twoElecPhotonParams[1][i] = gammafitpars[i];
        }
        std::cout << "###############" << std::endl;
        
        std::cout << "nglnL = " << twoElecPhotonnglnL[1] << std::endl;
        
      }
      
      if (icluster_good==fPi0fitDcye) {
        
        if (flgpi0fit) {
          
          double tmpSeedPi0Params[fiTQun_shared::nPi0Params];
          double tmpPi0Params[fiTQun_shared::nPi0Params];
          
          // 1st pi0 fit
          bool constrainPi0Mass = false;
          if (flgpi0fit == 2) constrainPi0Mass = true;
          tmpnglnL = thefit->DoPi0Fit(ie,ie,constrainPi0Mass,pi0SeedType,tmpPi0Params,PCdumflg,0,tmpSeedPi0Params);
          
          for (int ipar=0; ipar<fiTQun_shared::nPi0Params; ipar++) {
            pi0Params[0][ipar] = tmpSeedPi0Params[ipar];
            pi0Params[1][ipar] = tmpPi0Params[ipar];
          }
          
        }
        fitquninfo_.fqproctime[4] = time(NULL)-tcurrevt;
        
        if (flgpi0fit2) {
          
          double tmpSeedPi0Params[fiTQun_shared::nPi0Params];
          double tmpPi0Params[fiTQun_shared::nPi0Params];
          
          int tmpPi0SeedType2 = Pi0SeedType2;
          if ((Pi0SeedType2 == pi0SeedType)&&(flgpi0fit)) {
            tmpPi0SeedType2 += 100000000;
            for (int ipar=0; ipar<fiTQun_shared::nPi0Params; ipar++) {
              tmpSeedPi0Params[ipar] = pi0Params[0][ipar];
            }
          }
          // pi0 fit 2
          bool constrainPi0Mass = false;
          if (flgpi0fit2 == 2) constrainPi0Mass = true;
          tmpnglnL = thefit->DoPi0Fit(ie,ie,constrainPi0Mass,tmpPi0SeedType2,tmpPi0Params,PCdumflg,1,tmpSeedPi0Params);
          
          for (int ipar=0; ipar<fiTQun_shared::nPi0Params; ipar++) {
            pi0Params[2][ipar] = tmpPi0Params[ipar];
          }
          
        }
        fitquninfo_.fqproctime[5] = time(NULL)-tcurrevt;
        
      }
      
      if (icluster_good==0) {
        
        if (flgMRfit) {
          timeevt = time(NULL);
          thefit->DoMRFit(flgMRConvfit,flgMRfit);
          std::cout << "MR fit total: " << time(NULL)-timeevt << "s" << std::endl << std::endl;
        }
        fitquninfo_.fqproctime[6] = time(NULL)-tcurrevt;
        
        int tmpTArr[10];
        if (flgMoreMRfit) {
          timeevt = time(NULL);
          thefit->DoMoreMRFit(flgMoreMRfit,flgModeSeqFit,doMSFit,tmpTArr);
          std::cout << "MR fit improvement total: " << time(NULL)-timeevt << "s" << std::endl << std::endl;
        }
        fitquninfo_.fqproctime[7] = fitquninfo_.fqproctime[6]+tmpTArr[0];
        fitquninfo_.fqproctime[8] = fitquninfo_.fqproctime[6]+tmpTArr[1];
        fitquninfo_.fqproctime[9] = fitquninfo_.fqproctime[6]+tmpTArr[2];
	
	if (flgPDK_MuGamma){
          PDK_MuGammafit->SetTimeWindow(0);
          if (useSubEvents && fitqunse_.cluster_ncand && noGoodClusters==0) {
            // Shift fiTQun time window and offset
            PDK_MuGammafit->SetTimeWindow(icluster_good,fitqunse_.cluster_tstart[icluster],
                                  fitqunse_.cluster_tend[icluster]);
          } 
          else {
            flgCluster=false;
            PDK_MuGammafit->SetTimeWindow(icluster_good);
          }
          std::cout << "Starting PDK_MuGamma fitter" << std::endl;
	  PDK_MuGammafit->DoFit(nevt, icluster_good, PDK_MuGammaSeedType, PDK_MuGammaFitType);
	}
      }
      
      if (icluster_good==0) {
        if (flgQTdist) {
#ifndef NOSKLIBRARIES
          double wallVal = apdstnt_.wall;
          double EneVal = apdstnt_.evis;
          
          bool fNoAPInfo = false;
          if (apdstnt_.nring==0) {// use truth info!
            fNoAPInfo = true;
            wallVal = apmcnt_.wallv;
            EneVal = apmcnt_.pmomv[iVctTmp];
          }
          
          int iErang;// evis range!
          for (iErang=0; iErang<nEvisRang-1; iErang++) {
            if (EneVal<EvisRang[iErang+1]) break;
          }
          
          int iFV = 0;
          if (wallVal>500.) iFV = 1;
          if (wallVal>1000.) iFV = 2;
          
          TH2D *hTQdist_tmp = NULL;
          if ((fNoAPInfo && apdstnt_.nhitac<16 && wallVal>200.) || (apdstnt_.nhitac<16 && apdstnt_.wall>200. && apdstnt_.evis>100. && apdstnt_.nring==1 && apdstnt_.ip[0]==2 && apntmue_.nmue==0)) {
            hTQdist_tmp = hTQdist[iErang][iFV];
          }
          if (pidtmp!=11) hTQdist_tmp = NULL;// for truth seeding, only plot true electrons!
          thefit->GetQTdist_1Re(d_vtx_for_QT,hTQdist_tmp);
#endif
        }
      }
      
      icluster_good++;
      
    } // End icluster loop
    
    fitquntwnd_.fqntwnd = icluster_good;


	if(doDecayESearch)
	{
	thefit->SetTimeWindow(fitquntwnd_.fqntwnd, tmptwl, tmptwr);
	fitquntwnd_.fqtwnd_iclstr[fitquntwnd_.fqntwnd]=-2;		//only for ignored event
	fitquntwnd_.fqntwnd++;
	thefit->DoDecayESearch();	//the three paremeter will be updated. the information will contain previous peak (not only electron).
	}
    
    fitquninfo_.fqproctime[0] = time(NULL)-tcurrevt;
    std::cout << "Processing time for event # " << nevt << " = " << fitquninfo_.fqproctime[0] << "s" << std::endl;
    
    nevt++;
    
#ifndef NOSKLIBRARIES
    // Set atmospheric MC run/subrun numbers
    if (isATMPD && isMC) {
      skhead_.nrunsk = mcrun;
      skhead_.nsubsk = mcsubrun;
    }
#endif
    
    if (theROOTOut) theROOTOut->Fill();
    
#ifndef NOSKLIBRARIES
    int fillstatus = outZBSFileName.size();
    fillfqzbsbank_(&fillstatus);
#endif
    
  }
  
  if (theROOTOut) {
    if (flgQTdist) {
      theROOTOut->GetFile()->cd();
      for (int iErang=0; iErang<nEvisRang; iErang++) {
        for (int iFV=0; iFV<3; iFV++) {
          hTQdist[iErang][iFV]->Write();
        }
      }
    }
    
    theROOTOut->SaveTree();
    
    delete theROOTOut;
  }
  
#ifndef NOSKLIBRARIES
  int flun = 10;
  skclosef_(&flun);
  if (outZBSFileName.size() >= 1) {
    flun = 20;
    skclosef_(&flun);
  }
#endif
  
  std::cout << "Total time elapsed: " << time(NULL)-timer << "s" << std::endl;
  
  std::cout << "runfiTQun finished" << std::endl;
  
}