MadNsID.cxx

Go to the documentation of this file.

//
// MadNsID
//
// A class to calculate Niki Saoulidou's ANN based CC PID
//
// Created:  N. Saoulidou --Analysis first presented early 2004, Code implementation in May-2005
//
// Converted to this format: M. Kordosky, 12 August 2005
//
// $Author: niki $ 
//
// $Revision: 1.22 $
// 
// $Name:  $
//
// $Id: MadNsID.cxx,v 1.22 2007/04/25 15:48:53 niki Exp $
//
// $Log: MadNsID.cxx,v $
// Revision 1.22  2007/04/25 15:48:53  niki
// changes for new 7 var ANN without total pulse height in input variables
//
// Revision 1.21  2007/03/19 23:50:45  niki
// New Cedar Daikon ANN Weight files and new 7-Input variable ANN
//
// Revision 1.20  2007/03/15 19:44:54  niki
// Small changes in order to be able to call methods from NCUtils
//
// Revision 1.19  2007/03/14 16:17:58  niki
// Fixed a the beam type definition for LE010
//
// Revision 1.18  2007/02/17 00:32:13  niki
// UPdated ann weights files
//
// Revision 1.17  2007/02/13 01:09:10  niki
// Reading Cedar/Carrot weights for ANN separation
//
// Revision 1.16  2007/02/09 21:33:57  niki
// Josh pointed out "unprotected" occurences of strip index = -! (Thanks
// Josh!) All such occurences are now protected against.
//
// Revision 1.15  2007/02/01 22:01:19  rhatcher
// Changes to ROOT headers means we need an explicit #include "TClonesArray.h"
//
// Revision 1.14  2007/01/30 03:17:21  rhatcher
// No longer user BeamType "class" found in MCReweight, but rather the
// BeamType "namespace" found in Conventions package.
// Also covert to using "Detector::" rather than "DetectorType::", though
// for now they are synonyms.
//
// Revision 1.13  2007/01/29 23:05:13  boehm
// Adjusting how the data files are loaded into the two PIDs.  Now the code looks first in the SRT_PRIVATE_CONTEXT and then checks the SRT_PUBLIC_CONTEXT.  This further detachs the CC-pids from the Mad Framework.
//
// Revision 1.12  2006/10/16 19:53:30  kordosky
// cope with CEDAR strip indices set to -1. allow user to force beam type (useful for 150cm data). Open files with TFile::Open() to allow natural dcache access. IMPORTANT: Big change to pans. Tacking on a spill marker, looks like event=-1 with no muon or shower energy, pass and is_fid set to fail. Beam variables and nevent are correct for these entries. Modification intended to ease POT counting
//
// Revision 1.11  2006/05/08 21:27:44  boehm
// Primary change is the addition of functions to calculate the PID's from NtpSt objects.  Both PIDs may now be generated without relying on the rest of the Mad infrastructure.  Mad based code though should not be changed.
//
// Also removed the if(!ntpTrack.fit.pass) return;  from MadDpID.cxx as it does not reflect the current state of the CC analysis. This was removed in the branched version but remained in the development version.
//
// Revision 1.10  2006/02/13 00:34:22  niki
// update my ANN weight files to version 18_2
//
// Revision 1.9  2005/12/14 23:19:58  niki
//  Fix a bug (ntrklike was zero for a given track(?) and I was dividing by
// it).
//
//  Thanks to David P. and Mike K. for catching this!
//
// Revision 1.8  2005/12/12 15:34:32  niki
//  Revised the ANN weight files in order to reflect training with R1_18
//  Addded a few comments to explain the use of priors, fiducials e.t.c
//
// Revision 1.7  2005/10/07 11:19:12  kordosky
// deal with BeamMonSpill::StatusBits.beam_type==0 without asserting. Clearly, it can be zero!
//
// Revision 1.6  2005/10/06 18:40:52  kordosky
// random additions
//
// Revision 1.5  2005/10/06 15:03:26  kordosky
// various improvements, pids finally integrated into pan construction routine, still not fully debugged though
//
// Revision 1.4  2005/09/14 13:32:45  kordosky
// DpID and NsID classes now fully working and validated with MadTestAnalysis by comparing nannpid vs. annpid and npid0 vs. pid0.  pdf files and weights stored in data subdirectory. FD pdfs needed for DpID.
//
//


#include <fstream>
#include <iostream>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include "TClonesArray.h"
#include "TSystem.h"

#include "Mad/MadNsID.h"
#include "Mad/MadQuantities.h"

MadNsID::MadNsID() : weights_read(false), weight_fname("") {
  cache_det=Detector::kUnknown;
  cache_beam=BeamType::kUnknown;
  cache_fid = -1;
  cache_prior = -1;

}

