NCExtractionMDA.cxx

Go to the documentation of this file.

#include "NCUtils/Extraction/NCExtractionMDA.h"
#include "NCUtils/NCType.h"

#include "MessageService/MsgService.h"
#include "Conventions/Detector.h"
#include "Conventions/BeamType.h"
#include "Conventions/ReleaseType.h"
#include "Conventions/SimFlag.h"
#include "Util/UtilString.h"
#include "AnalysisNtuples/ANtpDefaultValue.h"

#include "NCUtils/Cuts/NCAnalysisCutsNC.h"

#include "TDirectory.h"
#include "TLeaf.h"
#include "TObjArray.h"

#include <algorithm>
#include <cassert>

#include <iostream>
#include <fstream>
using std::cout;
using std::endl;

#include <fstream>
using std::ifstream;

using std::string;

ClassImp(NCExtractionMDA)

CVSID("$Id: NCExtractionMDA.cxx,v 1.8 2009/09/16 13:51:16 bckhouse Exp $");

#include "NCUtils/Extraction/MicroDSTMaker.h"
REGISTER_NCEXTRACTION(NCExtractionMDA, MDA)

//----------------------------------------------------------------------
NCExtractionMDA::NCExtractionMDA(NCAnalysisCuts* cuts,
                                 const Registry& r) :
  NCExtraction(cuts, r),
  fIsFilled(false),
  fThreshCut(ANtpDefaultValue::kDouble)
{
  //Get macro config parameters
  int tmpi;
  const char* tmps;

  // Need this key for some reason
  assert(r.KeyExists("MDAMCPath"));

  if (r.Get("TreeName",   tmps)) fTreeName   = tmps;
  if (r.Get("MDAMCPath",  tmps)) fFilePath   = tmps;
  if (r.Get("TrainTrueE", tmpi)) fTrainTrueE = tmpi;
}

