MadTestAnalysis.cxx

Go to the documentation of this file.

#ifndef madtestanalysis_cxx
#define madtestanalysis_cxx
#include <iostream>
#include "Mad/MadTestAnalysis.h"
#include "TClonesArray.h"
#include "TMath.h"
#include "Validity/VldContext.h"
#include "Conventions/Detector.h"
#include "Registry/Registry.h"
#include "Mad/MadAnalysis.h"
#include "MCReweight/MCReweight.h"
#include "MCReweight/GnumiInterface.h"
#include "MCReweight/NuParent.h"
#include "Mad/SpillInfo.h"
#include "Mad/SpillInfoBlock.h"
#include "BeamDataUtil/BDSpillAccessor.h"
#include "BeamDataUtil/BeamMonSpill.h"
#include "SpillTiming/SpillTimeFinder.h"
#include "Riostream.h"
#include "Mad/ScanList.h"

using namespace std;

//********************************************************
MadTestAnalysis::MadTestAnalysis(TChain *chainSR,TChain *chainMC,
                                 TChain *chainTH,TChain *chainEM)
{

  if(!chainSR) {
    record = 0;
    strecord = 0;
    emrecord = 0;
    mcrecord = 0;
    threcord = 0;
    Clear();
    whichSource = -1;
    isMC = true;
    isTH = true;
    isEM = true;
    Nentries = -1;
    return;
  }
  
  InitChain(chainSR,chainMC,chainTH,chainEM);
  whichSource = 0;
  fLikeliFile = NULL;
  for(int i=0;i<6;i++) fLikeliHist[i] = NULL;

}

//********************************************************
MadTestAnalysis::MadTestAnalysis(JobC *j,string path,int entries)
{
  record = 0;
  strecord = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  Clear();
  isMC = true;
  isTH = true;
  isEM = true;
  Nentries = entries;
  whichSource = 1;
  jcPath = path;
  fJC = j;
  fLikeliFile = NULL;
  for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
}

//********************************************************
MadTestAnalysis::~MadTestAnalysis()
{
  if(fLikeliFile) fLikeliFile->Close();
}

//*********************************************************
Bool_t MadTestAnalysis::PassTruthSignal(Int_t mcevent)
{
  if(!LoadTruth(mcevent)) return false;
  if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true;
  return false;
}

//*********************************************************
Bool_t MadTestAnalysis::PassAnalysisCuts(Int_t event)
{
  if(!LoadEvent(event)) return false;
  //  if(PassBasicCuts(event)) return true;
  Float_t pidcut=-0.4;
  if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;}
  
  if(!(PID(event,0)>pidcut)) return false; 
  return true;
}

//******************************************
Bool_t MadTestAnalysis::PassBasicCuts(Int_t event)
{ 
  if(!LoadEvent(event)) return false;
   if(ntpEvent->ntrack==0) return false;
   Int_t track = -1;
   LoadLargestTrackFromEvent(event,track);  
   if(!PassTrackCuts(track)) return false;  
   if(!IsFidVtxEvt(ntpEvent->index)) return false;  
//  if(!IsFidVtx(track)) return false;
  return true;
}
Bool_t MadTestAnalysis::PassFid(Int_t event)
{ 
  if(!LoadEvent(event)) return false;
  if(!IsFidVtxEvt(ntpEvent->index)) return false;  
  return true;
}
//******************************************
// ANN VARIABLES
//
AnnInputBlock MadTestAnalysis::AnnVar(Int_t event) 
{
AnnInputBlock anninput;
  // Initialize
anninput.aTotrk       =0.;
anninput.aTotstp      =0.;
anninput.aTotph       =0.;
anninput.aTotlen      =0.;
anninput.aPhperpl     =0.;
anninput.aPhperstp    =0.;


anninput.aTrkpass     =0.;
anninput.aTrkph       =0.;
anninput.aTrklen      =0.;
anninput.aTrkphperpl  =0.;
anninput.aTrkphper    =0.;
anninput.aTrkplu      =0.;
anninput.aTrkplv      =0.;
anninput.aTrkstp      =0.;

anninput.aShwph       =0.;
anninput.aShwstp      =0.;
anninput.aShwdig      =0.;
anninput.aShwpl       =0.;
anninput.aShwphper    =0.;
anninput.aShwphperpl  =0.;
anninput.aShwphperdig =0.;
anninput.aShwphperstp =0.;
anninput.aShwplu      =0.;
anninput.aShwplv      =0.;
anninput.aPh3         =0.;
anninput.aPh6         =0.;
anninput.aPhlast      =0.;
anninput.aPhcommon    =0.;
//  
  
// Calculate variables needed for ANN (and a few more)
  
  LoadEvent(event) ;  
  int track_index   = -1;
  LoadLargestTrackFromEvent(event,track_index);
  int shower_index  = -1;
  LoadLargestShowerFromEvent(event,shower_index);
  LoadTrack(track_index);
  LoadShower(shower_index);
 
anninput.aTotrk       =ntpEvent->ntrack;
anninput.aTotstp      =ntpEvent->nstrip;
anninput.aTotph       =ntpEvent->ph.sigcor;
anninput.aTotlen      =ntpEvent->plane.end-ntpEvent->plane.beg+1; // prepei na to alla3w
anninput.aPhperpl     =ntpEvent->ph.sigcor/ntpEvent->plane.n;
anninput.aPhperstp    =ntpEvent->ph.sigcor/ntpEvent->nstrip;

if(track_index!=-1) {
anninput.aTrkpass     =ntpTrack->fit.pass;
anninput.aTrkph       =ntpTrack->ph.sigcor;
anninput.aTrklen      =ntpTrack->plane.ntrklike;
if(ntpTrack->plane.ntrklike>0)anninput.aTrkphperpl  =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
if(ntpEvent->ph.sigcor>0)     anninput.aTrkphper    =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
anninput.aTrkplu      =ntpTrack->plane.nu;
anninput.aTrkplv      =ntpTrack->plane.nv;
anninput.aTrkstp      =ntpTrack->nstrip;
}
if(shower_index!=-1) {
anninput.aShwph       =ntpShower->ph.sigcor;
anninput.aShwstp      =ntpShower->nstrip;
anninput.aShwdig      =ntpShower->ndigit;
anninput.aShwpl       =ntpShower->plane.n;
anninput.aShwphper    =ntpShower->ph.sigcor/ntpEvent->ph.sigcor;
anninput.aShwphperpl  =ntpShower->ph.sigcor/ntpShower->plane.n;
anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit;
anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip;
anninput.aShwplu      =ntpShower->plane.nu;
anninput.aShwplv      =ntpShower->plane.nv; 
}

      Int_t ssind,ssplind;
      Double_t ssphtot;
      Bool_t foundshw,foundtrk;
      Int_t planes=ntpEvent->plane.beg;
     
      Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
      ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;
// loop over strips to compute ph3 ph6 phlast phcommon
     for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
       Int_t stp_index=((ntpEvent->stp)[evss]);
       LoadStrip(stp_index);    
       ssind=ntpStrip->strip;
       ssplind=ntpStrip->plane;
       ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;       
         
       foundshw=false;
       foundtrk=false;
   
       Double_t phstrips=0;
       Double_t phstript=0;       
       // shower strips
       Int_t sshwind,sshwplind;
       Double_t sshwphtot;
       
       if(shower_index!=-1) {
       for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
        Int_t stp_indexs=((ntpShower->stp)[jj]);
        LoadStrip(stp_indexs);  
        
        if(!foundshw){
          sshwind=ntpStrip->strip;
          sshwplind=ntpStrip->plane;
          sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;             
            if(sshwind==ssind && sshwplind==ssplind) {
             foundshw=true;
             phstrips=sshwphtot;
           }
         }
       }
      }      
       // tracks strips
       Int_t strkind,strkplind;
       Double_t strkphtot;
      if(track_index!=-1) {
       for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){
        Int_t stp_indext=((ntpTrack->stp)[jj]);
        LoadStrip(stp_indext);  
         if(!foundtrk){
          strkind=ntpStrip->strip;
          strkplind=ntpStrip->plane;
          strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;   
           if(strkind==ssind && strkplind==ssplind) {
            foundtrk=true;
            phstript=strkphtot;
           }
         }
       }
      }
       
       if(foundshw && foundtrk) {
        hitcommon=hitcommon+1;
        phcommon=phcommon+phstrips+phstript;
       }  
       if(!foundshw && ! foundtrk) {
         hitnowere=hitnowere+1; 
         phnowere=phnowere+ssphtot; 
        }                
       if(ssplind>=planes && ssplind<=(planes+3)){
        ph3=ph3+ssphtot;
       }
       else if(ssplind>(planes+3) && ssplind<=(planes+6)){
        ph6=ph6+ssphtot; 
       }
       else {
        phlast=phlast+ssphtot;
       }     
     } // end of looping over event strips....    