bool MadNsID::ReadWeights(const char* file){
  // read a text file containing the ANN weights
  // returns true if file was correctly read
  // returns false otherwise
  weights_read = false;

  std::ifstream weightfile(file);
  if(!weightfile){ 
    std::cout<<"MadNsID::ReadWeights(), cannot read: "<<file<<std::endl;
    return false;
  }
  std::cout<<"Reading MadNsID weights from: "<<file<<std::endl;
  
  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;

  // read weights
  for(Int_t k=0;k<all;k++){  
    Double_t var;
    weightfile >> var ;
    if(!weightfile.good()){
      std::cerr<<"MadNsID::ReadWeights() bad input, k="<<k<<std::endl;
      weightfile.close();
      return false;
    }
    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();  
  weights_read=true;
  return true;

}


Double_t MadNsID::CalcPID(const AnnInputBlock& anninput, 
                          const Detector::Detector_t& det){
  // calculate the ANN PID based on the information in the input block
  // code does a little detector specific stuff...
  if(!weights_read) assert( ReadWeights(weight_fname.c_str()) );

  Double_t rin[hidneuron];
  Double_t out[hidneuron];
  for(Int_t i=0; i<hidneuron; i++){rin[i]=0; out[i]=0;}
  
  
  if(det==Detector::kFar){  
  
      out[0]  =  (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/2000.;
      out[1]  =  (anninput.aPhperpl/10000.); 
   // out[1]  =  (anninput.aTotph/50000.); 
      out[2]  =  (anninput.aTotlen/40);
      out[3]  =  (anninput.aPh3/(anninput.aTotph+1.));
      out[4]  =  (anninput.aPh6/(anninput.aTotph+1.));
      out[5]  =  (anninput.aPhlast/(anninput.aTotph+1.));
      out[6]  =  ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
 /*     
    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==Detector::kNear){
       out[0]  =  (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/5000.;
       out[1]  =  (anninput.aPhperpl/10000.); 
     //out[1]  =  (anninput.aTotph/100000.); 
       out[2]  =  (anninput.aTotlen/40);
       out[3]  =  (anninput.aPh3/(anninput.aTotph+1.));
       out[4]  =  (anninput.aPh6/(anninput.aTotph+1.));
       out[5]  =  (anninput.aPhlast/(anninput.aTotph+1.));
       out[6]  =  ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
/*       
    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.);


     
  // 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]);      
  }             

  Double_t prob =constanto[0];  
  for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i];
  
  
  if(anninput.aTotlen>50) prob=0.;

  return prob;

}

bool MadNsID::CalcVars(MadQuantities* mq, Int_t event, 
                       AnnInputBlock& anninput){

  // calculates ann input variables for this event and records them in
  // anninput.  returns false in case of problems, true otherwise
  

  // 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)
  
  const NtpSREvent* ntpEvent = mq->GetEvent(event) ;  
  if(!ntpEvent) return false; // perhaps ought to assert here?
  int track_index   = -1;
  const NtpSRTrack* ntpTrack   = mq->GetLargestTrackFromEvent(event,track_index);
  int shower_index  = -1;
  const NtpSRShower* ntpShower = mq->GetLargestShowerFromEvent(event,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]);

    if(stp_index!=-1){
 
    const NtpSRStrip* ntpStrip =  mq->GetStrip(stp_index);      
    if(!ntpStrip) continue;
    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]);
        
       if(stp_indexs!=-1){

        ntpStrip = mq->GetStrip(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]);

        if(stp_indext!=-1){

        ntpStrip = mq->GetStrip(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 if stp_index != -1

  } // end of looping over event strips....    

//
  anninput.aPh3       =ph3;
  anninput.aPh6       =ph6;
  anninput.aPhlast    =phlast;
  anninput.aPhcommon  =phcommon; 
  
  return true;


}

bool MadNsID::CalcVars(NtpSREvent* ntpEvent, NtpSRTrack *ntpTrack,
                       NtpSRShower *ntpShower, NtpStRecord* st,
                        AnnInputBlock& anninput)
{                                                                            
                                                                                
  // calculates ann input variables for this event and records them in
  // anninput.  returns false in case of problems, true otherwise
                                                                                
  // 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)
                                                                                
  if(!ntpEvent) return false; // perhaps ought to assert here?
                                                                                
  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(ntpTrack != 0) {
    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(ntpShower != 0) {
    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 i = 0; i < ntpEvent->nstrip;i++)
  {
    Int_t stp_index=((ntpEvent->stp)[i]);

    if(stp_index!=-1){

    NtpSRStrip* ntpStrip = dynamic_cast<NtpSRStrip *>((*st->stp)[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(ntpShower != 0) {
      for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
        Int_t stp_indexs=((ntpShower->stp)[jj]);
       
        if(stp_indexs!=-1){

        ntpStrip = dynamic_cast<NtpSRStrip *>((*st->stp)[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(ntpTrack != 0) {
      for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){
        Int_t stp_indext=((ntpTrack->stp)[jj]);
       
        if(stp_indext!=-1){

         ntpStrip = dynamic_cast<NtpSRStrip *>((*st->stp)[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 if strip index !=-1
  } // end of looping over event strips....
//
  anninput.aPh3       =ph3;
  anninput.aPh6       =ph6;
  anninput.aPhlast    =phlast;
  anninput.aPhcommon  =phcommon;
                                                                                
  return true;
}

std::string MadNsID::GetWeightFileName(const Detector::Detector_t& det,
                                       const BeamType::BeamType_t& beam, 
                                       Int_t fid, Int_t prior)
{
  // M. Kordosky Sept 09, 2005
  //
  // Choose a file of ANN weights according to the detector, beam, 
  // fiducial volue and prior 
  // 
  //
  // N Saoulidou December 12 2005
  // 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
                                                                                                                         

  //
  std::string priv="";
  std::string pub="";
                                                                                
  std::string base="";
                                                                                
  priv=getenv("SRT_PRIVATE_CONTEXT");
  pub = getenv("SRT_PUBLIC_CONTEXT");
                                                                                
  if(priv == "" && pub == ""){
    std::cerr<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
    assert(false);
  }
                                                                                
  priv += "/Mad/data/";
  pub  += "/Mad/data/";
                                                                                
  struct stat ss;
  if(stat(priv.c_str(), &ss) == -1) {
   std::cout<<"Data not present in SRT_PRIVATE_CONTEXT trying SRT_PUBLIC"<<std::endl;
   if(stat(pub.c_str(), &ss) == -1) {
     std::cout<<"Data not present in SRT_PUBLIC - doesn't seem to exist"<<std::endl;
      assert(false);
   }else{ base = pub; }
  }else { base = priv; }
  
  if(det==Detector::kFar){  
    if(prior==1){ 
      if(fid==1) base+="weights_farle_cedar_daikon7v2.dat";    
      if(fid==2) base+="weights_farle_cedar_daikon7v2.dat";      
    }
    if(prior==2){ 
      if(fid==1) base+="weights_farle_cedar_daikon7v2.dat";    
      if(fid==2) base+="weights_farle_cedar_daikon7v2.dat";
    }
    if(prior==3){ 
      if(fid==1) base+="weights_farle_cedar_daikon7v2.dat";    
      if(fid==2) base+="weights_farle_cedar_daikon7v2.dat";
    }
  }
  else if(det==Detector::kNear){
    if(beam==BeamType::kL010z185i || beam==BeamType::kL000z200i || beam==BeamType::kL010z170i || beam==BeamType::kL010z200i)  base+="weights_nearle_cedar_daikon7v2.dat";    
    if(beam==BeamType::k100)  base+="weights_nearpme_cedar_daikon7v2.dat";    
    if(beam==BeamType::k250)  base+="weights_nearphe_cedar_daikon7v2.dat";    
  }

  if(stat(base.c_str(),&ss) == -1) {
    std::cerr<<"The file '"<<base<<"' doesn't seem to exist"<<std::endl;
    assert(false);
  }

  return base;
}

bool MadNsID::ChooseWeightFile(const Detector::Detector_t& det,
                               const BeamType::BeamType_t& beam, 
                               Int_t fid, Int_t prior){
  // look for the proper weight file
  // return true if we already have read it
  // or if we haven't read it but can find it
  if(det==cache_det && beam==cache_beam && 
     prior==cache_prior && fid==cache_fid){
    // no action needed
    return true;
  }
  
  weight_fname = GetWeightFileName(det,beam,fid,prior);
  // check if this file or directory exists
  Long_t id,size,flags,modtime;
  if(gSystem->GetPathInfo(weight_fname.c_str(),&id,&size,&flags,&modtime) == 1) {
    return false;
  }
  // check that it is a file and not a directory
  if(flags>1) return false;
  else {
    cache_det=det; cache_beam=beam; cache_prior=prior; cache_fid=fid;
    weights_read=false;
    return true;
  }
}

void MadNsID::SetWeightFile(const char* file){
  weight_fname = file;
}

bool MadNsID::GetPID(MadQuantities* mq, Int_t event, 
                     const Detector::Detector_t & det, Double_t& pid){
  AnnInputBlock ib;
  bool result=true;
  result = CalcVars(mq, event, ib);
  if(result) pid = CalcPID(ib,det);
  return result;
}

bool MadNsID::GetPID(NtpSREvent* ntpEvent, NtpSRTrack *ntpTrack,
                     NtpSRShower *ntpShower, NtpStRecord* st,
                     const Detector::Detector_t& det, Double_t& pid)
{
  AnnInputBlock ib;
  bool result=true;
  result = CalcVars(ntpEvent, ntpTrack, ntpShower, st, ib);
  if(result) pid = CalcPID(ib,det);
  return result;
}

bool MadNsID::CompareAnnBlocks(MadQuantities* mq, Int_t event, const AnnInputBlock& ext){
  AnnInputBlock ib;
  CalcVars(mq, event, ib);
  if( !(ib==ext) ){
    std::cout<<"AnnInputBlocks disagree!!!"<<std::endl;
    std::cout<<"Event: "<<event<<std::endl;
    std::cout<<"This block: "<<std::endl;
    std::cout<<ib<<std::endl;
    std::cout<<"Other block: "<<std::endl;
    std::cout<<ext<<std::endl;
    return false;
  }
  else return true;
}

---