// A program to compute a PDF correlation cosine between two observables X and Y,
// evaluated for N PDF sets (with values for each observable provided in text
// files file1 and file2)
// Author: P. Nadolsky, 2007
// Compilation: g++ -o corr corr.c
// Usage: corr file1 file2

using namespace std;
#include <iostream>
#include <fstream>
#include <math.h>


//--------------------------------------------------------------------------
int cntnbk(const char *str)
/*Counts non-blank charachters in a string str until a first blank character
  or the end of the string is met*/
{int m=0;
 char tem;
 while (str[m] !=' ' && (tem=str[m++],tem)) ;

 if (str[m-1]=='\0')
   return m-1;  //m non-blank symbols counted
 else
   return m; //we do not count the end-of-string symbol
}

//----------------------------------------------------------------------

int main(int argc, char *argv[])
{
 char fname1[40], fname2[40];
 const int mpdf=100;
 int npdf=0, ipdf1=0,ipdf2=0;
 float xsec1[mpdf], xsec2[mpdf], err1, err2, corr;


  if (argc < 2 )
    { 
      cout << "corr: Computation of PDF errors and correlation cosines for "
	   <<endl << "two physical observables X and Y" << endl;
      cout << "Author: P. Nadolsky, 2007" << endl;
      cout << endl<< "Usage: corr file1 file2"<<endl<<
	"where file1 and file2 are text files with values of X and Y"<<endl
	   <<"for N_PDF eigenvector sets (with the first values given for" << endl
	   << "the central PDF set)"<<endl;
      cout << endl << "When publishing results obtained with this program, please" << endl
      << "acknowledge its source, " << endl << 
	"\"Implications of CTEQ global analysis for collider observables\"" <<endl
	   <<"by P. M. Nadolsky, H.-L. Lai, Q.-H. Cao, J. Huston, J. Pumplin, D. Stump," << endl
	   << "W.-K. Tung, and C.-P. Yuan, published in Phys. Rev. D78, 013004 (2008)."<<endl; 
      cout << endl << "STOP: too few arguments" << endl;
      exit(1); 
    }

  strcpy(fname1,argv[1]);
  strcpy(fname2,argv[2]);

  ifstream fin1(argv[1]), fin2(argv[2]);
  
  if (fin1==NULL || fin2==NULL)
    {
      cout << "STOP: error reading files" << endl; 
      return 1;
    }
  
  for (ipdf1=0; ipdf1<=mpdf-1; ipdf1++)
    xsec1[ipdf1]=xsec2[ipdf1]=0.0;


  ipdf1=0;
  while (fin1.good())
    {
      fin1>>xsec1[ipdf1];
      ipdf1++;
    }

  ipdf2=0;
  while (fin2.good())
    {
      fin2>> dec>>xsec2[ipdf2];
      ipdf2++;
    }

  //Drop trailing zeros at the end of cross sections
  while(xsec1[ipdf1-1]==0.0) ipdf1--;
  while(xsec2[ipdf2-1]==0.0) ipdf2--;

  if (ipdf1 == ipdf2)
    {
      cout << "Total number of extreme PDF sets: "<<ipdf1<<endl;
      npdf=ipdf1-1;
    }
  else
    {
      cout << "STOP: numbers of PDF sets do not match: "
	   <<ipdf1<<" and "<<ipdf2<<endl;
      return 1;
    }

  err1=err2=corr=0.0;

  for (ipdf1=1; ipdf1 <=npdf/2; ipdf1++)
    {
      err1 += pow( xsec1[2*ipdf1] - xsec1[2*ipdf1-1] ,2);
      err2 += pow( xsec2[2*ipdf1] - xsec2[2*ipdf1-1] ,2);
      corr += (xsec1[2*ipdf1] - xsec1[2*ipdf1-1])*
	      (xsec2[2*ipdf1] - xsec2[2*ipdf1-1]);
    }

  err1 = 0.5*sqrt(err1); err2=0.5*sqrt(err2); 
  corr = 0.25*corr/err1/err2;
  
  cout<<"   File1     File2   X   dX     Y    dY  cos(phi(X,Y)):"<<endl;
  cout << fname1<<" "<<fname2<<" "<<xsec1[0]<<" "<<err1<<" "
       <<xsec2[0]<<" "<<err2<<" "<<corr<<endl;

  fin1.close(); fin2.close();

  return 0;
} /*main -> */