FRIB_e21062/armory/AnalysisLibrary.h

#ifndef Analysis_Library_h
#define Analysis_Library_h

#include <TF1.h>
#include <TGraph.h>
#include <TSpectrum.h>
#include <TMath.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <string>

std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){

  std::vector<std::string> output;

  size_t pos;
  do{
    pos = tempLine.find(splitter); // fine splitter
    if( pos == 0 ){ //check if it is splitter again
      tempLine = tempLine.substr(pos+1);
      continue;
    }

    std::string secStr;
    if( pos == std::string::npos ){
      secStr = tempLine;
    }else{
      secStr = tempLine.substr(0, pos+shift);
      tempLine = tempLine.substr(pos+shift);
    }

    //check if secStr is begin with space
    while( secStr.substr(0, 1) == " "){
      secStr = secStr.substr(1);
    };
    
    //check if secStr is end with space
    while( secStr.back() == ' '){
      secStr = secStr.substr(0, secStr.size()-1);
    }

    output.push_back(secStr);
    //printf(" |%s---\n", secStr.c_str());
    
  }while(pos != std::string::npos );

  return output;
}

std::vector<std::vector<double>> combination(std::vector<double> arr, int r){
  
  std::vector<std::vector<double>> output;
  
  int n = arr.size();
  std::vector<int> v(n);
  std::fill(v.begin(), v.begin()+r, 1);
  do {
    //for( int i = 0; i < n; i++) { printf("%d ", v[i]); }; printf("\n");
    
    std::vector<double> temp;
    for (int i = 0; i < n; ++i) { 
      if (v[i]) {
        //printf("%.1f, ", arr[i]); 
        temp.push_back(arr[i]);
      }
    }
    //printf("\n");
    
    output.push_back(temp);
    
  } while (std::prev_permutation(v.begin(), v.end()));
  
  return output;
}

double* sumMeanVar(std::vector<double> data){
  
  int n = data.size();
  double sum = 0;
  for( int k = 0; k < n; k++) sum += data[k];
  double mean = sum/n;
  double var = 0;
  for( int k = 0; k < n; k++) var += pow(data[k] - mean,2);
  
  static double output[3];
  output[0] = sum;
  output[1] = mean;
  output[2] = var;
  
  return output;
}

double Rsquared(std::vector<double> x, std::vector<double> y){
  
  if ( x.size() != y.size() ) return 0;
  
  int n = x.size();
  double  * qX = sumMeanVar(x);
  
  double meanX = qX[1];
  double varX = qX[2];
  
  double  * qY = sumMeanVar(y);
  double meanY = qY[1];
  double varY = qY[2];
  
  double sumXY = 0;
  for( int k = 0; k < n; k++) sumXY += x[k] * y[k];

  double r2 = (sumXY - n * meanX * meanY)/ sqrt( varX * varY );
  
  return r2;
} 

double*  fitSlopeIntercept(std::vector<double> dataX, std::vector<double> dataY){
  
  double * smvY = sumMeanVar(dataY);
  double sumY = smvY[0];
  double meanY = smvY[1];

  double * smvX = sumMeanVar(dataX);
  double sumX = smvX[0];
  double meanX = smvX[1];
  double varX = smvX[2];
  
  int n = dataX.size();
  double sumXY = 0;
  for( int j = 0; j < n; j++) sumXY += dataX[j] * dataY[j];

  double slope = ( sumXY - sumX * sumY/n ) / varX;
  double intercept = meanY - slope * meanX;
  
  static double output[2];
  output[0] = slope;
  output[1] = intercept;
  
  return output;
  
}



