/*
  Fit function shown in the paper

  DESY-20-080

  Measurement of exclusive pi+pi- and rho0 meson photoproduction at HERA

  This is example code how to parse and read the "stat" and "syst"
  files with the full uncertainty breakdown corresponding to table 8
  (the fit depicted in figure 11)

  The code is optimized for execution with the Root analysis framework.

  Copyright: 2020 for the benefit of the H1 collaboration

  When using this code, please include a proper citation 
  to our work in your publication.

*/

#include


void d20_080t8_Wfit(char const *name="d20-080t8_syst.csv",
		char const *nameStat="d20-080t8_statCorr.csv") {
  // systematic (offset variations)
  vector<vector<double> > syserrUp,syserrDown;
  // parameter and stat error
  vector<double> par,dparStat;
  // papameter names
  vector<TString> nameUp,nameDown;

  // loop over "syst" file
  ifstream infile(name);
  while(!infile.eof()) {
    char line[10000];
    infile.getline(line,9999);
    line[9999]=0;
    // skip comment lines
    if(line[0]=='#') continue;

    // replace ",' by blanks for easy parsing
    for(int i=0;line[i];i++) {
      if(line[i]==',') line[i]=32;
    }
    // parse a single line, arguments separated by blanks
    istringstream sline(line);
    string name;
    // first entry is not used (parameter number)
    sline>>name;
    // second entry: parameter name
    sline>>name;
    if(name=="parName") {
       // header line of the CSV
       // this is used to determine the number of systematic variations
       string chan,val,stat,tpos,tneg;
       sline>>chan>>val>>stat>>tpos>>tneg;
       // this loop reads all systematic parameter names and counts them
       while(!sline.eof()) {
          string up,down;
          sline>>up>>down;
          if(sline.fail()) continue;
          nameUp.push_back(up);
          nameDown.push_back(down);
       }
       // resize array of systematic shifts
       syserrUp.resize(nameUp.size());
       syserrDown.resize(nameDown.size());
    } else {
       // data line of the CSV
       // channel (for table 8 always elastic)
      string chan;
      sline>>chan;
      // read central value, stat error, total positive and negative error
      double x,stat=0.,up=0.,dn=0.;
      sline>>x>>stat>>up>>dn;
      if(!sline.fail()) {
         // parameter is valid -> store it
	par.push_back(x);
	dparStat.push_back(stat);
	int ierr=0;
        // in this loop, read all systematic up and down shifts
        double sysUp=0.,sysDown=0.;
	while(!sline.eof()) {
	  double eup=0.,edn=0.;
          // shift resultig from up and down variation
	  sline>>eup>>edn;
	  if(sline.fail()) break;
	  syserrUp[ierr].push_back(eup);
	  syserrDown[ierr++].push_back(edn);
          // the following is for determining the total up and down
          //  systematic uncertainty (not used for plotting, only for
          //  consistency check)
	  if(eup<edn) {
	    double t=eup; eup=edn;edn=t;
	  }
          // the if clauses reject cases where both the up and down variation
          // result in shifts in ths same direction
	  if(eup>0.)
	    sysUp =TMath::Hypot(sysUp,eup);
	  if(edn<0.)
	    sysDown =TMath::Hypot(sysDown,edn);
	}
	double totUp =TMath::Hypot(sysUp,stat);
	double totDown =TMath::Hypot(sysDown,stat);
	// verify total uncertainty in table
	cout<<name<<" "<<x<<" +/-"<<stat<<" (stat)  "
            <<"+"<<up<<"-"<<dn<<" (tot,fromFile)  "
            <<"+"<<sysUp<<"-"<<sysDown<<" (syst)  "
            <<"+"<<totUp<<"-"<<totDown<<" (tot,calculated)\n";
      }
    }
  }
  int nsys=syserrUp.size();
  int npar=syserrUp[nsys-1].size();
  //cout<<"number of sys errors: "<<nsys<<" parameters: "<<npar<<"\n";

  // read stat correlations (ststCorr file)
  ifstream infile2(nameStat);
  TMatrixDSym rhoij(npar),Vstat(npar);
  while(!infile2.eof()) {
     // read one line
     char line[10000];
     infile2.getline(line,9999);
     line[9999]=0;
     // ship comment lines
     if(line[0]=='#') continue;
     // replace comma by blank for simple parsing
     for(int i=0;line[i];i++) {
        if(line[i]==',') line[i]=32;
     }
     istringstream sline(line);
     int ipar,jpar;
     string iname,jname;
     // try to parse all lines with their numeric values
     // this will automatically skip the CSV header line
     sline>>ipar>>jpar>>iname>>jname;
     if(!sline.fail()) {
        string dummy;
        sline>>dummy;
        sline>>dummy;
        double rho;
        sline>>rho;
        if(sline.fail()) continue;
        // only one half f teh matrix is stored, so set both Vij and Vji
        rhoij(ipar,jpar)=rho;
        rhoij(jpar,ipar)=rho;
        Vstat(ipar,jpar)=rho*dparStat[ipar]*dparStat[jpar];
        Vstat(jpar,ipar)=rho*dparStat[ipar]*dparStat[jpar];
     }
  }

  // determine independent eigenvectors, these can be added quadratically later
  TMatrixDSymEigen EStat(Vstat);

  vector<double> Wi,fi,fstatP,fstatM,ftotP,ftotM;

  // loop over the W range to plot
  // logarithmic spacing
  for(double W=2.;W<300.;W*=1.2) {
    // offset variation of stat eigenvalues and sys errors
    for(int ivar=0;ivar<1+2*par.size()
	  +syserrUp.size()+syserrDown.size();ivar++) {
      // ivar=0: central value
      double A=par[0];
      double fR=par[1];
      double deltaR=par[2];
      double deltaP=par[3];
      if(ivar>0.) {
         // ivar>0: statistical of systematic variation
	if(ivar<=2*par.size()) {
           // statistical eigenvector variation
	  TVectorD eval=EStat.GetEigenValues();
	  TMatrixD evec=EStat.GetEigenVectors();
	  int k=(ivar-1)%(2*par.size());
	  // iv: variant number 
	  int iv=k%par.size();
	  double d=sqrt(eval[iv]);
          // positive variations: k=0..N-1
          // negative variations: k=N..2*N-1
	  if(k>=par.size()) d=-d;
          // modify central value
	  A += d*evec(0,iv);
	  fR += d*evec(1,iv);
	  deltaR += d*evec(2,iv);
	  deltaP += d*evec(3,iv);
	} else {
	  // systematic offset up/down variations
	  int k=ivar-1-2*par.size();
	  if(k<syserrUp.size()) {
             // k=0..Nsys-1
             //  Up variations

             // modify central value
             A += syserrUp[k][0];
             fR += syserrUp[k][1];
             deltaR += syserrUp[k][2];
             deltaP += syserrUp[k][3];
	  } else {
             // k=Nsys..2*Nsys-1
             // down variations

             // bring k to range 0..Nsys-1
             k-=syserrUp.size();

             // modify central value
             A += syserrDown[k][0];
             fR += syserrDown[k][1];
             deltaR += syserrDown[k][2];
             deltaP += syserrDown[k][3];
	  }
	}
      }
      // parametrisations
      double W0=40;
      double f=A*(fR*pow(W/W0,-deltaR)+pow(W/W0,deltaP));

      if(ivar==0) {
         // central value
         // store Wi[] and fi[]
         //tes up and down envelopes to zero
         Wi.push_back(W);
         fi.push_back(f);
         fstatP.push_back(0.);
         fstatM.push_back(0.);
         ftotP.push_back(0.);
         ftotM.push_back(0.);
      } else {
         // error contribution, add quadratically, depending on size
         int ii=fi.size()-1;
         double df=f-fi[ii];
         if(ivar<=2*par.size()) {
            // statistical error band
            if(df<0.) {
               fstatM[ii]= -hypot(fstatM[ii],df);
            } else {
               fstatP[ii]= hypot(fstatP[ii],df);
            }
         }
         // total error band
         if(df<0.) {
            ftotM[ii]= -hypot(ftotM[ii],df);
         } else {
            ftotP[ii]= hypot(ftotP[ii],df);
         }
      }
    }
  }
  // print calculated Wi[] and predicted cross section fi[]
  // and size of statistical and total envelope
  for(int i=0;i<fi.size();i++) {
    cout<<Wi[i]<<" "<<fi[i]<<" "<<fstatM[i]<<" "<<fstatP[i]
	<<" "<<ftotM[i]<<" "<<ftotP[i]<<"\n";
  }
}