//-----------------------------------------------------------------------
bool NCExtractionMDA::ReadCalibInfoFromFile(NCEventInfo& evtInfo)
{

  // Try $SRT_LOCAL, $SRT_PRIVATE_CONTEXT then $SRT_PUBLIC_CONTEXT
  TString basedir(getenv("SRT_LOCAL"));
  if(basedir=="") basedir=getenv("SRT_PRIVATE_CONTEXT");
  // SRT_PRIVATE_CONTEXT gets set to '.' instead of empty, for some reason
  if(basedir==".") basedir=getenv("SRT_PUBLIC_CONTEXT");

  TString dataDir = basedir+"/NCUtils/data/";

  Detector::Detector_t detType = (Detector::Detector_t)evtInfo.header->detector;

  fTrainTrueE=false;

  if (fTrainTrueE){
    fCoeffFileName = dataDir +
      "Mda_Coeff_Near_L010z185i_CedarDaikon_TTE.dat";

    fHistDefFileName = dataDir +
      "Mda_HistDef_Near_L010z185i_CedarDaikon_TTE.dat";
  }else{
    fCoeffFileName = dataDir +
      WhichMdaCoeffHistDef(true, detType, fReleaseType);

    fHistDefFileName = dataDir +
      WhichMdaCoeffHistDef(false, detType, fReleaseType);
  }

  //Setup chain for P Vector here, since we only run this once.
  //Could go in constructor but would have to multiply constructor
  //arguments, cleaner this way.
  fChain = new TChain(fTreeName);

  MSG("NCExtractionMDA",Msg::kInfo) << "Trying to add " << fFilePath
                                    << " to MDA chain" << endl;
  fChain->Add(fFilePath);


  //Reads the SAS coefficient input file
  //(containing quadratic, linear and constant coeffs)
  //and fills the necessary matrices and vectors
  //Also fills a vector

  MSG("NCExtractionMDA",Msg::kDebug)<<"Reading Calibration Info! "
                                    << " Should be done once per file. "
                                    << endl;

  static const Int_t kLineSize=4000;

  ifstream sasInput;

  MSG("NCExtractionMDA",Msg::kInfo)<< "MDA Calib File: "
                                   << fCoeffFileName
                                   << endl;

  sasInput.open(fCoeffFileName);
  if (sasInput.fail()) {

    MAXMSG("NCExtractionMDA",Msg::kInfo, 20)<< "SAS coefficient file "
                                            << "not found! "
                                            << "Not calculating MDA PID..."
                                            << endl;
    return false;
  }

  Char_t oneLine[kLineSize];

  string sepaCol=" ";
  vector<string> oneVec;

  //Read the first line and get the Threshold cut value
  sasInput.getline(oneLine,kLineSize);
  string topLine=oneLine;
  UtilString::StringTok(oneVec,topLine,sepaCol);
  // The file has 0=CC, 1=NC - NCUtils convention is the other way round
  fThreshCut = 1 - atof((oneVec.at(2)).c_str());

  MSG("NCExtractionMDA",Msg::kInfo) << "MDA Threshold cut: "
                                    << fThreshCut
                                    << endl;

  //Reset the line vector
  oneVec.clear();

  //Read the second line (mean values) and get the number
  //of variables
  sasInput.getline(oneLine,kLineSize);
  topLine=oneLine;
  UtilString::StringTok(oneVec,topLine,sepaCol);

  //Read the mean values into a vector
  TVectorD meanTemp(oneVec.size()-2);
  for(UInt_t col=2; col < oneVec.size(); col++){
    meanTemp(col-2)=(atof((oneVec.at(col)).c_str()));
  }

  for(Int_t m=0; m < meanTemp.GetNoElements(); m++){
    MSG("NCExtractionMDA",Msg::kDebug)<< "Mean:"
                                      << m
                                      << "="
                                      << meanTemp(m)
                                      << endl;
  }
  fMeanVec.push_back(meanTemp);

  //Declare matrices and vectors for coefficient reading
  TMatrixD quadTemp(oneVec.size()-2, oneVec.size()-2);

  TVectorD linTemp(oneVec.size()-2);

  TVectorD constTemp(oneVec.size()-2);

  vector <string> tempClassVec;

  //Count the rows for each matrix
  Int_t rowCounter=0;

  //Reset the number of classes counter
  fNClass=0;

  //Reset the ifstream pointer to the beginning of the file
  // sasInput.seekg(0,ios::beg);

  //Reset the line vector
  oneVec.clear();

  while (! sasInput.eof() ){
    sasInput.getline(oneLine,kLineSize);
    topLine=oneLine;
    UtilString::StringTok(oneVec,topLine,sepaCol);

    string coeffType;

    for(UInt_t col=0; col < oneVec.size(); col++){

      if (col==0) {

        tempClassVec.push_back(oneVec.at(col));

      }else if (col==1){

        coeffType=oneVec.at(col);

        if(coeffType!="_LINEAR_" && coeffType !="_CONST_"){
          if (fVarNameVec.size()<(oneVec.size()-2))
            fVarNameVec.push_back(oneVec.at(col));
          fTypeCoeffVec.push_back("_QUAD_");
        }else fTypeCoeffVec.push_back(oneVec.at(col));

      }else{

        if (coeffType=="_LINEAR_"){

          linTemp(col-2)=(atof((oneVec.at(col)).c_str()));

        }else if (coeffType=="_CONST_"){

          constTemp(col-2)=(atof((oneVec.at(col)).c_str()));

        }else{
          quadTemp(rowCounter,(col-2))
            =(atof((oneVec.at(col)).c_str()));
        }
      }
    }
    oneVec.clear();
    rowCounter++;

    if (coeffType=="_CONST_"){
      fNClass++;
      rowCounter=0;
      fTypeClassVec.push_back(tempClassVec);
      fQuadCoeff.push_back(quadTemp);
      fLinCoeff.push_back(linTemp);
      fConstCoeff.push_back(constTemp);

      tempClassVec.clear();

    }else continue;
  }

  if(sasInput) sasInput.close();

  MSG("NCExtractionMDA",Msg::kDebug)<<"Completed file Reading" << endl;

  // TODO: Replace this with NCEventInfo's SetChainBranches or whatever it's called
  //Associate branches with chain
  fChain->SetBranchAddress("header.", &evtInfo.header);
  fChain->SetBranchAddress("event.",  &evtInfo.event);
  fChain->SetBranchAddress("shower.", &evtInfo.shower);
  fChain->SetBranchAddress("track.",  &evtInfo.track);

  //Recover ANtp style variable names from SAS format.
  //Best to do this outside main loop
  std::vector < std::string> varNameVecDot(fVarNameVec.size());

  for(UInt_t i = 0; i < fVarNameVec.size();i++){
    string varNameTemp=fVarNameVec.at(i);

    //Replace "_" with "." in the ntuple variable names
    if (varNameTemp.length() > varNameTemp.find_first_of("_"))
      varNameTemp.replace(varNameTemp.find_first_of("_"),1,".");

    //Handle var names with additional "Info_"
    if (varNameTemp.length() > varNameTemp.find("Info_"))
      varNameTemp.replace(varNameTemp.find_first_of("_"),1,".");

    varNameVecDot.at(i)=(varNameTemp);

  }

  // Loop once over all chain leaves and get index vector for the
  // relevant variables

  // Get complete list of variables.
  TObjArray *fListLeaves;

  fListLeaves=fChain->GetListOfLeaves();

  //Get Iterator on list of leaves
  TIter *leafIter= new TIter(fListLeaves);
  TLeaf* thisLeaf;
  Int_t iterIndex=0;

  leafIter->Reset();

  fLeafIndexVec.resize(varNameVecDot.size());

  // Note: Order of variables as read into fVarNameVec must be
  // preserved or P vector will not match calib coefficients
  while ((thisLeaf = dynamic_cast<TLeaf*>(leafIter->Next()))){

    string varName=thisLeaf->GetName();

    for(UInt_t i = 0; i < varNameVecDot.size();i++){

      if(varName == varNameVecDot.at(i)){
        fLeafIndexVec.at(i)=iterIndex;
      }
    }
    iterIndex++;
  }

  return true;
}

