LEMAmbyE50S491.cxx

Go to the documentation of this file.

#include "MCNNAnalysis/LEMAmbyE50S491.h"

// Constructor

LEMAmbyE50S491::LEMAmbyE50S491()
{

  // initialize min and max vectors (for normalisation)
      fVmin[0] = -1;
      fVmax[0] = 1;
      fVmin[1] = -1;
      fVmax[1] = 1;
      fVmin[2] = -1;
      fVmax[2] = 0.999999940395355;
      fVmin[3] = -1;
      fVmax[3] = 0.999999940395355;

      // initialize input variable types
      fType[0] = 'F';
      fType[1] = 'F';
      fType[2] = 'F';
      fType[3] = 'F';

      // initialize constants
      Initialize();

      // initialize transformation
      InitTransform();
}//LEMAmbyE50S491


inline void LEMAmbyE50S491::Initialize()
{
   
   // build network structure
   fLayers = 3;
   fLayerSize[0] = 5; fWeights[0] = new double[5]; 
   fLayerSize[1] = 10; fWeights[1] = new double[10]; 
   fLayerSize[2] = 1; fWeights[2] = new double[1]; 
   // weight matrix from layer 0 to 1
   fWeightMatrix0to1[0][0] = -0.152318321120594;
   fWeightMatrix0to1[1][0] = 0.140038941840388;
   fWeightMatrix0to1[2][0] = 0.150830788727814;
   fWeightMatrix0to1[3][0] = -0.152255144703905;
   fWeightMatrix0to1[4][0] = 0.137836159069207;
   fWeightMatrix0to1[5][0] = -0.121184535878213;
   fWeightMatrix0to1[6][0] = 0.139033891203744;
   fWeightMatrix0to1[7][0] = -0.108956700708734;
   fWeightMatrix0to1[8][0] = 0.349134725553892;
   fWeightMatrix0to1[0][1] = 0.205967629818834;
   fWeightMatrix0to1[1][1] = -0.134852396728545;
   fWeightMatrix0to1[2][1] = -0.178592543312265;
   fWeightMatrix0to1[3][1] = 0.186278178092904;
   fWeightMatrix0to1[4][1] = -0.125205805564187;
   fWeightMatrix0to1[5][1] = 0.16269157723916;
   fWeightMatrix0to1[6][1] = -0.130902963928656;
   fWeightMatrix0to1[7][1] = 0.552392965698308;
   fWeightMatrix0to1[8][1] = 2.00605701594573;
   fWeightMatrix0to1[0][2] = -0.351796753103481;
   fWeightMatrix0to1[1][2] = 0.119495150480259;
   fWeightMatrix0to1[2][2] = 0.234526141140616;
   fWeightMatrix0to1[3][2] = -0.259501720569135;
   fWeightMatrix0to1[4][2] = 0.0971029232485688;
   fWeightMatrix0to1[5][2] = -0.311255332162166;
   fWeightMatrix0to1[6][2] = 0.110460663357879;
   fWeightMatrix0to1[7][2] = -1.42763468116269;
   fWeightMatrix0to1[8][2] = 0.962545781664626;
   fWeightMatrix0to1[0][3] = 0.118397441928664;
   fWeightMatrix0to1[1][3] = -0.0814482182643213;
   fWeightMatrix0to1[2][3] = -0.11415203218872;
   fWeightMatrix0to1[3][3] = 0.119649439353558;
   fWeightMatrix0to1[4][3] = -0.0734675726864178;
   fWeightMatrix0to1[5][3] = 0.0952023767112997;
   fWeightMatrix0to1[6][3] = -0.0783740766209056;
   fWeightMatrix0to1[7][3] = 0.347262763001096;
   fWeightMatrix0to1[8][3] = -0.270506518649011;
   fWeightMatrix0to1[0][4] = -0.0842070177217706;
   fWeightMatrix0to1[1][4] = 0.0156252738586454;
   fWeightMatrix0to1[2][4] = 0.0481085654335284;
   fWeightMatrix0to1[3][4] = -0.0567634251254012;
   fWeightMatrix0to1[4][4] = 0.010838022295969;
   fWeightMatrix0to1[5][4] = -0.0359937871007552;
   fWeightMatrix0to1[6][4] = 0.0136216961169865;
   fWeightMatrix0to1[7][4] = -1.81344309339823;
   fWeightMatrix0to1[8][4] = -2.01702270861376;
   // weight matrix from layer 1 to 2
   fWeightMatrix1to2[0][0] = 0.11386357600865;
   fWeightMatrix1to2[0][1] = -0.0654124264943447;
   fWeightMatrix1to2[0][2] = -0.0759398815499141;
   fWeightMatrix1to2[0][3] = 0.0835250008640337;
   fWeightMatrix1to2[0][4] = -0.0682978470479161;
   fWeightMatrix1to2[0][5] = 0.0451858425686448;
   fWeightMatrix1to2[0][6] = -0.0663313873242806;
   fWeightMatrix1to2[0][7] = -0.446499096504961;
   fWeightMatrix1to2[0][8] = 0.58210274049323;
   fWeightMatrix1to2[0][9] = 0.302181334928035;

}//L4Init