//
anninput.aPh3       =ph3;
anninput.aPh6       =ph6;
anninput.aPhlast    =phlast;
anninput.aPhcommon  =phcommon; 
return anninput;
}
//******************************************
Double_t MadTestAnalysis::Sigmoid(Double_t x)
{
  double sig;
      if(x>37.){
         sig = 1.;
      }
      else if(x<-37.){
         sig = 0.;
      }
      else {
         sig = 1./(1.+exp(-x));
      }
      return sig;
}
 
//**************************************************
Double_t MadTestAnalysis::GetAnnPid(AnnInputBlock anninput,Int_t det, Int_t tar,Int_t fid, Int_t prior)
{
   // DET IS THE DETECTOR 1 = NEAR 2 = FAR
   // TAR IS THE ENERGY   1 = LE 2 = PME 3 = PHE
   // FID IS THE FIDUCIAL REGION FOR FAR , 2 = DAVIDS  1 = MINE (MORE STRICT)
   // PRIOR IS THE FAR TRAINED METHOD    , 1 = NO OSCILLATIONS , 2 = DM2 0.0025 SINTHETA2 = 0.95 
   //                                      3 = DM2 0.002 SINTHETA2 = 0.95
      
      int    inneuron  = 13;
      int    hidneuron = 15;
      double var;
     
      double out[15]; // input neurons
      double rin[15]; // first hidden layer
     
      double weight1[13][15];// weights first layer
      double constant1[15];   // constant first layer
      
      double weighto[15];// weights second  (output) layer
      double constanto[1];   // constant second (output) layer  
      double prob;
      
// Initialize    

       for(Int_t i=0;i<hidneuron;i++) {
         out[i]=0.;
         rin[i]=0.; 
         constant1[i]=0.;        
       }
            
       for(Int_t i=0;i<inneuron;i++) {
        for(Int_t j=0;j<hidneuron;j++){
          weight1[i][j]=0.;
        }
       }

        constanto[0]=0.;
        for(Int_t i=0;i<hidneuron;i++){
          weighto[i]=0.;
        }
      Int_t all  = inneuron*hidneuron+hidneuron+hidneuron+1;
      Int_t all1 = inneuron*hidneuron+hidneuron;
      Int_t n1  =0;
      Int_t nw1 =0;
      Int_t no  =0;
      Int_t nwo =0;
      ifstream weightfile;
      
//   INPUT VARIABLES
      if(det==2){  
       if(prior==1){ 
        if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farnooscilnew12.dat");    
        if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farnooscilnewfid12.dat");        
       }
       if(prior==2){ 
        if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc1new12.dat");    
        if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc1newfid12.dat");
       }
       if(prior==3){ 
        if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc2new12.dat");    
        if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc2newfid12.dat");
       }
      out[0]  = anninput.aTrkphperpl/2000.;
      out[1]  = anninput.aTotrk*anninput.aTrkpass;
      out[2]  = (anninput.aTotstp/100.);
      out[3]  = (anninput.aTotph/50000.);   
      out[4]  = (anninput.aTotlen/40);
      out[5]  = (anninput.aPhperpl/5000.);      
      out[6]  = (anninput.aPhperstp/1000.);
      out[10] = ((anninput.aTrkplv-anninput.aShwplv)+20.)/40.;
      out[11] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
      }
      else if(det==1){
       if(tar==1)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle16.dat");    
       if(tar==2)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme16.dat");    
       if(tar==3)  weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe16.dat");    
       out[0]  = anninput.aTrkphperpl/5000.;
       out[1]  = anninput.aTotrk*anninput.aTrkpass;
       out[2]  = (anninput.aTotstp/100.);
       out[3]  = (anninput.aTotph/100000.);   
       out[4]  = (anninput.aTotlen/40);
       out[5]  = (anninput.aPhperpl/10000.);    
       out[6]  = (anninput.aPhperstp/2000.);  
       out[10] = ((anninput.aTrkplv-anninput.aShwplv)+10.)/20.;
       out[11] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
      }
      out[7]  = (anninput.aPh3/(anninput.aTotph+1.));
      out[8]  = (anninput.aPh6/(anninput.aTotph+1.));
      out[9]  = (anninput.aPhlast/(anninput.aTotph+1.));
      out[12] = (anninput.aShwphperdig/800.);
      
// INPUT FILE
        
     for(Int_t k=0;k<all;k++){  
       weightfile >> var ;
       if(k<all1){
         if(k%(inneuron+1)==0){
          nw1=0;
          constant1[n1]=var;      
          n1=n1+1;         
         }
         else {
           weight1[nw1][n1-1]=var;         
           nw1=nw1+1;      
         }
        }
        else if(k>=all1){
          if((k-all1)%(hidneuron+1)==0){
          nwo=0;
          constanto[no]=var;      
          no=no+1;        
         }
         else {
           weighto[nwo]=var;       
           nwo=nwo+1;      
         }
        }          
       }
       
     weightfile.close();  
// end of reading the weights
   
       // first layer         
       for(Int_t i=0;i<hidneuron;i++){
        rin[i]=constant1[i];
        for(Int_t j=0;j<inneuron;j++){
         rin[i]=rin[i]+weight1[j][i]*out[j];
        }
       }
       // output                        
        for(Int_t i=0;i<hidneuron;i++){  
         out[i]=Sigmoid(rin[i]);         
        }               
        
        prob=constanto[0];      
        for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i];
        
                        
        if(anninput.aTotlen>40) prob=0.;
                
 return prob;
}
//*************************************************
//SP likelihood analysis
Float_t MadTestAnalysis::PID(Int_t event, Int_t method)
{


  if(!LoadEvent(event)) return -10;
  Int_t track = -1;
  if(!LoadLargestTrackFromEvent(event,track)) return -10;
  if(method!=0) return -10;
  
  Float_t ProbNC = 1.;
  Float_t ProbCC = 1.;
  Float_t PidCC = 0.;
  

  Float_t var1=eventSummary->plane.n;
  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
    var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
  }
  Float_t var2=0;
  Float_t var3=0;


  Float_t phtrack=ntpTrack->ph.sigcor;
  Float_t phevent=ntpEvent->ph.sigcor;

  if (phevent>0) {var2=phtrack/phevent;}

  if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor)
                              /float(ntpTrack->plane.n);


  Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
  Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
  Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));

  if (prob1==0) {prob1=0.0001;}
  if (prob2==0) {prob2=0.0001;}
  if (prob3==0) {prob3=0.0001;}
 
  ProbCC=prob1*prob2*prob3;

  prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
  prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
  prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
 
  if (prob1==0) {prob1=0.0001;}
  if (prob2==0) {prob2=0.0001;}
  if (prob3==0) {prob3=0.0001;}
 
  ProbNC=prob1*prob2*prob3;

  PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));

  return PidCC;
}