std::vector<std::vector<double>> FindMatchingPair(std::vector<double> enX, std::vector<double> enY){

   //output[0] = fitEnergy;
   //output[1] = refEnergy;

   int nX = enX.size();
   int nY = enY.size();
   
   std::vector<double> fitEnergy;
   std::vector<double> refEnergy;
   
   if( nX > nY ){
      
      std::vector<std::vector<double>> output = combination(enX, nY);
      
      double * smvY = sumMeanVar(enY);
      double sumY = smvY[0];
      double meanY = smvY[1];
      double varY = smvY[2];
      
      double optRSquared = 0;
      double absRSqMinusOne = 1;
      int maxID = 0;
      
      for( int k = 0; k < (int) output.size(); k++){
        
        double * smvX = sumMeanVar(output[k]);
        double sumX = smvX[0];
        double meanX = smvX[1];
        double varX = smvX[2];
        
        double sumXY = 0;
        for( int j = 0; j < nY; j++) sumXY += output[k][j] * enY[j];
        
        double rSq = abs(sumXY - sumX*sumY/nY)/sqrt(varX*varY);
        
        //for( int j = 0; j < nY ; j++){ printf("%.1f, ", output[k][j]); }; printf("| %.10f\n", rSq);
        
        if( abs(rSq-1) < absRSqMinusOne ) {
          absRSqMinusOne = abs(rSq-1);
          optRSquared = rSq;
          maxID = k;
        }
      }
      
      fitEnergy = output[maxID];
      refEnergy = enY;
      
      printf(" R^2 : %.20f\n", optRSquared);      
      
      //calculation fitting coefficient
      //double * si = fitSlopeIntercept(fitEnergy, refEnergy);
      //printf( " y = %.4f x + %.4f\n", si[0], si[1]);
      
    }else if( nX < nY ){
    
      std::vector<std::vector<double>> output = combination(enY, nX);
      
      
      double * smvX = sumMeanVar(enX);
      double sumX = smvX[0];
      double meanX = smvX[1];
      double varX = smvX[2];
      
      double optRSquared = 0;
      double absRSqMinusOne = 1;
      int maxID = 0;
      
      for( int k = 0; k < (int) output.size(); k++){
        
        double * smvY = sumMeanVar(output[k]);
        double sumY = smvY[0];
        double meanY = smvY[1];
        double varY = smvY[2];
        
        double sumXY = 0;
        for( int j = 0; j < nX; j++) sumXY += output[k][j] * enX[j];
        
        double rSq = abs(sumXY - sumX*sumY/nX)/sqrt(varX*varY);
        
        //for( int j = 0; j < nX ; j++){ printf("%.1f, ", output[k][j]); }; printf("| %.10f\n", rSq);
        
        if( abs(rSq-1) < absRSqMinusOne ) {
          absRSqMinusOne = abs(rSq-1);
          optRSquared = rSq;
          maxID = k;
        }
      }
      
      fitEnergy = enX;
      refEnergy = output[maxID];
      printf(" R^2 : %.20f\n", optRSquared);   
    
    }else{
      fitEnergy = enX;
      refEnergy = enY;
      
      //if nX == nY, ther could be cases that only partial enX and enY are matched. 
      
    }
    
    
    printf("fitEnergy = ");for( int k = 0; k < (int) fitEnergy.size() ; k++){ printf("%7.2f, ", fitEnergy[k]); }; printf("\n");
    printf("refEnergy = ");for( int k = 0; k < (int) refEnergy.size() ; k++){ printf("%7.2f, ", refEnergy[k]); }; printf("\n");
    
    std::vector<std::vector<double>> haha;
    haha.push_back(fitEnergy);
    haha.push_back(refEnergy);
    
    return haha;
   
}


std::vector<std::vector<double>> LoadCorrectionParameters(TString corrFile, bool show=false){

  printf("  load correction parameters : %s", corrFile.Data());
  std::ifstream file;
  file.open(corrFile.Data());
  
  std::vector<std::vector<double>> corr;
  corr.clear();
  
  std::vector<double> detCorr;
  detCorr.clear();
  
  if( file.is_open() ){  
    
    while( file.good() ){
    
      std::string line;
      getline(file, line);
    
      if( line.substr(0,1) == "#" ) continue;
      if( line.substr(0,2) == "//" ) continue;
      if( line.size() == 0 ) continue;

      //printf("%s\n", line.c_str());
      std::vector<std::string> temp = SplitStr(line, " ");
      
      detCorr.clear();
      for( int i = 0; i < (int) temp.size() ; i++){
        detCorr.push_back(std::stod(temp[i]));
      }
      corr.push_back(detCorr);
    }
    
    file.close();
    
    printf(".... done\n");
    if( show ){
      printf("===== correction parameters \n");
      for( int i = 0; i < (int) corr.size(); i++){
        printf("row : %2d | ", i );
        int len = (int) corr[i].size();
        for( int j = 0; j < len - 1 ; j++){
          printf("%14.6f, ", corr[i][j]);
        }
        printf("%14.6f\n", corr[i][len-1]);
      }
    }
  
  }else{
       printf(".... fail\n");
       std::vector<double> temp = {0, 1};
       for( int i = 0; i < 36; i++){
         corr.push_back(temp);
      }
  }
  
  return corr;
}

#endif