inline double LEMAmbyE50S491::GetMvaValue__( const std::vector<double>& inputValues ) const
{
   if (inputValues.size() != (unsigned int)fLayerSize[0]-1) {
      std::cout << "Input vector needs to be of size " << fLayerSize[0]-1 << std::endl;
      return 0;
   }

   for (int l=0; l<fLayers; l++)
      for (int i=0; i<fLayerSize[l]; i++) fWeights[l][i]=0;

   for (int l=0; l<fLayers-1; l++)
      fWeights[l][fLayerSize[l]-1]=1;

   for (int i=0; i<fLayerSize[0]-1; i++)
      fWeights[0][i]=inputValues[i];

   // layer 0 to 1
   for (int o=0; o<fLayerSize[1]-1; o++) {
      for (int i=0; i<fLayerSize[0]; i++) {
         double inputVal = fWeightMatrix0to1[o][i] * fWeights[0][i];
         fWeights[1][o] += inputVal;
      }
      fWeights[1][o] = ActivationFnc(fWeights[1][o]);
   }
   // layer 1 to 2
   for (int o=0; o<fLayerSize[2]; o++) {
      for (int i=0; i<fLayerSize[1]; i++) {
         double inputVal = fWeightMatrix1to2[o][i] * fWeights[1][i];
         fWeights[2][o] += inputVal;
      }
   }

   return fWeights[2][0];
}//GetMVAvalue

double LEMAmbyE50S491::ActivationFnc(double x) const {
   // hyperbolic tan
   return tanh(x);
}



// Clean up
inline void LEMAmbyE50S491::Clear() 
{
   // nothing to clear
}
   inline double LEMAmbyE50S491::GetMvaValue( const std::vector<double>& inputValues ) const
   {
      // classifier response value
      double retval = 0;

      // classifier response, sanity check first
   
         if (IsNormalised()) {
            // normalise variables
            std::cout << "Doing some normalization I probably shouldn't" << std::endl;
            std::vector<double> iV;
            int ivar = 0;
            for (std::vector<double>::const_iterator varIt = inputValues.begin();
                 varIt != inputValues.end(); varIt++, ivar++) {
               iV.push_back(NormVariable( *varIt, fVmin[ivar], fVmax[ivar] ));
            }
            Transform( iV, -1 );
            retval = GetMvaValue__( iV );
         }
         else {
            std::vector<double> iV;
            int ivar = 0;
            for (std::vector<double>::const_iterator varIt = inputValues.begin();
                 varIt != inputValues.end(); varIt++, ivar++) {
               iV.push_back(*varIt);
            }
            Transform( iV, -1 );
            retval = GetMvaValue__( iV );
         }
     

      return retval;
   }//cleanup