//********************************************************
Float_t MadTestAnalysis::RecoAnalysisEnergy(Int_t event){

  if(!LoadEvent(event)) return 0;
  int largestTrack = -1;
  LoadLargestTrackFromEvent(event,largestTrack);
  float emu = RecoMuEnergy(0,largestTrack);
  float ehad = RecoShwEnergy(event);
  float ereco = emu + ehad; 
  return ereco;

}
//***********************************************
void MadTestAnalysis::MakeMyFile(std::string tag,int tarpos){

  std::string savename = "DPHistos_" + tag + ".root";
  TFile save(savename.c_str(),"RECREATE"); 
  save.cd();
  
  //#declare lots of histos:

  TH1F *evlength_nc = new TH1F("Event length - nc",
                               "Event length - nc",
                               60,0,600);
  evlength_nc->SetXTitle("Event length (planes)");
  TH1F *evlength_cc = new TH1F("Event length - cc",
                               "Event length - cc",
                               60,0,600);
  evlength_cc->SetXTitle("Event length (planes)");


  TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
                               "Track ph frac - nc",
                               50,0,1);
  phfrac_nc->SetXTitle("pulse height fraction");


  TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
                               "Track ph frac - cc",
                               50,0,1);
  phfrac_cc->SetXTitle("pulse height fraction");


  TH1F *phplane_nc = new TH1F("ph per plane (nc)",
                              "ph per plane (nc)",
                               50,0,5000);
  phplane_nc->SetXTitle("pulse height per plane");
  TH1F *phplane_cc = new TH1F("ph per plane (cc)",
                              "ph per plane (cc)",
                               50,0,5000);
  phplane_cc->SetXTitle("pulse height per plane");



  for(int i=0;i<Nentries;i++){

    
    if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
                            << "% done" << std::endl;
   if(!this->GetEntry(i)) continue;

    Int_t nevent = eventSummary->nevent;
    if(nevent==0) continue;
    
    Int_t trkcount=0;
    Int_t shwcount=0;

        
    //fill tree once for each reconstructed event
    for(int j=0;j<eventSummary->nevent;j++){

      if(!LoadEvent(j)) continue; //no event found

      Int_t ntrack = 0; 
      ntrack=ntpEvent->ntrack;
      Int_t nshower = 0;
      nshower=ntpEvent->nshower;

      Float_t totph=0;


      for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
        LoadTrack(ii);

        totph+=ntpTrack->ph.sigcor;
        LoadTHTrack(ii);

      }

      trkcount+=ntrack;


       for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
        LoadShower(ii);

        totph+=ntpShower->ph.sigcor;
        }

       shwcount+=nshower;
 


      Int_t cc_nc = 0;
      Int_t flavour = 0;
      Int_t detector = -1;
      Int_t is_fid = 0;
      Double_t neu_e  = -1; 
      Double_t weight = -1; 
      
      //get detector type:
      if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
        detector = 2;
        //fiducial criteria on vtx for far detector
        if (evlength_cc->GetNbinsX()==60) {
          evlength_cc->SetBins(50,0,500);
          evlength_nc->SetBins(50,0,500);
        }
        is_fid = 1;
        if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
           (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
           sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
                +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0;
      }
      else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
        detector = 1;   

        if (evlength_cc->GetBinLowEdge(50)>200) {
          evlength_cc->SetBins(60,0,300);
          evlength_nc->SetBins(60,0,300);
        }

        //fiducial criteria on vtx for near detector
        is_fid = 1;
        if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 ||
           sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) +
                ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) is_fid = 0;
      }
      
      //check for a corresponding mc event
 
      Int_t mcevent=-1;

      if(LoadTruthForRecoTH(j,mcevent)) {
        //true info for tree:

        cc_nc             = IAction(mcevent);
        flavour           = Flavour(mcevent);
        neu_e             = TrueNuEnergy(mcevent);  
                          
   // tarpos   1 = LE 2 = ME 3 = HE     
          if(tarpos==2)  weight=1;
          if(tarpos==1)  weight=1;      
          if(tarpos==3)  weight=1; 
      }

      int track_index   = -1;
      LoadLargestTrackFromEvent(j,track_index);

      if (track_index!=-1) {

        Int_t trkpass=ntpTrack->fit.pass;

        if (is_fid && trkpass) {

        Float_t evlength=0;
        Float_t phfrac=0;
        Float_t phplane=0;

        evlength=ntpEvent->plane.n;
        if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
          evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
        }

        Float_t phtrack=ntpTrack->ph.sigcor;
        Float_t phevent=ntpEvent->ph.sigcor;
        if (phevent>0) {phfrac=phtrack/phevent;}

        Float_t trkplane=ntpTrack->plane.n;
        if (trkplane>0) {phplane=phtrack/trkplane;}
       
        if(flavour==2 && cc_nc==1) {
          evlength_cc->Fill(evlength,weight);
          phfrac_cc->Fill(phfrac,weight);
          phplane_cc->Fill(phplane,weight);      
        }

        else if(cc_nc==0) {
          evlength_nc->Fill(evlength,weight);
          phfrac_nc->Fill(phfrac,weight);
          phplane_nc->Fill(phplane,weight);      
        }


        }
      }
    }
  }
  //#close file  
  save.Write();
  save.Close();
}
//************************************************
Float_t MadTestAnalysis::SingleTimeFrame(Int_t snarlentry,Int_t nentries){

 Int_t tf=0,tfbef=0,tfaft=0;
   
   
   if(this->GetEntry(snarlentry))                                tf      = ntpHeader->GetTimeFrame();
   if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1))   tfaft   = ntpHeader->GetTimeFrame();   
   if(snarlentry>0        && this->GetEntry(snarlentry-1))       tfbef   = ntpHeader->GetTimeFrame();
 
 Float_t singletimeframe=0;
 if(snarlentry<(nentries-1) && snarlentry>0 ){
  if(tf!=tfaft && tf!=tfbef) singletimeframe=1;
 }
 if(snarlentry==(nentries-1)) {
  if(tf!=tfbef) singletimeframe=1; 
 }  
 if(snarlentry==0) {
  if(tf!=tfaft) singletimeframe=1;
 }

 this->GetEntry(snarlentry);
 return singletimeframe;  
 
}
//********************************************************
void MadTestAnalysis::CreatePAN(std::string tag, Int_t scanfilter){

  std::string savename = "PAN_" + tag + ".root";
  TFile save(savename.c_str(),"RECREATE"); 
  save.cd();
  
  GnumiInterface gnumi;
  if(!gnumi.Status()) {
    cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
         << " Will not be filling NuParent info"
         << endl;
    cout << "Set environmental variable $GNUMIAUX to point to the "
         << "directory containing the gnumi files to fix this!"
         << endl;

  }

 SpillInfo pppinfo;
  
  //PAN Quantities 
  //Truth:
  Float_t true_enu;       // true neutrino energy (GeV)
  Float_t true_pxnu;      // true neutrino momentum-x (GeV)
  Float_t true_pynu;      // true neutrino momentum-y (GeV)
  Float_t true_pznu;      // true neutrino momentum-z (GeV)
  Float_t true_etgt;       // true target energy (GeV)
  Float_t true_pxtgt;      // true target momentum-x (GeV)
  Float_t true_pytgt;      // true target momentum-y (GeV)
  Float_t true_pztgt;      // true target momentum-z (GeV)
  Float_t true_emu;       // true muon energy
  Float_t true_eshw;      // true shower energy
  Float_t true_x;         // true x
  Float_t true_y;         // true y
  Float_t true_q2;        // true q2
  Float_t true_w2;        // true w2
  Float_t true_dircosneu; // true muon dircos w.r.t neutrino
  Float_t true_dircosz;   // track z-direction cosine
  Float_t true_vtxx;      // true x vtx
  Float_t true_vtxy;      // true y vtx
  Float_t true_vtxz;      // true z vtx

  //For Neugen:
  Int_t flavour;          // true flavour: 1-e 2-mu 3-tau
  Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
  Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
  Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
  Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
  Int_t had_fs;           // hadronic final state, number between 200-300
  
  //For Beam Reweighting:
  NuParent *nuparent = new NuParent();

  //Reco Quantities
  Int_t detector;         // Near = 1, Far = 2;
  Int_t run;              // run number
  Int_t snarl;            // snarl number
  Int_t nevent;           // number of events in snarl
  Int_t event;            // event index
  Int_t subrun;           // subrun number     
  Int_t trigsrc;          // trigger source    
  Int_t mcevent;          // mc event index
  Int_t ntrack;           // number of tracks in event
  Int_t nshower;          // number of showers in event

  Int_t is_fid;           // pass fid cut. true = 1, false = 0
  Int_t is_cev;           // fully contained. true = 1, false = 0 
  Int_t is_mc;            // is a corresponding mc event found
  Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
  Float_t pid0;           // pid parameter (usual method)
  Float_t pid1;           // pid parameter (alternative method 1)
  Float_t pid2;           // pid parameter (alternative method 2)
  Float_t annpid;        // ANN PID
  Int_t pass;             // pass analysis cuts. true = 1, false = 0


  Float_t npid0;           // pid parameter (usual method)
  Float_t npid1;           // pid parameter (alternative method 1)
  Float_t npid2;           // pid parameter (alternative method 2)
  Double_t nannpid;        // ANN PID


  Float_t reco_enu;       // reco neutrino energy
  Float_t reco_emu;       // reco muon energy
  Float_t reco_eshw;      // reco shower energy  (shw.ph.gev/1.23)
  Float_t reco_eshw_sqrt; // reco shower energy using sqrt method
  Float_t reco_qe_enu;    // reco qe neutrino energy
  Float_t reco_qe_q2;     // reco qe q2
  Float_t reco_y;         // reco y
  Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
  Float_t reco_dircosz;   // reco track vtx z-dircos
  Float_t reco_vtxx;      // reco x vtx
  Float_t reco_vtxy;      // reco y vtx
  Float_t reco_vtxz;      // reco z vtx

  Float_t evtlength;      // reco event length
  Float_t evtph;          // reco event ph
  Float_t trklength;      // reco track length
  Float_t trkmomrange;    // reco track momentum from range
  Float_t trkqp;          // reco track q/p
  Float_t trkeqp;         // reco track q/p error
  Float_t trkphfrac;      // reco pulse height fraction in track
  Float_t trkphplane;     // reco average track pulse height/plane
  Float_t trkds;          // track DS 
  Float_t rawph;          // Raw ph
// Additional track info
  Float_t trkendz;  // track end Z position
  Float_t trkendx;  // track end X position
  Float_t trkendy;  // track end Y position
  Float_t trkendu;  // track end U position
  Float_t trkendv;  // track end V position
  Float_t trkvtxz;  // track begin Z position
  Float_t trkvtxx;  // track begin X position
  Float_t trkvtxy;  // track begin Y position
  Float_t trkvtxu;  // track begin U position
  Float_t trkvtxv;  // track begin V position

  // SCAN DECISION   // <- define an instance of ScanList class here
  Int_t scandecision;
  ScanList scan;
  
// EVENT TIMING
   Double_t evttimemin;  //MIN STRIP TIME OF EVENT
   Double_t evttimemax;  //MAX STRIP TIME OF EVENT
   Double_t snarlt;

// WEIGHT FOR THE PME BEAM
   Double_t w_le_me;  
// BEAM INFO
   Float_t  pot,horn,tar,vvpos,hhpos,magn;    
// Beam INFO DATABASE
   Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb;
   Double_t timedb;
// TIME INFO DATABASE
   Double_t neartdb,fartdb,neardifftdb;

// IS SNARL CUT IN PIECES 
   Float_t singletimeframe;
// Beam Info Block    
   SpillInfoBlock *spinfo = new SpillInfoBlock();  
   spinfo->Zero(); 
   singletimeframe=-999;
    
  //Trees
  TTree *tree = new TTree("pan","pan");
  //Truth
  tree->Branch("true_enu",&true_enu,"true_enu/F",512000);
  tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000);
  tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000);
  tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000);
  tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000);
  tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000);
  tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000);
  tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000);
  tree->Branch("true_emu",&true_emu,"true_emu/F",512000);
  tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000);
  tree->Branch("true_x",&true_x,"true_x/F",512000);
  tree->Branch("true_y",&true_y,"true_y/F",512000);
  tree->Branch("true_q2",&true_q2,"true_q2/F",512000);
  tree->Branch("true_w2",&true_w2,"true_w2/F",512000);
  tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000);
  tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000);
  tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000);
  tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000);
  tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000);
  //Neugen
  tree->Branch("flavour",&flavour,"flavour/I",512000);
  tree->Branch("process",&process,"process/I",512000);
  tree->Branch("nucleus",&nucleus,"nucleus/I",512000);
  tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000);
  tree->Branch("initial_state",&initial_state,"initial_state/I",512000);
  tree->Branch("had_fs",&had_fs,"had_fs/I",512000);
  //NuParent
  tree->Branch("nuparent","NuParent",&nuparent,8000,1);
  //Reco
  tree->Branch("detector",&detector,"detector/I",512000);
  tree->Branch("run",&run,"run/I",512000);
  tree->Branch("snarl",&snarl,"snarl/I",512000);
  tree->Branch("event",&event,"event/I",512000);
  tree->Branch("mcevent",&mcevent,"mcevent/I",512000);
  tree->Branch("ntrack",&ntrack,"ntrack/I",512000);
  tree->Branch("nshower",&nshower,"nshower/I",512000);
  tree->Branch("subrun",&subrun,"subrun/I",512000);   
  tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000); 
  
  tree->Branch("is_fid",&is_fid,"is_fid/I",512000);
  tree->Branch("is_cev",&is_cev,"is_cev/I",512000);
  tree->Branch("is_mc",&is_mc,"is_mc/I",512000);
  tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000);
  tree->Branch("pid0",&pid0,"pid0/F",512000);
  tree->Branch("pid1",&pid1,"pid1/F",512000);
  tree->Branch("pid2",&pid2,"pid2/F",512000);
  tree->Branch("annpid",&annpid,"annpid/F",512000);  
  tree->Branch("pass",&pass,"pass/I",512000);

  tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000);
  tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000);
  tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000);
  tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000);
  tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000);
  tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000);
  tree->Branch("reco_y",&reco_y,"reco_y/F",512000);
  tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000);
  tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000);
  tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000);
  tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000);
  tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000);

  tree->Branch("evtlength",&evtlength,"evtlength/F",512000);
  tree->Branch("evtph",&evtph,"evtph/F",512000);
  tree->Branch("trklength",&trklength,"trklength/F",512000);
  tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000);
  tree->Branch("trkqp",&trkqp,"trkqp/F",512000);
  tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000);
  tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000);
  tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000);
  tree->Branch("rawph",&rawph,"rawph/F",512000);    // <- !! added this
  
  tree->Branch("w_le_me",&w_le_me,"w_le_me/D",512000); 
  tree->Branch("trkendz",&trkendz,"trkendz/F",512000);
  tree->Branch("trkendu",&trkendu,"trkendu/F",512000);
  tree->Branch("trkendv",&trkendv,"trkendv/F",512000);
  tree->Branch("trkendx",&trkendx,"trkendx/F",512000);
  tree->Branch("trkendy",&trkendy,"trkendy/F",512000);

  tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000);
  tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000);
  tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000);
  tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000);
  tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000);
  tree->Branch("trkds",  &trkds,"trkds/F",512000);
  
  tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
  tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
  
  // SCAN DECISIONS 
  tree->Branch("scandecision",&scandecision,"scandecision/I");

  //BEAM INFO 
  tree->Branch("pot",    &pot,    "pot/F",512000);
  tree->Branch("horn",   &horn,   "horn/F",512000);
  tree->Branch("tar",    &tar,    "tar/F",512000);
  tree->Branch("vvpos",  &vvpos,  "vvpos/F",512000);
  tree->Branch("hhpos",  &hhpos,  "hhpos/F",512000);
  tree->Branch("magn",   &magn,   "magn/F",512000);
  tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000);

  // BEAM INFO DB
  tree->Branch("potdb",     &potdb,      "potdb/F",512000);
  tree->Branch("horndb",    &horndb,     "horndb/F",512000);
  tree->Branch("tardb",     &tardb,      "tardb/F",512000);
  tree->Branch("vvposdb",   &vvposdb,    "vvposdb/F",512000);
  tree->Branch("hhposdb",   &hhposdb,    "hhposdb/F",512000);
  tree->Branch("vvwidthdb", &vvwidthdb,  "vvwidthdb/F",512000);
  tree->Branch("hhwidthdb", &hhwidthdb,  "hhwidthdb/F",512000);
  tree->Branch("timedb",    &timedb,     "timedb/D");

  // TIME DB  
  tree->Branch("neartdb",      &neartdb,    "neartdb/D");
  tree->Branch("fartdb",       &fartdb,     "fartdb/D");
  tree->Branch("neardifftdb",  &neardifftdb,"neardifftdb/D"); 
  tree->Branch("snarlt",       &snarlt,     "snarlt/D");
  
  // new pid variables
  tree->Branch("npid0",&npid0,"npid0/F",512000);
  tree->Branch("npid1",&npid1,"npid1/F",512000);
  tree->Branch("npid2",&npid2,"npid2/F",512000);
  tree->Branch("nannpid",&nannpid,"nannpid/D",512000);  
  

  

  // MAIN LOOP (over tree entries)

  for(int i=0;i<Nentries;i++){
 
    if(i%100==0) std::cout << float(i)*100./float(Nentries) 
                            << "% done" << std::endl;

    if(!this->GetEntry(i)) continue;

    nevent = eventSummary->nevent;
    if(nevent==0) continue;
    
    run     = ntpHeader->GetRun();
    snarl   = ntpHeader->GetSnarl();
    subrun  = ntpHeader->GetSubRun();   
    trigsrc = ntpHeader->GetTrigSrc(); 
    
    Double_t snarltime=ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);       
    Int_t    month    =eventSummary->date.month; 
    Int_t    year     =eventSummary->date.year;

    // SCAN INFO  
    
    scandecision=scan.GetDecision(run,snarl);

    // only pass through scanned snarls if scanfilter=1
    if (scanfilter==0 || scandecision>0) {  

      // Get Beam Monitorint Info if NOT MC: 
      
      if(ntpHeader->GetVldContext().GetSimFlag()!=4){
        pppinfo.GetSpillInfo(month,year,snarltime,*spinfo); 
        pot   =spinfo->GetPot();   
        horn  =spinfo->GetHorn();
        hhpos =spinfo->GetHpos();
        vvpos =spinfo->GetVpos();
        tar   =spinfo->GetTgtpos();
        magn  =spinfo->GetMagnet();   
        
        potdb       = -999;
        horndb      = -999; 
        tardb       = -999; 
        vvposdb     = -999; 
        hhposdb     = -999; 
        vvwidthdb   = -999; 
        hhwidthdb   = -999; 
        timedb      = -999; 
        neartdb     = -999; 
        fartdb      = -999; 
        neardifftdb = -999; 
        
        // BEAM INFO DB
        Detector::Detector_t detectort;
        SimFlag::SimFlag_t simflag;
        VldTimeStamp  timestamp;
        VldContext    cx;
        cx        = ntpHeader->GetVldContext(); 
        detectort = ntpHeader->GetVldContext().GetDetector();
        simflag   = ntpHeader->GetVldContext().GetSimFlag();
        timestamp = ntpHeader->GetVldContext().GetTimeStamp();
        
        BDSpillAccessor &spillAccess=BDSpillAccessor::Get();
        const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp);
        if(spillInfoDB) { 
          if(spillInfoDB->fTortgt !=0.)      potdb     = spillInfoDB->fTortgt;
          else if(spillInfoDB->fTrtgtd !=0.) potdb     = spillInfoDB->fTrtgtd;
          else if(spillInfoDB->fTor101 !=0.) potdb     = spillInfoDB->fTor101;
          else potdb = -999;
          
          horndb    = spillInfoDB->fHornCur;
          tardb     = (int)spillInfoDB->GetStatusBits().target_in;
          vvposdb   = spillInfoDB->fTargBpmY[0];
          hhposdb   = spillInfoDB->fTargBpmX[0];
          vvwidthdb = spillInfoDB->fProfWidY;
          hhwidthdb = spillInfoDB->fProfWidX;
          timedb    = spillInfoDB->SpillTime();
        }
        
        // TIME DB
        SpillTimeFinder& spillfinder = SpillTimeFinder::Instance();
        fartdb                       = spillfinder.GetTimeOfNearestSpill(cx);
        const SpillTimeND& spnd      = spillfinder.GetNearestSpill(cx);
        neartdb                      = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9);
        neardifftdb                  = spillfinder.GetTimeToNearestSpill(cx);
        snarlt                       = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
        
      }
      singletimeframe=SingleTimeFrame(i,Nentries);
      
      Int_t trkcount=0;
      Int_t shwcount=0;
      Float_t totph=0;
      
      // EVENT LOOP - fill tree once for each reconstructed event

      for(int j=0;j<nevent;j++){ 

        if(!LoadEvent(j)) continue; //no event found

        totph = 0;

        //set event index:
        event = j;
        ntrack = ntpEvent->ntrack;
        nshower = ntpEvent->nshower;
        
        // Initialize everything
        nuparent->Zero();     
        true_enu = 0; true_emu = 0; true_eshw = 0; 
        true_pxnu = 0; true_pynu = 0; true_pznu = 0;
        true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
        true_dircosneu = 0; true_dircosz = 0;
        true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
        flavour = 0; process = 0; nucleus = 0; cc_nc = 0;
        initial_state = 0; had_fs = 0;
        true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;     
        detector = -1; mcevent = -1;
        is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0; 
        pid0 = -999; pid1 = -999; pid2 = -999; pass = 0;
        reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0; 
        reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0; 
        reco_dircosz = 0;
        reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999;
        evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1;
        trkmomrange = -999; trkqp = -999; trkeqp = -999;
        trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.;   
        w_le_me=-1; trkendz=999; trkendu=999; trkendv=999; trkendx=999;
        trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999;
        trkds=-1; evttimemin=-1; evttimemax=-1;
        evtph=0.;
        annpid=-999.;
        npid0 = -999; npid1 = -999; npid2 = -999;
        nannpid=-999;

        rawph=ntpEvent->ph.raw;
 
        //get detector type:
        if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
          detector = 2;
          totph=eventSummary->ph.sigcor;
          //fiducial criteria on vtx for far detector
          is_fid = 1;
          if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
        }
        else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
          detector = 1;
        
          //get total charge in trk/shw
          for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
            LoadTrack(ii);
            totph+=ntpTrack->ph.sigcor;
          }
          trkcount+=ntrack;
          
          for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
            LoadShower(ii);
            totph+=ntpShower->ph.sigcor;
          }
          shwcount+=nshower;
          
          //fiducial criteria on vtx for near detector
          is_fid = 1;
          if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
        }
        
        //check for a corresponding mc event      
        if(LoadTruthForRecoTH(j,mcevent)) {
          Float_t *vtx = TrueVtx(mcevent);
          Float_t *nu3mom = TrueNuMom(mcevent);
          Float_t *targ4vec = Target4Vector(mcevent);
          //true info for tree:
          is_mc             = 1;
          nucleus           = Nucleus(mcevent);
          flavour           = Flavour(mcevent);
          initial_state     = Initial_state(mcevent);
          had_fs            = HadronicFinalState(mcevent);
          process           = IResonance(mcevent); 
          cc_nc             = IAction(mcevent);
          true_enu          = TrueNuEnergy(mcevent); 
          true_pxnu         = nu3mom[0];
          true_pynu         = nu3mom[1];
          true_pznu         = nu3mom[2];
          true_emu          = TrueMuEnergy(mcevent); 
          true_eshw         = TrueShwEnergy(mcevent); 
          true_x            = X(mcevent);
          true_y            = Y(mcevent);
          true_q2           = Q2(mcevent);
          true_w2           = W2(mcevent);
          true_dircosz      = TrueMuDCosZ(mcevent);
          true_dircosneu    = TrueMuDCosNeu(mcevent);
          true_vtxx         = vtx[0];
          true_vtxy         = vtx[1];
          true_vtxz         = vtx[2];
          true_etgt         = targ4vec[3];
          true_pxtgt        = targ4vec[0];
          true_pytgt        = targ4vec[1];
          true_pztgt        = targ4vec[2];
          
          if(gnumi.Status()) {
            if(isST) gnumi.GetParent(strecord,mcevent,*nuparent);
            else     gnumi.GetParent(mcrecord,mcevent,*nuparent);
          }
          
          delete [] vtx;
          delete [] nu3mom;
          delete [] targ4vec;
        }
        
        
        if(PassBasicCuts(event)) {pass_basic=1;} // make this reflect pass/fail flag
        
        //reco info for tree:
        reco_vtxx         = ntpEvent->vtx.x;
        reco_vtxy         = ntpEvent->vtx.y;
        reco_vtxz         = ntpEvent->vtx.z;
        evtlength         = ntpEvent->plane.end-ntpEvent->plane.beg+1;
        evtph             = ntpEvent->ph.sigcor; 
        
        // different definition for neardet - accounts for 
        // uninstrumented planes
   
        // ANN PID // 
        
        AnnInputBlock anninput=AnnVar(event);
        //      AnnInputBlock anninput2=anninput;
        //      anninput2.aTotrk=1000;
        //      if(!nsid.CompareAnnBlocks(this,event,anninput2)) assert(false);

        Int_t target=0;
        if(ntpHeader->GetVldContext().GetSimFlag()!=4){
          if(tar<900) target=1;
          if(tar>900 && tar<40000.) target=2;
          if(tar>40000 && tar <100000) target=3;
        }
        if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){
          Int_t mod=(int)(run/1e6);
          if(mod==14) target=1; 
          if(mod==16) target=2;
          if(mod==18) target=3; 
        }
        if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){
          target=0; // FOR THE MOMENT IN THE ABSENCE OF ANY FAR PME PHE MC 
        }                       
        // SET FIDUCIAL AND TARGET DEFAULT IS DAVIDS FIDUCIAL AND NO OSCILLATION TRAINING
        
        Int_t fid    =2;
        Int_t prior = 1;
        
        if(!PassFid(event)) { 
          annpid=-999;
          nannpid=-999;
        }
        else {
          Detector::Detector_t det=
            ntpHeader->GetVldContext().GetDetector();
          // ANN PID //           
          AnnInputBlock anninput=AnnVar(event);  
          if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
            annpid=GetAnnPid(anninput,detector,target,fid,prior);
            if(!nsid.GetPID(this,event,det,nannpid)) nannpid=-999;
            
            if(!((reco_vtxz>1. &&reco_vtxz<5.)&&sqrt((reco_vtxx-1.4885)*(reco_vtxx-1.4885)+(reco_vtxy-0.1397)*(reco_vtxy-0.1397))<1)) {
              annpid=-999.; nannpid=-999;
            }
          } 
          
          if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
            annpid=GetAnnPid(anninput,detector,target,fid,prior);  
            if(!nsid.GetPID(this,event,det,nannpid)) nannpid=-999;
            
            if(fid==1) if(!(((reco_vtxz>1. &&reco_vtxz<14.)||(reco_vtxz>17. &&reco_vtxz<27.))&& sqrt(reco_vtxx*reco_vtxx+reco_vtxy*reco_vtxy)<3.5)) {
              annpid=-999.; nannpid=-999;
            }
          }
        }
        
        Double_t timemin=9.*10e30; 
        Double_t timemax=-9999999; 
        Double_t trgtime=eventSummary->trigtime;

        //  Find max min time of events by loading strips for ND 
        if(detector==1){
          for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
            Int_t stp_index=((ntpEvent->stp)[evss]);
            LoadStrip(stp_index); 
            Double_t striptime=ntpStrip->time1;
            if(striptime<=timemin) timemin=striptime;
            if(striptime>=timemax) timemax=striptime;
          }
          evttimemax=timemax-trgtime;
          evttimemin=timemin-trgtime;  
        }        
        
        if(detector==2){  
          for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
            Int_t stp_index=((ntpEvent->stp)[evss]);
            LoadStrip(stp_index);
            Double_t striptime1=ntpStrip->time1;
            Double_t striptime0=ntpStrip->time0;
            Double_t striptime = 0;
            if(striptime1>0 && striptime0<0)  striptime=striptime1;
            if(striptime0>0 && striptime1<0)  striptime=striptime0;
            if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
            if(striptime<=timemin) timemin=striptime;
            if(striptime>=timemax) timemax=striptime;
          }
          evttimemax=timemax-trgtime;
          evttimemin=timemin-trgtime;
        }
  
        //index will -1 if no track/shower present in event
        int track_index   = -1;
        LoadLargestTrackFromEvent(j,track_index);
        int shower_index  = -1;
        LoadLargestShowerFromEvent(j,shower_index);
        
        //methods should all be protected against index = -1
        reco_emu          = RecoMuEnergy(0,track_index);
        reco_eshw         = RecoShwEnergy(shower_index);
        reco_eshw_sqrt    = RecoShwEnergySqrt(event);
        reco_enu          = reco_emu + reco_eshw;
        reco_qe_enu       = RecoQENuEnergy(track_index);
        if (reco_enu>0) {
          reco_y          = reco_eshw/reco_enu;
        }
        reco_qe_q2        = RecoQEQ2(track_index);
        reco_dircosz      = RecoMuDCosZ(track_index);
        reco_dircosneu    = RecoMuDCosNeu(track_index);
        is_cev            = IsFidAll(track_index);
        
        if(track_index!=-1){ //check track is present before using ntpTrack
          trklength         = ntpTrack->plane.end-ntpTrack->plane.beg+1; 
          trkmomrange       = ntpTrack->momentum.range;
          trkqp             = ntpTrack->momentum.qp;
          trkeqp            = ntpTrack->momentum.eqp;
          trkendz           = ntpTrack->end.z;
          trkendu           = ntpTrack->end.u;
          trkendv           = ntpTrack->end.v;
          trkendx           = ntpTrack->end.x;
          trkendy           = ntpTrack->end.y;
          
          trkvtxz           = ntpTrack->vtx.z;
          trkvtxu           = ntpTrack->vtx.u;
          trkvtxv           = ntpTrack->vtx.v;
          trkvtxx           = ntpTrack->vtx.x;
          trkvtxy           = ntpTrack->vtx.y;
          trkds             = ntpTrack->ds;               
          Float_t phtrack=ntpTrack->ph.sigcor;
          Float_t phevent=ntpEvent->ph.sigcor;
          if (phevent>0) {trkphfrac=phtrack/phevent;}
          if(ntpTrack->plane.n>0) {
            trkphplane      = ntpTrack->ph.sigcor/ntpTrack->plane.n;
          }
        }
        else {
          trklength         = 0;
          trkmomrange       = 0;
          trkqp             = 0;
          trkeqp            = 0;
          trkphfrac         = 0;
          trkphplane        = 0;
        }

        // only calculate likelihood PID for events that pass event 
        // fiducial cuts and track quality cuts 

        if(PassBasicCuts(event) && PassFid(event)){
          pid0              = PID(event,0);
          pid1              = PID(event,1);
          pid2              = PID(event,2);
          npid0 = dpid.CalcPID(this,event,0);
          npid1 = dpid.CalcPID(this,event,1);
          npid2 = dpid.CalcPID(this,event,2);
          
          if(PassAnalysisCuts(event)) pass = 1;
        }
        else{ 
          pid0            = -999.;
          pid1            = -999.;
          pid2            = -999.;
          pass              =  0;
          npid0           = -999.;
          npid1           = -999.;
          npid2           = -999.;

        }
        tree->Fill();

      } // end of EVENT LOOP
    } //  end of SCANINFO loop 
  } // end of NTUPLE ENTRIES loop

  delete nuparent;
  delete spinfo;
  
  save.cd();
  tree->Write();
  save.Write();
  save.Close();
}

