trimvar.C

Go to the documentation of this file.

//............................................
// This script will be usefull for trimming 
// caf files to acess some variable for 
// further Training and Testing using TMVA
//
// Author: Biswaranjan Behera 
// email:  bbehera@fnal.gov
//...........................................
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/SAMQuerySource.h"
#include "CAFAna/Core/SAMProjectSource.h"
#include "CAFAna/Core/WildcardSource.h"
#include "CAFAna/Core/IFileSource.h"
#include "CAFAna/Core/FileListSource.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/OscillatableSpectrum.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/OscCurve.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/EventList.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"

#include "StandardRecord/SRSpillTruthBranch.h"

#include "TLorentzVector.h"
#include <cmath>
#include <iostream>
#include <functional>
#include <list>
#include <memory>
#include <set>
#include <string>
#include <vector>

//From Root
#include "TStyle.h"
#include "TFile.h"
#include "TROOT.h"
#include "TTree.h"
#include "TVector3.h"
#include "TH1.h"

#include "StandardRecord/StandardRecord.h"

using namespace ana;

void trimvar(const std::string outputFilename, const std::string setname)
{

  // Nominal RHC MC defname
  const std::string trainingSet = "muonid_training_sample";
  const std::string testingSet  = "muonid_testing_sample";

  // const std::string samQuery = "defname: " + defname;
  std::string samQuery = ""; //"defname: " + defname;
  if (setname == "train"){
    cout << "Making Training Trim file from dataset: " << trainingSet << endl;
    samQuery = "dataset_def_name_newest_snapshot " + trainingSet;
  }
  else if (setname == "test"){
    cout << "Making Testing Trim File from dataset: " << testingSet << endl;
    samQuery = "dataset_def_name_newest_snapshot " + testingSet;
  }
  else{
    cerr << "Setname should be either \"train\" or \"test\". Exiting." << endl;
    return;
  }

  TFile outputFile(outputFilename.c_str(), "RECREATE");
  
  TTree *sigTree = new TTree("sigTree", "Signal Events");
  TTree *bkgTree = new TTree("bkgTree", "Background Events");
  
  //...............................................................................
  
  // Six variable for MuonID
  int fismuon = -5.0;
  int ismuon  = -5.0;
  
  float  fScatLL           = -999.0;
  float  fDedxLL           = -999.0;
  float  fAvededxlast10cm  = -999.0;
  float  fAvededxlast20cm  = -999.0;
  float  fAvededxlast30cm  = -999.0;
  float  fAvededxlast40cm  = -999.0;
  
  //-------------------------------------------------------------------
  
  // Define the branches on the Tree
  
  // Six Muon PID variable
  sigTree->Branch("DedxLL",            &fDedxLL,                   "DedxLL/F");
  sigTree->Branch("ScatLL",            &fScatLL,                   "ScatLL/F");
  sigTree->Branch("Avededxlast10cm",   &fAvededxlast10cm,  "Avededxlast10cm/F");
  sigTree->Branch("Avededxlast20cm",   &fAvededxlast20cm,  "Avededxlast20cm/F");
  sigTree->Branch("Avededxlast30cm",   &fAvededxlast30cm,  "Avededxlast30cm/F");
  sigTree->Branch("Avededxlast40cm",   &fAvededxlast40cm,  "Avededxlast40cm/F");

  bkgTree->Branch("DedxLL",            &fDedxLL,                   "DedxLL/F");
  bkgTree->Branch("ScatLL",            &fScatLL,                   "ScatLL/F");
  bkgTree->Branch("Avededxlast10cm",   &fAvededxlast10cm,  "Avededxlast10cm/F");
  bkgTree->Branch("Avededxlast20cm",   &fAvededxlast20cm,  "Avededxlast20cm/F");
  bkgTree->Branch("Avededxlast30cm",   &fAvededxlast30cm,  "Avededxlast30cm/F");
  bkgTree->Branch("Avededxlast40cm",   &fAvededxlast40cm,  "Avededxlast40cm/F");
  
  //---------------------------------------------------------------------
  
  //  Begin looping over the selected files from the SAM Database.
  
  int stride = -1;
  int offset = -1;
  int limit  = -1;

  if(getenv("CAFANA_STRIDE")){
    stride = atoi(getenv("CAFANA_STRIDE"));
    if(stride > 1 && getenv("CAFANA_OFFSET")){
      offset = atoi(getenv("CAFANA_OFFSET"));
    }
  }


  SAMQuerySource* filesrc = new SAMQuerySource(samQuery, stride, offset, limit);
  
  int Nfiles = filesrc->NFiles();

  Progress* prog = 0;
  
  int fileIdx = -1;
  bool foundFile = true;
  TFile* file1 = filesrc->GetNextFile();
  if(!file1) foundFile = false; // out of files
  assert(!file1->IsZombie());
  
  // Swapped from while loop to do-while, so we don't lose the first file
  do {
      //Start with some file management.
      ++fileIdx;

      if(!prog) prog = new Progress(TString::Format("Filling signal and background trees from %d files", Nfiles).Data());

      if(Nfiles > 1 && prog) prog->SetProgress((fileIdx+1.)/Nfiles);

      // Define the recTree object
      TTree *recTree = (TTree*)file1->Get("recTree");

      TTree *spillTree = (TTree*)file1->Get("spillTree");

      // Every time a GetEntry(i) is done, the data in the tree is copied to this object
      caf::StandardRecord* sr = 0;
      recTree->SetBranchAddress("rec", &sr);

      caf::SRSpill* sp = 0;
      spillTree->SetBranchAddress("spill", &sp);

      //--------------------------------------------------------
      // Standard Spill Cuts
      //--------------------------------------------------------
      
      bool ismc = false;
      if(sp->ismc) ismc = true;

      //cosmic trigger is 02,
      //skip beam cuts
      bool istrigger =false;
      if(sp->trigger == 2) istrigger = true;

      // docdb 12474 slide 11
      bool goodspill = false;
      if(sp->spilltimesec==0 && sp->deltaspilltimensec==0 && sp->widthx==0)
        goodspill = bool(sp->isgoodspill); // decafs have spill info removed 
                                      // and logic applied in decaf step 
                                      // so isgoodspill is equivalent to 
                                      // tightbeamqualitycuts

      if(std::abs(sp->deltaspilltimensec) > 0.5e9) continue;

      if(sp->spillpot < 2e12) continue;
      
      if(sp->hornI < -202 || sp->hornI > -198) continue;
      
      if(sp->posx < -2.00 || sp->posx > +2.00) continue;
      
      if(sp->posy < -2.00 || sp->posy > +2.00) continue;
      
      if(sp->widthx < 0.57 || sp->widthx > 1.58) continue;
      
      if(sp->widthy < 0.57 ||  sp->widthy > 1.58) continue;

      bool gdspill = false;
      if(sp->dcmedgematchfrac==0 && sp->fracdcm3hits==0 && sp->nmissingdcmslg==0)
        gdspill =  bool(sp->isgoodspill);
      
      // fracdcm3hits > 0.7 for releases before 14-11-25,
      // ND specific cut in cases where lights were left on
      // nmissingdcms counts 0 hit dcms. No livegeo.
      if(sp->det == caf::kNEARDET  &&
      (sp->fracdcm3hits > 0.45  ||
       sp->nmissingdcms > 0) )
        continue;
      
      // cut out events which have the isEventIncomplete DAQ header flag
      // set.
      if ( sp->eventincomplete ) continue;

      // nmissingdcmslg counts 0 hit dcms among those
      // that livegeo think are live.
      if( sp->det == caf::kFARDET &&
       sp->nmissingdcmslg > 0 )
        continue;
      

      // FD specific cut to detect out of sync dcms.
      // Only works if cosmics are overlaid. so we
      // don't apply it on MC by default, although
      // I assume you could safely apply it on cry.
      if( sp->det == caf::kFARDET &&
       !sp->ismc &&
       sp->dcmedgematchfrac <= 0.2)
        continue;


      //--------------------------------------------------------

      //Looping over entries
      Int_t numberOfEntries = recTree->GetEntriesFast();

      int i = 0;
      //Begin filling the branches with variables and cuts applied.
      for (Int_t event = 0; event < numberOfEntries; ++event) {
        recTree->GetEntry(event);
    
        if(sr->vtx.nelastic == 0)      continue;
        if(sr->mc.nnu < 1)             continue;
        if(sr->trk.kalman.ntracks < 1) continue;

        int nkals = sr->trk.kalman.ntracks;
    
        //bool isquality   = false;
    
        // bool isContained = false;
        bool isquality   = false;
        bool iscontprong  = true;
        // bool isfiducial  = false;
    
        for( unsigned int i = 0; i < sr->vtx.elastic[0].fuzzyk.nshwlid; ++i ) {
          TVector3 start = sr->vtx.elastic[0].fuzzyk.png[i].shwlid.start;
          TVector3 stop  = sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop;
          if( std::min( start.X(), stop.X() ) < -180.0 ||
              std::max( start.X(), stop.X() ) >  180.0 ||
              std::min( start.Y(), stop.Y() ) < -180.0 ||
              std::max( start.Y(), stop.Y() ) >  180.0 ||
              std::min( start.Z(), stop.Z() ) <   20.0 ||
              std::max( start.Z(), stop.Z() ) > 1525.0 ){
                iscontprong = false;
          }
        }


        if(sr->trk.kalman.ntracks > 0  &&
          sr->slc.nhit           > 20 &&
          sr->slc.ncontplanes    > 4)
            isquality    = true;
        else
            isquality    = false;    

        if (isquality && iscontprong){
          i++;  

          if( nkals < 1 ) continue;
        
          // loop over track
          for (int itrk = 0; itrk < nkals; ++itrk){
            
            if( sr->trk.kalman.tracks[itrk].rempid < 0 ) continue;        
            
            fDedxLL   = sr->trk.kalman.tracks[itrk].dedxllh2017;
            
            fScatLL   = sr->trk.kalman.tracks[itrk].scatllh2017;        
            
            if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
            if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
            if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
            if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;

            if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30)
              fAvededxlast10cm = 30;
            else
              fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
            
            if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm > 30)
              fAvededxlast20cm = 30;
            else
              fAvededxlast20cm = sr->trk.kalman.tracks[itrk].avedEdxlast20cm;
            
            if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm > 30)
              fAvededxlast30cm = 30;
            else
              fAvededxlast30cm = sr->trk.kalman.tracks[itrk].avedEdxlast30cm;
            
            if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30)
              fAvededxlast40cm = 30;
            else
              fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
            
            // If true muon, fill the signal tree. Otherwise fill in into the bkgd tree.
            if(abs(sr->trk.kalman.tracks[itrk].truth.pdg) == 13)
              sigTree->Fill();
            else
              bkgTree->Fill();

          } // end of trk loop
          
        }
        else continue;
        // end of pre selected if statement   
    
    
      } // end of event loop  
     
  } while( (file1 =(TFile*)filesrc->GetNextFile()) );
  // end of while loop of genie files

  if(prog){
    prog->Done();
    delete prog;
  }
  
  outputFile.Write();
  outputFile.Close();
  
  std::cout << " Congratulation! You did it succesfully..." << std::endl;
  
} // end of void