//_______________________________________________________________________
inline void LEMAmbyE50S491::InitTransform_1()
{

   // Normalization transformation, initialisation
   fMin_1[0][0] = 0.0228846855462;
   fMax_1[0][0] = 0.98898011446;
   fMin_1[1][0] = 0.00716135371476;
   fMax_1[1][0] = 0.998298585415;
   fMin_1[2][0] = 0.00716135371476;
   fMax_1[2][0] = 0.998298585415;
   fMin_1[0][1] = 0.0475189238787;
   fMax_1[0][1] = 1;
   fMin_1[1][1] = 0;
   fMax_1[1][1] = 1;
   fMin_1[2][1] = 0;
   fMax_1[2][1] = 1;
   fMin_1[0][2] = 0.019999999553;
   fMax_1[0][2] = 0.980000019073;
   fMin_1[1][2] = 0.019999999553;
   fMax_1[1][2] = 0.980000019073;
   fMin_1[2][2] = 0.019999999553;
   fMax_1[2][2] = 0.980000019073;
   fMin_1[0][3] = 0.665335953236;
   fMax_1[0][3] = 15.999958992;
   fMin_1[1][3] = 0.632570564747;
   fMax_1[1][3] = 15.9997463226;
   fMin_1[2][3] = 0.632570564747;
   fMax_1[2][3] = 15.999958992;

}

//_______________________________________________________________________
inline void LEMAmbyE50S491::Transform_1( std::vector<double>& iv, int cls) const
{
if (cls < 0 || cls > 2) {
   if (2 > 1 ) cls = 2;
   else cls = 2;
}
   for (int ivar=0;ivar<4;ivar++) {
      double offset = fMin_1[cls][ivar];
      double scale  = 1.0/(fMax_1[cls][ivar]-fMin_1[cls][ivar]);
      iv[ivar] = (iv[ivar]-offset)*scale * 2 - 1;
   }
}

//_______________________________________________________________________
inline void LEMAmbyE50S491::InitTransform()
{
   InitTransform_1();
}

//_______________________________________________________________________
inline void LEMAmbyE50S491::Transform( std::vector<double>& iv, int sigOrBgd ) const
{
   Transform_1( iv, sigOrBgd );
}




//Get LEM4 PID value using a NueRecord:
Double_t LEMAmbyE50S491::GetLEMAmbyE50S491(NueRecord *nr){

  //LEM4 PID value
  Double_t lemAmby=-1;

  //Input variables
  Double_t val_mny=-1;
  Double_t val_mfqm=-1;
  Double_t val_fcc=-1;
  Double_t val_recoe=-1;

  //Define the variables (y_cut = 1.0 = no y cut)
  val_mny = LEMStandard::Getymean(nr,1.0);
  val_mfqm = LEMStandard::GetMeanFracQMatched(nr,1.0);
  val_fcc = LEMStandard::GetfracCCy(nr,1.0);
  val_recoe = nr->srevent.phNueGeV;

  //Preselection: return automatic -1 for out of range events.
  if (val_mny<0||val_mfqm<0||val_fcc<0||val_recoe<0) {
    lemAmby=-1;
  } else {
    //Define an input vector with the cariable names.
    std::vector<double> inputVec( 4 );

    inputVec[0]=val_mny;
    inputVec[1]=val_mfqm;
    inputVec[2]=val_fcc;
    inputVec[3]=val_recoe;
      
    //Calculate the LEM discriminant!
   lemAmby = LEMAmbyE50S491::GetMvaValue(inputVec);

  
  //Transform so that the PID value is between 0.0 to 1.0.
  //A1Testrun showed spread from -0.10 to 1.13
  //A2lemAmby = (lemAmby + 0.029)/1.170;
  //A3lemAmby -0.071 -- 1.134
  //A5(A3_9) -0.0598634--1.13229
  //E50S491 -0.0571659--1.05873

   Double_t lowval=-0.0571659;
   Double_t highval=1.05873;
   Double_t midpoint=(lowval+highval)/2.0;
   Double_t denom=(highval-lowval)*1.001; //1.022;
   //lemAmby = (( lemAmby - midpoint) / denom) + 0.5;
   lemAmby = (( lemAmby - lowval ) / denom); 


  }
  //return the PID:
  return lemAmby;

}

//For NueAnaReader, extract the NueRecord and use that function.
Double_t LEMAmbyE50S491::GetLEMAmbyE50S491(NueAnaReader *nueana){

  Double_t lemAmby=-1;
  lemAmby = GetLEMAmbyE50S491(nueana->nuerecord);
  return lemAmby;

}

---