#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TSystemDirectory.h"
#include "TSystemFile.h"
#include "TObjString.h"
#include "TGraphErrors.h"
#include "RooRealVar.h"
#include "RooHistPdf.h"
#include "RooPolyVar.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooPlot.h"
#include "RooDataSet.h"

#include "RooGaussian.h"
#include "RooFFTConvPdf.h"

#include <ROOT/RDataFrame.hxx>
#include "ROOT/RVec.hxx"
#include <ROOT/RLogger.hxx>
#include <ROOT/RSnapshotOptions.hxx>

#include <iostream>
#include <stdlib.h>
#include <vector>
#include <map>
#include <fstream>
#include <cassert>

using namespace std;
using namespace RooFit;
using namespace ROOT;
using namespace ROOT::VecOps;

#define RVecD ROOT::RVec<double>

//ECAL[2][5][6] cell
int ecalID(int iplane, int ix, int iy) {
  int ret = 0;
  ret = iplane * 30 + ix * 6 + iy;
  return ret;
}

//HCAL[4][3][3]
int hcalID(int iplane, int ix, int iy) {
  int ret = 0;
  ret = iplane * 9 + ix * 3 + iy;
  return ret;
}

/* Identify period from run number
 */
UShort_t GetPeriodFromRunN(const unsigned int &runN) {
  UShort_t period = 0;
  if (runN >= 12427 && runN <= 12510) period = 1;
  else if (runN >= 12511 && runN <= 12664) period = 2;
  else if (runN >= 12664 && runN <= 12777) period = 3;
  else if (runN >= 12778 && runN <= 12920) period = 4;
  else if (runN >= 12921 && runN <= 13089) period = 5;
  else if (runN >= 13090 && runN <= 13183) period = 6;
           
  if (period == 0) {
    cout << "ERROR: Period not identified for run " << runN << endl;
    cout << "ERROR: exiting... " << endl;
    exit(1);
  }

  return period;
}

/*The following function is used to define in-time energy for ECAL cells.
 */
ROOT::RVec<double> EcalEneInTime(RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma, int nsigma) {
  ROOT::RVec<double> ECAL_eneT;
  for (int is = 0; is < 2; is++) {
    for (int ix = 0; ix < 5; ix++) {
      for (int iy = 0; iy < 6; iy++) {
        auto id = ecalID(is, ix, iy);
        auto deltaT = t0[id] - tE[id];
        if (fabs(deltaT) < nsigma * tSigma[id])
          ECAL_eneT.push_back(ene[id]);
        else
          ECAL_eneT.push_back(0.);
      }
    }
  }
  return ECAL_eneT;
}

/*The following function is used to define in-time energy for HCAL cells.
 */
ROOT::RVec<double> HcalEneInTime(RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma, int nsigma) {

  ROOT::RVec<double> HCAL_eneT;
  for (int is = 0; is < 4; is++) {
    for (int ix = 0; ix < 3; ix++) {
      for (int iy = 0; iy < 3; iy++) {
        auto id = hcalID(is, ix, iy);
        auto deltaT = t0[id] - tE[id];
        if (fabs(deltaT) < nsigma * tSigma[id])
          HCAL_eneT.push_back(ene[id]);
        else
          HCAL_eneT.push_back(0.);
      }
    }
  }
  return HCAL_eneT;
}


/*the following function is used to define in-time energy for srd cells.
 */
ROOT::RVec<double> LYSOEneInTime(RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma, int nsigma) {

  //const int nsigma = 5.;

  ROOT::RVec<double> LYSO_eneT;
  for (int ix = 0; ix < 6; ix++) {
    auto deltaT = t0[ix] - tE[ix];
    if (fabs(deltaT) < nsigma * tSigma[ix])
      LYSO_eneT.push_back(ene[ix]);
    else
      LYSO_eneT.push_back(0.);
  }
  return LYSO_eneT;
}