//*********************************************************

void MadTestAnalysis::ReadPIDFile(std::string tag){

  fLikeliFile = new TFile(tag.c_str(),"READ");
 
  if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) { //if file is found
 
    fLikeliFile->cd();

    fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc");
    fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc");
    fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc");
    fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc");
    fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)");
    fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)");
    
    for(int i=0;i<6;i++){
      float integ = fLikeliHist[i]->Integral(); 
      fLikeliHist[i]->Scale(1./integ);
    }
  }
  else fLikeliFile = NULL;
}


bool MadTestAnalysis::InFidVol(const Detector::Detector_t& det,
                             const Float_t& x, const Float_t& y, 
                             const Float_t& z){
  // Mike Kordosky: 11 August 2005
  // Defines a fiducial volume used to pass events which fill the PDFs
  // Transcribed from MakeMyFile() above
  
  bool is_fid=false;
  
  if(det==Detector::kFar){
    
    //fiducial criteria on vtx for far detector
    /*
      if (evlength_cc->GetNbinsX()==60) {
      evlength_cc->SetBins(50,0,500);
      evlength_nc->SetBins(50,0,500);
      }
    */
    is_fid = true;
    if(z<0.5 || z>29.4 ||   //ends
       (z<16.5&&z>14.5) ||  //between SMs
       sqrt((x*x)           //radial cut
            +(y*y))>3.5) is_fid = false;
  }
  else if(det==Detector::kNear){
    
    /*
      if (evlength_cc->GetBinLowEdge(50)>200) {
      evlength_cc->SetBins(60,0,300);
      evlength_nc->SetBins(60,0,300);
      }
    */
    //fiducial criteria on vtx for near detector
    is_fid = true;
    if(z<1 || z>5 ||
       sqrt(((x-1.4885)*(x-1.4885)) +
            ((y-0.1397)*(y-0.1397)))>1) is_fid = false;
  }

  return is_fid;
  
}
    
#endif // #ifdef madtestanalysis_cxx

---