LEMAmbyE50S5111.cxx

Go to the documentation of this file.

#include "MCNNAnalysis/LEMAmbyE50S5111.h"

// Constructor

LEMAmbyE50S5111::LEMAmbyE50S5111()
{

      // 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;
      fVmin[4] = -1;
      fVmax[4] = 1;

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

      // initialize constants
      Initialize();

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


inline void LEMAmbyE50S5111::Initialize()
{

  // build network structure
   fLayers = 3;
   fLayerSize[0] = 6; fWeights[0] = new double[6]; 
   fLayerSize[1] = 12; fWeights[1] = new double[12]; 
   fLayerSize[2] = 1; fWeights[2] = new double[1]; 
   // weight matrix from layer 0 to 1
   fWeightMatrix0to1[0][0] = 0.457164863943827;
   fWeightMatrix0to1[1][0] = 0.0220928136491208;
   fWeightMatrix0to1[2][0] = -0.0220155740490428;
   fWeightMatrix0to1[3][0] = 0.0219686579306614;
   fWeightMatrix0to1[4][0] = 0.0220231333381128;
   fWeightMatrix0to1[5][0] = -0.0220049047643356;
   fWeightMatrix0to1[6][0] = 0.0221341437418853;
   fWeightMatrix0to1[7][0] = -0.0219962650033325;
   fWeightMatrix0to1[8][0] = 0.0220472731270565;
   fWeightMatrix0to1[9][0] = 0.0219508308069218;
   fWeightMatrix0to1[10][0] = 0.0220989327625126;
   fWeightMatrix0to1[0][1] = 3.28101780366812;
   fWeightMatrix0to1[1][1] = -0.0776435671584362;
   fWeightMatrix0to1[2][1] = 0.0776263099380779;
   fWeightMatrix0to1[3][1] = -0.0776158321238156;
   fWeightMatrix0to1[4][1] = -0.0776279984708945;
   fWeightMatrix0to1[5][1] = 0.0776239268662966;
   fWeightMatrix0to1[6][1] = -0.0776528050093853;
   fWeightMatrix0to1[7][1] = 0.0776219972306231;
   fWeightMatrix0to1[8][1] = -0.0776333912002757;
   fWeightMatrix0to1[9][1] = -0.0776118516500025;
   fWeightMatrix0to1[10][1] = -0.0776449347033747;
   fWeightMatrix0to1[0][2] = 1.23693441152423;
   fWeightMatrix0to1[1][2] = -0.13057738176211;
   fWeightMatrix0to1[2][2] = 0.130496678670272;
   fWeightMatrix0to1[3][2] = -0.130447662171314;
   fWeightMatrix0to1[4][2] = -0.130504576619991;
   fWeightMatrix0to1[5][2] = 0.13048553150881;
   fWeightMatrix0to1[6][2] = -0.130620567960954;
   fWeightMatrix0to1[7][2] = 0.130476504870904;
   fWeightMatrix0to1[8][2] = -0.130529798327316;
   fWeightMatrix0to1[9][2] = -0.130429037619701;
   fWeightMatrix0to1[10][2] = -0.130583775555715;
   fWeightMatrix0to1[0][3] = -0.199013184025842;
   fWeightMatrix0to1[1][3] = 0.0361593634871003;
   fWeightMatrix0to1[2][3] = -0.0361047740517259;
   fWeightMatrix0to1[3][3] = 0.0360716192096197;
   fWeightMatrix0to1[4][3] = 0.0361101163128128;
   fWeightMatrix0to1[5][3] = -0.0360972340227785;
   fWeightMatrix0to1[6][3] = 0.0361885765061547;
   fWeightMatrix0to1[7][3] = -0.0360911283636637;
   fWeightMatrix0to1[8][3] = 0.0361271766934121;
   fWeightMatrix0to1[9][3] = 0.0360590217410728;
   fWeightMatrix0to1[10][3] = 0.0361636884876753;
   fWeightMatrix0to1[0][4] = -0.639962278203516;
   fWeightMatrix0to1[1][4] = -0.0519357491034027;
   fWeightMatrix0to1[2][4] = 0.0519940946882203;
   fWeightMatrix0to1[3][4] = -0.0520295581992349;
   fWeightMatrix0to1[4][4] = -0.0519883823676268;
   fWeightMatrix0to1[5][4] = 0.0520021579314653;
   fWeightMatrix0to1[6][4] = -0.0519045490115962;
   fWeightMatrix0to1[7][4] = 0.0520086880586148;
   fWeightMatrix0to1[8][4] = -0.0519701438024577;
   fWeightMatrix0to1[9][4] = -0.0520430382975507;
   fWeightMatrix0to1[10][4] = -0.0519311289034354;
   fWeightMatrix0to1[0][5] = -2.65616566273247;
   fWeightMatrix0to1[1][5] = 0.0242893162510659;
   fWeightMatrix0to1[2][5] = -0.0243270284656713;
   fWeightMatrix0to1[3][5] = 0.0243499119561036;
   fWeightMatrix0to1[4][5] = 0.0243233397486465;
   fWeightMatrix0to1[5][5] = -0.0243322339928733;
   fWeightMatrix0to1[6][5] = 0.0242691172849299;
   fWeightMatrix0to1[7][5] = -0.0243364486492918;
   fWeightMatrix0to1[8][5] = 0.0243115571667766;
   fWeightMatrix0to1[9][5] = 0.0243586025722464;
   fWeightMatrix0to1[10][5] = 0.0242863265606441;
   // weight matrix from layer 1 to 2
   fWeightMatrix1to2[0][0] = 0.331026185122903;
   fWeightMatrix1to2[0][1] = -0.114384460634508;
   fWeightMatrix1to2[0][2] = 0.114382212447269;
   fWeightMatrix1to2[0][3] = -0.114380892573337;
   fWeightMatrix1to2[0][4] = -0.114382428341914;
   fWeightMatrix1to2[0][5] = 0.114381909256857;
   fWeightMatrix1to2[0][6] = -0.114385702040211;
   fWeightMatrix1to2[0][7] = 0.114381665049863;
   fWeightMatrix1to2[0][8] = -0.114383123779133;
   fWeightMatrix1to2[0][9] = -0.114380400108112;
   fWeightMatrix1to2[0][10] = -0.114384642740348;
   fWeightMatrix1to2[0][11] = 0.491114735426463;


}//L4Init

inline double LEMAmbyE50S5111::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 LEMAmbyE50S5111::ActivationFnc(double x) const {
   // hyperbolic tan
   return tanh(x);
}



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

      // classifier response, sanity check first
   
         if (IsNormalised()) {
            // normalise variables
            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 LEMAmbyE50S5111::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;
   fMin_1[0][4] = 0;
   fMax_1[0][4] = 0.998935520649;
   fMin_1[1][4] = 0;
   fMax_1[1][4] = 0.997940421104;
   fMin_1[2][4] = 0;
   fMax_1[2][4] = 0.998935520649;

}

//_______________________________________________________________________
inline void LEMAmbyE50S5111::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<5;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 LEMAmbyE50S5111::InitTransform()
{
   InitTransform_1();
}

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




//Get LEM4 PID value using a NueRecord:
Double_t LEMAmbyE50S5111::GetLEMAmbyE50S5111(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;
  Double_t val_ncmny=-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;
  val_ncmny = LEMStandard::GetNCymean(nr,1.0);

  //Preselection: return automatic -1 for out of range events.
  if (val_ncmny<0||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( 5 );

    inputVec[0]=val_mny;
    inputVec[1]=val_mfqm;
    inputVec[2]=val_fcc;
    inputVec[3]=val_recoe;
    inputVec[4]=val_ncmny;
      
    //Calculate the LEM discriminant!
   lemAmby = LEMAmbyE50S5111::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
  //E50S5111 -0.0997124--1.0242

   Double_t lowval=-0.0997124;
   Double_t highval=1.0242;
   Double_t midpoint=(lowval+highval)/2.0;
   Double_t denom=(highval-lowval)*1.001; //prev 1.022mod
   //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 LEMAmbyE50S5111::GetLEMAmbyE50S5111(NueAnaReader *nueana){

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

}

---