ROOT::RVec<double> VetoEneInTime(RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma, int nsigma) {

  //const int nsigma = 5.;

  ROOT::RVec<double> VETO_eneT;
  for (int ix = 0; ix < 6; ix++) {
    auto deltaT = t0[ix] - tE[ix];
    if (fabs(deltaT) < nsigma * tSigma[ix])
      VETO_eneT.push_back(ene[ix]);
    else
      VETO_eneT.push_back(0.);
  }
  return VETO_eneT;
}

/* HCAL pedestal cut function *****************************/
/**********************************************************/

//Cut on noise per cell
ROOT::RVec<double> HcalEnePed(RVecD &ene) {
  const double HCNoise[36] = {0.47,0.79,0.51,0.83,0.94,0.48,0.61,0.64,0.85,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  ROOT::RVec<double> HCAL_eneP;
  for (int imod = 0; imod < 4; imod++) {
    for (int ix = 0; ix < 3; ix++) {
      for (int iy = 0; iy < 3; iy++) {
        auto id = hcalID(imod, ix, iy);
        if (ene[id] > HCNoise[id]) 
          HCAL_eneP.push_back(ene[id]); 
        else
          HCAL_eneP.push_back(0.); 
      }
    }
  }
  return HCAL_eneP;
}

// FILTER FUNCTIONS --------------------------------------------------
//Trigger time filter
bool filterTriggerTime(double &masterT) {
  bool ret = ((masterT >= 59.5) && (masterT <= 85.5));
  return ret;
}

//From Anna's analysis
bool filterInvertedLYSO(RVecD &LYSO_eneT) {
  int nLYSO=0;
  double totLYSO=0;

  const double thrModule=1E-3;
  const double thrTotMin=2E-3;
  const double thrTotMax=250E-3;
  
  for (int ii=0;ii<6;ii++){
    totLYSO+=LYSO_eneT[ii];
    if (LYSO_eneT[ii]>thrModule)
      nLYSO++;
  }

  //condition of Anna to select electrons:
  //Multiplicity > 2
  //Total energy btween 2 MeV and 250 MeV
  bool isElectron=((nLYSO>2)&&(totLYSO>thrTotMin)&&(totLYSO<thrTotMax));

  return (!isElectron);
  
}

bool filterInvertedVETO(RVecD &VETO_eneT){
  const double vetoCut01 = 0.0076;
  double vetoCut23   = 0.013; // select MIP
  double vetoCut23UP = 0.038; // select MIP
  const double vetoCut45 = 0.0076;

  double eVETO01(0), eVETO23(0), eVETO45(0);
  // VETO 0
  eVETO01 += VETO_eneT[0]; eVETO01 += VETO_eneT[1];
  // VETO 1
  eVETO23 += VETO_eneT[2]; eVETO23 += VETO_eneT[3];
  // VETO 2
  eVETO45 += VETO_eneT[4]; eVETO45 += VETO_eneT[5];

  return ((eVETO01 < vetoCut01) && (eVETO23 > vetoCut23) && (eVETO23 < vetoCut23UP) && (eVETO45 < vetoCut45));
}

//Function to produce histograms, to be called on the result of the "Filter" call
std::vector<ROOT::RDF::RResultPtr<TH1>> doHisto(ROOT::RDF::RNode df, int icut, string cutName, int runN) {

  std::vector<ROOT::RDF::RResultPtr<TH1>> ret;
  /*
  df = df.Define("SRDtot_T", "Sum(SRD_eneT)");
  df = df.Define("SRD0_T", "SRD_eneT[0]");
  df = df.Define("SRD1_T", "SRD_eneT[1]");
  df = df.Define("SRD2_T", "SRD_eneT[2]");

  //ECAL
  ret.push_back(df.Histo1D( { Form("ECAL_T3_%i_%i", icut, runN), "EC_T3", 2000u, 0., 200. }, "EC_T3"));
  ret.push_back(df.Histo1D( { Form("ECAL_T4_%i_%i", icut, runN), "EC_T4", 2000u, 0., 200. }, "EC_T4"));
  ret.push_back(df.Histo1D( { Form("ECAL_T5_%i_%i", icut, runN), "EC_T5", 2000u, 0., 200. }, "EC_T5"));

  //Veto energy
  //ret.push_back(df.Define("eveto", "VETO_enePMT[0]+VETO_enePMT[1]").Histo1D( { Form("Veto01_E_%i", icut), "Veto01_E", 400u, 0., 0.1 }, "eveto"));
  //ret.push_back(df.Define("eveto", "VETO_enePMT[2]+VETO_enePMT[3]").Histo1D( { Form("Veto23_E_%i", icut), "Veto23_E", 400u, 0., 0.1 }, "eveto"));
  //ret.push_back(df.Define("eveto", "VETO_enePMT[4]+VETO_enePMT[5]").Histo1D( { Form("Veto45_E_%i", icut), "Veto45_E", 400u, 0., 0.1 }, "eveto"));
  ret.push_back(df.Define("eveto", "VETO_enePMT_T[0]+VETO_enePMT_T[1]").Histo1D( { Form("Veto01_E_T_%i", icut), "Veto01_E_T", 400u, 0., 0.1 }, "eveto"));
  ret.push_back(df.Define("eveto", "VETO_enePMT_T[2]+VETO_enePMT_T[3]").Histo1D( { Form("Veto23_E_T_%i", icut), "Veto23_E_T", 400u, 0., 0.1 }, "eveto"));
  ret.push_back(df.Define("eveto", "VETO_enePMT_T[4]+VETO_enePMT_T[5]").Histo1D( { Form("Veto45_E_T_%i", icut), "Veto45_E_T", 400u, 0., 0.1 }, "eveto"));

  //SRD
  ret.push_back(df.Histo1D( { Form("SRDtot_T_%i", icut), "SRDtot_T", 400u, 0., 0.2 }, "SRDtot_T"));
  ret.push_back(df.Histo1D( { Form("SRD0_T_%i", icut), "SRD0_T", 400u, 0., 0.2 }, "SRD0_T"));
  ret.push_back(df.Histo1D( { Form("SRD1_T_%i", icut), "SRD1_T", 400u, 0., 0.2 }, "SRD1_T"));
  ret.push_back(df.Histo1D( { Form("SRD2_T_%i", icut), "SRD2_T", 400u, 0., 0.2 }, "SRD2_T"));
  
  //HCAL no time cut
  ret.push_back(df.Histo1D( { Form("HCAL012_%i", icut), "HCAL012", 2000u, 0., 200. }, "HCAL012"));
  ret.push_back(df.Histo1D( { Form("HCAL01_%i", icut), "HCAL01", 2000u, 0., 200. }, "HCAL01"));
  ret.push_back(df.Histo1D( { Form("HCAL0_%i", icut), "HCAL0", 2000u, 0., 200. }, "ehcal0"));
  ret.push_back(df.Histo1D( { Form("HCAL1_%i", icut), "HCAL1", 2000u, 0., 200. }, "ehcal1"));
  ret.push_back(df.Histo1D( { Form("HCAL2_%i", icut), "HCAL2", 2000u, 0., 200. }, "ehcal2"));
  ret.push_back(df.Histo1D( { Form("HCAL0_%i_%i", icut, runN), "HCAL0", 2000u, 0., 200. }, "ehcal0"));

  ret.push_back(df.Histo1D( { Form("HC011_%i_%i", icut, runN), "HC011", 2000u, 0., 200. }, "HC011"));
  ret.push_back(df.Histo1D( { Form("HC111_%i_%i", icut, runN), "HC111", 2000u, 0., 200. }, "HC111"));
  ret.push_back(df.Histo1D( { Form("HC211_%i_%i", icut, runN), "HC211", 2000u, 0., 200. }, "HC211"));

  //HCAL0 - OOT
  ret.push_back(df.Define("ehcal0_OOT", "ehcal0-HC0_T7").Histo1D( { Form("HCAL0_OOT_%i", icut), "HCAL0_OOT", 2000u, 0., 200. }, "ehcal0_OOT"));
 */ 
  //Time cut - for stability check (run by run)
  for (int sigma = 3; sigma <= 7; ++sigma) {
    ret.push_back(df.Histo1D( { Form("HCAL012_T%i_%i_%i", sigma, icut, runN), Form("HCAL012_T%i", sigma), 2000u, 0., 200. }, Form("HC012_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HCAL01_T%i_%i_%i", sigma, icut, runN), Form("HCAL01_T%i", sigma), 2000u, 0., 200. }, Form("HC01_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HCAL0_T%i_%i_%i", sigma, icut, runN), Form("HCAL0_T%i", sigma), 2000u, 0., 200. }, Form("HC0_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HCAL1_T%i_%i_%i", sigma, icut, runN), Form("HCAL1_T%i", sigma), 2000u, 0., 200. }, Form("HC1_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HCAL2_T%i_%i_%i", sigma, icut, runN), Form("HCAL2_T%i", sigma), 2000u, 0., 200. }, Form("HC2_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HC011_T%i_%i_%i", sigma, icut, runN), Form("HC011_T%i", sigma), 2000u, 0., 200. }, Form("HC011_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HC111_T%i_%i_%i", sigma, icut, runN), Form("HC111_T%i", sigma), 2000u, 0., 200. }, Form("HC111_T%i", sigma)));
    ret.push_back(df.Histo1D( { Form("HC211_T%i_%i_%i", sigma, icut, runN), Form("HC211_T%i", sigma), 2000u, 0., 200. }, Form("HC211_T%i", sigma)));
  }
  
  ret.push_back(df.Histo1D( { Form("HCAL0_P_%i_%i", icut, runN), "HCAL0_P", 2000u, 0., 200. }, "HC0_P"));
  ret.push_back(df.Histo1D( { Form("HCAL1_P_%i_%i", icut, runN), "HCAL1_P", 2000u, 0., 200. }, "HC1_P"));
  ret.push_back(df.Histo1D( { Form("HCAL2_P_%i_%i", icut, runN), "HCAL2_P", 2000u, 0., 200. }, "HC2_P"));
  ret.push_back(df.Histo1D( { Form("HCAL012_T7_P_%i_%i", icut, runN), "HCAL012_T7_P", 2000u, 0., 200. }, "HC012_T7_P"));
  ret.push_back(df.Histo1D( { Form("HCAL01_T7_P_%i_%i", icut, runN), "HCAL01_T7_P", 2000u, 0., 200. }, "HC01_T7_P"));
  ret.push_back(df.Histo1D( { Form("HCAL0_T7_P_%i_%i", icut, runN), "HCAL0_T7_P", 2000u, 0., 200. }, "HC0_T7_P"));
  ret.push_back(df.Histo1D( { Form("HCAL1_T7_P_%i_%i", icut, runN), "HCAL1_T7_P", 2000u, 0., 200. }, "HC1_T7_P"));
  ret.push_back(df.Histo1D( { Form("HCAL2_T7_P_%i_%i", icut, runN), "HCAL2_T7_P", 2000u, 0., 200. }, "HC2_T7_P"));


  return ret;
}