void NCExtractionMDA::FillPVector(NCEventInfo& evtInfo)
{

  MSG("NCExtractionMDA",Msg::kDebug)<<"Filling P Vector" << endl;

  //Now use fVarNameVec to fill the variable values p-vector
  if(fVarNameVec.size() == 0) return;

  Double_t valueF;
  Int_t valueI;

  TVectorD valueTemp(fVarNameVec.size());

  //Associate branches with chain
  // TODO replcce with the NCEventInfo member function
  fChain->SetBranchAddress("header.", &evtInfo.header);
  fChain->SetBranchAddress("event.",  &evtInfo.event);
  fChain->SetBranchAddress("shower.", &evtInfo.shower);
  fChain->SetBranchAddress("track.",  &evtInfo.track);

  //Get complete list of variables.
  TObjArray *fListLeaves;
  fListLeaves=fChain->GetListOfLeaves();

  //Get relevant leaves from index vector
  TLeaf* thisLeaf;

  for(UInt_t i = 0; i < fLeafIndexVec.size();i++){
    thisLeaf = dynamic_cast<TLeaf*>(fListLeaves->At(fLeafIndexVec.at(i)));

    MSG("NCExtractionMDA",Msg::kDebug)<< i <<" - Found variable "
                                      << "ANtp: " << thisLeaf->GetName()
                                      << endl;


    if(strcmp(thisLeaf->GetTypeName(),"Float_t")==0
       || strcmp(thisLeaf->GetTypeName(),"float")==0
       || strcmp(thisLeaf->GetTypeName(),"Double_t")==0
       || strcmp(thisLeaf->GetTypeName(),"double")==0){

      valueF=thisLeaf->GetValue();

      if (valueF > 1e+8 || ANtpDefVal::IsDefault(valueF)
          || ANtpDefVal::IsDefault(static_cast<Float_t>(valueF))
          || !TMath::Finite(valueF) || TMath::IsNaN(valueF)
          || valueF==99999){
        MSG("NCExtractionMDA",Msg::kDebug) << "Found "
                                           << "default value"
                                           << endl;
        valueTemp(i)=(fMeanVec.at(0))(i);

      } else {

        valueTemp(i)=valueF;

      }//if(defaultF)
    }//if(float)
    else if(strcmp(thisLeaf->GetTypeName(),"Int_t")==0
            || strcmp(thisLeaf->GetTypeName(),"UInt_t")==0
            || strcmp(thisLeaf->GetTypeName(),"int")==0){
      valueI=static_cast<Int_t>(thisLeaf->GetValue());

      if (valueI > 1e+8 || ANtpDefVal::IsDefault(valueI)
          || !TMath::Finite(valueI) || TMath::IsNaN(valueI)
          || valueI==99999){
        MSG("NCExtractionMDA",Msg::kDebug) << "Found default "
                                           << "value"
                                           << endl;
        valueTemp(i)=(fMeanVec.at(0))(i);

      } else {

        valueTemp(i)=static_cast <Double_t> (valueI);

      }//if(default)
    }//if(int)
    MSG("NCExtractionMDA",Msg::kDebug) << "Final value:"
                                       <<thisLeaf->GetName() << ", "
                                       << valueTemp(i) << endl;


  }//for fLeafIndexVec

  varPVector.push_back(valueTemp);
  valueTemp.Clear();
}