//Get all outputs
std::vector<ROOT::RDF::RResultPtr<TH1>> allH;
std::vector<ROOT::RDF::RResultPtr<ULong64_t>> allN;
std::vector<std::string> allCutName;

void postCutProcess(int idx, string cutName, ROOT::RDF::RNode df, int run) {
  allN.push_back(df.Count());
  allCutName.push_back(cutName);
  auto ret = doHisto(df, idx, cutName, run);
  for (auto h : ret) {
    allH.push_back(h);
  }
}

void ana_hadrons(string fIn = "", string fOut = "") {

  
  //ROOT::EnableImplicitMT();

  string fnameIn, fnameOut;
  if (fIn == "") {
    fnameIn = "./electron-2025A.root";
  } else {
    fnameIn = fIn;
  }

  if (fOut == "") {
    fnameOut = "out.root";
  } else {
    fnameOut = fOut;
  }

  //Open the data-frame
  ROOT::RDataFrame d("tout", fnameIn.c_str());
  
  // BLINDING
  auto dBlind = d.Filter("isSignalBox == 0");

  // PHYSICS TRIGGER or BEAM TRIGGER
  //auto dTrigger = d.Filter("trg == 3 || trg == 2");

  //Define columns
  auto d2 = dBlind.Define("EC", "eecal0+eecal1");
  d2 = d2.Define("HCAL01", "ehcal0+ehcal1");
  d2 = d2.Define("HCAL012", "ehcal0+ehcal1+ehcal2");
  d2 = d2.Define("HC011", "HCAL_ene[4]");
  d2 = d2.Define("HC111", "HCAL_ene[13]");
  d2 = d2.Define("HC211", "HCAL_ene[22]");
  //d2 = d2.Define("SRD_eneT",      SRDEneInTime,   { "SRD_ene", "SRD_t0", "SRD_tE", "SRD_tSigma" });
  d2 = d2.Define("periodN", GetPeriodFromRunN, { "runN" });
  auto center = ecalID(1,2,2); 
  //d2 = d2.Define("EC122_A", [&](const ROOT::VecOps::RVec<Double_t>& ECAL_A) {return ECAL_A[center];}, {"ECAL_A"});
  
  auto runNumber = d2.Take<unsigned int>("runN");
  std::vector<unsigned int> runNumberVec = runNumber.GetValue();

  //Get runN int
  unsigned int runN = 0;  // Default value
  if (!runNumberVec.empty()) {
    runN = runNumberVec[0];
  } else {
    std::cerr << "Warning: runN could not be determined!\n";
  }


  //In-time energy selection
  //========================
  for (int sigma = 3; sigma <= 7; ++sigma) {
    d2 = d2.Define(Form("ECAL_eneT%d", sigma), 
                   [sigma](RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma) { 
                       return EcalEneInTime(ene, t0, tE, tSigma, sigma);
                   }, 
                   { "ECAL_ene", "ECAL_t0", "ECAL_tE", "ECAL_tSigma" });

    d2 = d2.Define(Form("HCAL_eneT%d", sigma), 
                   [sigma](RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma) { 
                       return HcalEneInTime(ene, t0, tE, tSigma, sigma);
                   }, 
                   { "HCAL_ene", "HCAL_t0", "HCAL_tE", "HCAL_tSigma" });
 
   
    d2 = d2.Define(Form("EC_T%d", sigma), Form("Sum(ECAL_eneT%d)", sigma));
    d2 = d2.Define(Form("HC012_T%d", sigma), Form("Sum(Take(HCAL_eneT%d,27))", sigma));
    d2 = d2.Define(Form("HC01_T%d", sigma), Form("Sum(Take(HCAL_eneT%d,18))", sigma));
    d2 = d2.Define(Form("HC0_T%d", sigma), Form("Sum(Take(HCAL_eneT%d,9))", sigma));
    d2 = d2.Define(Form("HC1_T%d", sigma), Form("Sum(Take(HCAL_eneT%d,18))-Sum(Take(HCAL_eneT%d,9))", sigma, sigma));
    d2 = d2.Define(Form("HC2_T%d", sigma), Form("Sum(Take(HCAL_eneT%d,27))-Sum(Take(HCAL_eneT%d,18))", sigma, sigma));
    d2 = d2.Define(Form("HC3_T%d", sigma), Form("Sum(Take(HCAL_eneT%d,-9))", sigma));
    
    d2 = d2.Define(Form("HC011_T%d", sigma), Form("HCAL_eneT%d[4]", sigma));
    d2 = d2.Define(Form("HC111_T%d", sigma), Form("HCAL_eneT%d[13]", sigma));
    d2 = d2.Define(Form("HC211_T%d", sigma), Form("HCAL_eneT%d[22]", sigma));
  }

  for (int sigma = 2; sigma <= 7; ++sigma) {
    d2 = d2.Define(Form("VETO_enePMT_T%d", sigma), 
                   [sigma](RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma) { 
                       return VetoEneInTime(ene, t0, tE, tSigma, sigma);
                   }, 
                   { "VETO_enePMT", "VETO_t0", "VETO_tE", "VETO_tSigma" });
    d2 = d2.Define(Form("LYSO_ene_T%d", sigma), 
                   [sigma](RVecD &ene, RVecD &t0, RVecD &tE, RVecD &tSigma) { 
                       return LYSOEneInTime(ene, t0, tE, tSigma, sigma);
                   }, 
                   { "LYSO_ene", "LYSO_t0", "LYSO_tE", "LYSO_tSigma" });
 }
  
  d2 = d2.Define("HCAL_eneP", HcalEnePed, { "HCAL_ene"});
  d2 = d2.Define("HCAL_eneT7_P", HcalEnePed, { "HCAL_eneT7"});
  
  d2 = d2.Define("HC0_P", "Sum(Take(HCAL_eneP,9))");
  d2 = d2.Define("HC1_P","Sum(Take(HCAL_eneP,18))-Sum(Take(HCAL_eneP,9))");
  d2 = d2.Define("HC2_P","Sum(Take(HCAL_eneP,27))-Sum(Take(HCAL_eneP,18))");
  d2 = d2.Define("HC0_T7_P","Sum(Take(HCAL_eneT7_P,9))");
  d2 = d2.Define("HC01_T7_P", "Sum(Take(HCAL_eneT7_P,18))");
  d2 = d2.Define("HC012_T7_P","Sum(Take(HCAL_eneT7_P,27))");
  d2 = d2.Define("HC1_T7_P","Sum(Take(HCAL_eneT7_P,18))-Sum(Take(HCAL_eneT7_P,9))");
  d2 = d2.Define("HC2_T7_P", "Sum(Take(HCAL_eneT7_P,27))-Sum(Take(HCAL_eneT7_P,18))");
  d2 = d2.Define("HC011_P","HCAL_eneP[4]");
  d2 = d2.Define("HC111_P","HCAL_eneP[13]");
  d2 = d2.Define("HC211_P","HCAL_eneP[22]");

  
  //d2 = d2.Define("VETO_enePMT_T", VetoEneInTime,  { "VETO_enePMT", "VETO_t0", "VETO_tE", "VETO_tSigma" });

  //Cutflow
  //=======

  //postCutProcess(0, "AllEvents",d2, runN);
  
  auto dF1 = d2.Filter(filterTriggerTime, { "masterT" });
  //postCutProcess(1, "Trigger Time", dF1, runN);

  auto dF2 = dF1.Filter("EC_T5 < 5.");
  //postCutProcess(2, "EC_T5 < 5 ", dF2, runN);

  auto dF3 = dF2.Filter(filterInvertedLYSO, { "LYSO_ene_T5" }); //Anna was using a 5-sigma cut for LYSO
  //postCutProcess(3, "Inverted SRD ", dF3, runN);

  auto dF4 = dF3.Filter(filterInvertedVETO, { "VETO_enePMT_T4" });
  //postCutProcess(4, "MIP Veto cut", dF4, runN);

  auto dF5 = dF4.Filter("HC012_T7 > 60.");
  postCutProcess(5, "HC012_T7 > 60 ", dF5, runN);

  //Check HCAL0
  auto dF6 = dF5.Filter("HC01_T7>60.");
  postCutProcess(6, "HC01_T7 > 60 ", dF6, runN);
  
  auto dF7 = dF6.Filter("HC0_T7>60.");
  postCutProcess(7, "HC0_T7 > 60 ", dF7, runN);

  //Check HCAL1
  auto dF8 = dF5.Filter("HC0_T7<5.");
  //postCutProcess(8, "HC0_T7 < 5 on dF5", dF8, runN);
  
  auto dF9 = dF8.Filter("HC1_T7>60.");
  postCutProcess(9, "HC1_T7 > 60", dF9, runN);
  
  //Check HCAL2
  auto dF10 = dF5.Filter("HC01_T7<10.");
  //postCutProcess(10, "HC01_T7 < 10 on dF5", dF10, runN);

  auto dF11 = dF10.Filter("HC2_T7>60.");
  //postCutProcess(11, "HC2_T7 > 60", dF11, runN);

  /******************************************************/
 
  //Selection muons
  auto dF12 = dF2.Filter("Sum(Take(HCAL_ene,-9))<3.5");
  //postCutProcess(12, "HC3 cut", dF12, runN);
  
  //MUON in HC2-HC1 - strict
  auto dF13 = dF12.Filter("HC211_T7>2. && HC211_T7<6. && HC111_T7>2. && HC111_T7<6.");
  //postCutProcess(13, "MIP in HC1 & HC2 - strict", dF13, runN);
 
  //MUON in HC2-HC1 
  auto dF14 = dF12.Filter("HC211>1. && HC211<6. && HC111>1. && HC111<6.");
  postCutProcess(14, "MIP in HC1 & HC2", dF14, runN);

  //MUON in HC2-HC0
  auto dF15 = dF12.Filter("HC011>1. && HC011<6. && HC211>1. && HC211<6. ");
  postCutProcess(15, "MIP in HC0 & HC2", dF15, runN);

  //MUON in HC0-HC1
  auto dF16 = dF12.Filter("HC011>1. && HC011<6. && HC111>1. && HC111<6.");
  postCutProcess(16, "MIP in HC0 & HC1", dF16, runN);
  
  //MUON in HC2-HC1 T7
  auto dF17 = dF12.Filter("HC211_T7>1. && HC211_T7<6. && HC111_T7>1. && HC111_T7<6.");
  //postCutProcess(17, "MIP in HC1 & HC2", dF17, runN);

  auto dF18 = dF17.Filter(filterInvertedLYSO, { "LYSO_ene_T5" }); //Anna was using a 5-sigma cut for LYSO
  //postCutProcess(18, "Inverted SRD ", dF18, runN);

  /*
  //Print
  cout << "TotalEvents: " << (allN[0].GetValue()) << endl;

  cout << endl;
  cout << "Ordered cutflow" << endl;
  printf("Cut:|CutName\t\t\t   |NumEvents\t\t|Eff\t\t\t|EffRel\t\t\n");
  for (int ii = 1; ii < allN.size(); ii++) {
    double eff = (allN[ii - 1].GetValue() - allN[ii].GetValue()) / (1. * allN[0].GetValue());
    double effRel = (allN[ii - 1].GetValue() - allN[ii].GetValue()) / (1. * allN[ii - 1].GetValue());

    printf("%.2i: |%-30s|%llu\t\t|%f\t\t|%f\n", ii, allCutName[ii].c_str(), allN[ii].GetValue(), eff, effRel);
  }

  cout << endl;
  cout << "Final eff: " << allN.back().GetValue() / (1. * allN[0].GetValue()) << std::endl;
  */
  //Write to the output
  TFile *fout = new TFile(fnameOut.c_str(), "recreate");
  fout -> mkdir("Cutflow");
  fout -> cd("Cutflow");
  for (auto h : allH) {
    h->Write();
  }
  fout->Close();

  //Save a snapshot 
  ROOT::RDF::RSnapshotOptions opts;
  opts.fMode = "update";
  dF18.Snapshot("muons",fnameOut.c_str(),"",opts);

}