//-----------------------------------------------------------------------
double NCExtractionMDA::GetIdProbability(NCEventInfo& evtInfo, int)
{
  //Necessary for data running because Prepare now only runs for MC
  if (!fIsFilled) fIsFilled = ReadCalibInfoFromFile(evtInfo);
  assert(fIsFilled);
  //If event length is larger than 40 planes, classify as CC
  if ( (evtInfo.event->endPlane - evtInfo.event->begPlane+1) > 40 )
    return 1.0;

  //Fill P vector of variables for event
  FillPVector(evtInfo);

  VecDouble_t probClass(fNClass);
  VecDouble_t discrimValue(fNClass);
  //Calculate discriminant function Q for each class
  for (Int_t n=0; n < fNClass; n++){
    //Create an auxiliary copy of the Pvector
    TVectorD varPClone=varPVector.at(0);

    //Note that varPClone is modified by the "*=" operator
    Double_t quadTerm=((varPClone)*=(fQuadCoeff.at(n)))*(varPVector.at(0));

    Double_t linTerm=(fLinCoeff.at(n))*(varPVector.at(0));

    Double_t constTerm=(fConstCoeff.at(n))(0);

    discrimValue.at(n)=quadTerm+linTerm+constTerm;

    MSG("NCExtractionMDA",Msg::kDebug)<< "Q("
                                      << (fTypeClassVec.at(n)).at(0)
                                      <<") = "
                                      << discrimValue.at(n) <<endl;
  }

  TVectorD tempDiff(fNClass-1);
  VecTVectorD discrimDiff;

  //Calculate discriminator difference vectors
  //P(a) = e^a/(e^(a)+e^(b)+e^(c)+...) = 1/(1+e^(b-a)+e^(c-a)+...)
  //Avoids e^X where X is big => machine-dep floating exception.

  for (Int_t n=0; n < fNClass; n++){

    Int_t diffCount=0;

    for (Int_t m=0; m < fNClass; m++){
      if(m!=n){
        tempDiff(diffCount)=(discrimValue.at(m)-discrimValue.at(n));
        diffCount++;
      } else continue;
    }

    if (TMath::Abs(tempDiff.Min())>700. || TMath::Abs(tempDiff.Max())>700){
      MSG("NCExtractionMDA",Msg::kDebug)<<"Bypass floating exception "
                                        << endl;
      break;
    }else discrimDiff.push_back(tempDiff);

    tempDiff.Zero();
  }

  if(discrimDiff.size()!=static_cast<UInt_t>(fNClass)){
    varPVector.clear();
    probClass.clear();

    return ANtpDefVal::kFloat;
  }

  double pidNC = 0;

  //Calculate probability value for each class
  for (Int_t n=0; n < fNClass; n++){

    Double_t denomSum=0.0;

    for (Int_t m=0; m < (fNClass-1); m++){
      denomSum+=TMath::Exp((discrimDiff.at(n))(m));
    }

    probClass.at(n)=(1.0/(1.0+denomSum));

    MSG("NCExtractionMDA",Msg::kDebug) << "P("
                                       << (fTypeClassVec.at(n)).at(0)
                                       <<") = "
                                       << probClass.at(n) <<endl;


    if((fTypeClassVec.at(n)).at(0)=="ncu")
      pidNC=probClass.at(n);

  }

  varPVector.clear();
  probClass.clear();

  MSG("NCExtractionMDA",Msg::kDebug)<<"Calculating Prob(NC) done" << endl;

  //NCUtils assumes NC has low values of separation parameter
  //and CC high values, so this makes it happy.
  if (pidNC < 0.000001) return 1.0;
  else return (1.0-pidNC);

}

TString NCExtractionMDA::WhichMdaCoeffHistDef(Bool_t isCoeff,
                                              Detector::Detector_t detType,
                                              ReleaseType::Release_t mcType)
{
  assert(mcType != ReleaseType::kUnknown);

  assert(!ReleaseType::IsCarrot(mcType));

  if(detType == Detector::kFar){
    if(ReleaseType::IsDaikon(mcType)){
      TString s = isCoeff ? "Coeff" : "HistDef";
      return "Mda_"+s+"_Far_L010z185i_CedarDaikon.dat";
    }
    return "";
  }

  if(detType == Detector::kNear){
    if(ReleaseType::IsDaikon(mcType)){
      TString s = isCoeff ? "Coeff" : "HistDef";
      return "Mda_"+s+"_Near_L010z185i_CedarDaikon.dat";
    }
    return "";
  }

  return "";
}

---