/* Flow analysis through the cumulant method. 
   Version with multiparticle correlations. */

/* The input is a file 'events.dat' 
   following the structure of the file generated by 'detector.cc', i.e.: 
   - first, the number of events; 
   - then, for each event, the (detected) particle multiplicity, 
   and for each particle, the transverse momentum and the azimuth. */

/* The program then reconstructs integrated elliptic flow (V2) and 
   both differential elliptic (v2) and quadrupolar flows.
   Numbered equations refer to the paper 
   [A]: N.Borghini, P.M.Dinh, J.-Y.Ollitrault, Phys. Rev. C64 (2001) 054901
        [=nucl-th/0105040]. */

/* This program will give you first the INTEGRATED elliptic flow V2 
   reconstructed using the lowest order formulas of the paper, 
   up to order 6 (section II of paper A):
    . order 2 (in the paper: v2{2}) corresponds to the "standard" flow analysis
     (two-particle nonflow correlations neglected);
    . order 4 (v2{4}) eliminates two-particle nonflow correlations;
    . order 6 (v2{6}) eliminates also four-particle nonflow correlations. */

/* Once integrated flow is reconstructed, the program evaluates 
   the DIFFERENTIAL flow v'_{nm} of the particles
   using the lowest order formulas of the paper (section III):
   v'_nm{m+1} is the standard flow analysis (nonflow correlations neglected);
   v'_nm{m+3} is the next order. 
   m=1 yields elliptic flow v'2, m=2 yields v'4. */


/* The results are stored in the output file 'cumulants.txt' */

#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <sys/time.h>
#include "distributions.h"

/* */

#define Pi (2.*acos(0.))
#define rad2deg (180./Pi)
#define sqr(x) ((x)*(x))


/* PARAMETERS USED IN THE RECONSTRUCTION */
/* kr*kt is the number of points at which one computes the 
   generating function [cf Eq.(B1)]
   kr possible values for the radius 
   kt possible values for the angle: should satisfy kt>2*kr */
#define kr 3
#define kt 7

/* zmax is the maximum value of |z| for which you tabulate the generating 
   functions G_n and C_mn, defined by Eqs.(5) and (7). 
   If zmax is too small, then G_n and C_mn are too close to 1, 
   and you may have large numerical errors when computing, 
   especially if you do not use double precision. 
   If zmax is too large, then the power series expansion of C_mn, 
   which is truncated to order 2 in this program, is not accurate. 
   Try several values of zmax, to check the accuracy of the interpolation 
   procedure. */
#define zmax 1.5
#define r0 (zmax/sqrt(kr))

#define Nbins 60 /* number of bins for differential flow.
  Check that it is consistent with the number in distributions.cc! */

/* */

FILE *input,*output,*plotvdiff0,*plotvdiff1;

double prod1(complex<double> z, double phi, double w)
{
 double temp;
  
 temp=1.+2.*w*(real(z)*cos(phi)+imag(z)*sin(phi));
 return temp;
}

main()
{
  const complex<double> i(0,1); 

  /* Variable declaration */
  long neve; /* Number of events */
  long rmult;  /* Detected multiplicity in an event */
  long rpmult[Nbins]; /* Detected multiplicity per bin */
  double rmultmean=0.; /* Average detected multiplicity */
  float phi, pT; /* particle azimuth and transverse momentum */
  int binpT;

  /* Interpolation points [see Eq. (I-B1)]: z[k1][k2] = rz[k1] e^(i arg[k2]).*/
  complex<double> z[kr][kt]; 
  double rz[kr], arg[kt];

  /* Generating functions: 
     - G2 generating function for integrated flow; 
     - G2mean is the generating function G2 averaged over events;
     - g22, generating function for v'2 with respect to v2.
     - g42, generating function for v'4 with respect to v2. */
  double G2[kr][kt], G2mean[kr][kt];
  complex <double> g22[Nbins][kr][kt], g42[Nbins][kr][kt];

  double w2[Nbins]; /* particle weight */
  complex <double> sum22[Nbins][kr][kt], sum42[Nbins][kr][kt];

  /* Cnrad is the generating function of cumulants for integrated flow, 
     averaged over angles (Eq.(B3)); 
     Dmrad (m=1 or 2) is the generating function of cumulants for differential 
     flow, averaged over angles (Eq.(B7)).
     cn and dmn are the cumulants used for the reconstruction of integrated 
     and differential flows resp. (c_n{.} and d_mn/n{.}) */
  double C2rad[kr], c2[kr];
  double D22rad[Nbins][kr], D42rad[Nbins][kr], d22[Nbins][kr], d42[Nbins][kr];

  int mflag2[3]; /* tells whether the reconstruction to order k is applicable
	            or not. */

  /* V2[k] are the integrated flow estimates obtained from the cumulant to 
     order 2(k+1), i.e. v_2{2k+2} in the Ref. */
  double V2[kr]; 
  double chi2sq; 

  /* v22[], v42[], are the differential flow estimates in each bin.
     v22[0] is v'2/2{2} (Ref.[A]), computed with V2[0]: 
     this is a TWO-particle estimates; 
     v22[1] is v'2/2{4}, computed with V2[1]: this is a FOUR-particle estimate;
     v42[0], v42[1] are v'4/2{3}, v'4/2{5}, computed with V2[1]: 
     these are MULTIparticle estimates.  */
  double v22[Nbins][2], v42[Nbins][2]; 
  /* err22[], err42[] are the respective statistical error bars in each bin. */
  double err22[Nbins][2], err42[Nbins][2]; 
  double vd22av[2], vd42av[2];

  double temp, err; 

  struct timeval Time;
  int debut, fin;

  /* Initializations */
  gettimeofday (&Time, 0);
  debut=Time.tv_sec;

  for(int k1=0; k1<kr; k1++)
  {
    rz[k1]=r0*sqrt(k1+1.);
    C2rad[k1]=0.;
    for(int bin=0; bin<Nbins; bin++) D22rad[bin][k1]=D42rad[bin][k1]=0.;
    for(int k2=0; k2<kt; k2++)
    {
      /* z[][] are the coordinates of the points where you need to 
         compute the generating functions G and Dm, see Eq.(B1) */
      z[k1][k2]=rz[k1]*exp(i*(arg[k2]=2.*k2*Pi/kt));

      /* Gmean is the generating function averaged over events: <G> */
      G2mean[k1][k2]=0.;
      /* gmn is the generating function gm=exp(im psi) G_n */
      for(int bin=0; bin<Nbins; bin++) g22[bin][k1][k2]=g42[bin][k1][k2]=0.;
    }
  }

  for(int bin=0; bin<Nbins; bin++) 
  {
    rpmult[bin]=0;
    w2[bin]=1.; /* unit weights for the moment... */
  }

  /* Scan of the data */
  printf("Reading the data file\n");
  printf("For each event, building the generating function\n");
  input=fopen("events.dat", "r");
  fscanf(input, "%i\n", &neve); /* Read the number of events */
  printf("----------\n");
  printf("Number of events in the run: %i\n", neve);

  output=fopen("cumulants.txt", "w");
  fprintf(output, "CUMULANT ANALYSIS OF THE FLOW\n\n");
  fprintf(output, "version with the generating function of correlations\n\n"); 
  fprintf(output, "Available statistics:\n  number of events Neve=%i\n", neve);

  for(long ne=0; ne<neve; ne++) /* Loop over events */
  {
    if (floor((ne+1.)*10./neve)==((ne+1.)*10./neve))
      printf("%.f percent of events processed\n", (ne+1.)*100./neve);

    for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++)
      for(int bin=0; bin<Nbins; bin++) sum22[bin][k1][k2]=sum42[bin][k1][k2]=0.;

    fscanf(input, "%i", &rmult); /* Read the event multiplicity */
    rmultmean+=((double) rmult); /* and compute the average multiplicity. */

    /* RECONSTRUCTION OF THE FLOW FROM THE GENERATING FUNCTION */
    /* Initialization of the generating function G2 */
    for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++) G2[k1][k2]=1.;

    for(long np=0; np<rmult; np++) /* Loop over the event particles */
    {
      fscanf(input, "%f %f", &pT, &phi); 
      /* Read the particle momentum & azimuth. */
      binpT=pTbin(pT);
      rpmult[binpT]+=1; /* Increment the particle number in the relevant bin */

      for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++)
      {
        /* Compute the generating functions G2 (see Eq.(5)) at the points z[][],
           with (default) weight = 1/M. */
        G2[k1][k2]*=prod1(z[k1][k2], 2.*phi, w2[binpT]/rmult); /* Note the 2. */

        sum22[binpT][k1][k2]+=
          exp(2.*i*((double) phi))/prod1(z[k1][k2], 2.*phi, w2[binpT]/rmult);
        sum42[binpT][k1][k2]+=
          exp(4.*i*((double) phi))/prod1(z[k1][k2], 2.*phi, w2[binpT]/rmult);
      }
    } /* End of the loop over particles */

    /* Sum over all events of G at the points (xz,yz) */
    for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++) 
      G2mean[k1][k2]+=G2[k1][k2];
     
    for(int bin=0; bin<Nbins; bin++) /* Loop over differential bins */
    {
      for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++)
      {
        /* Differential v'2 with respect to integrated v2 */
        g22[bin][k1][k2]+=sum22[bin][k1][k2]*G2[k1][k2]; 

        /* Differential v'4 with respect to integrated v2 */
        g42[bin][k1][k2]+=sum42[bin][k1][k2]*G2[k1][k2]; 
      }
    } /* End of the loop over bins */
  } /* End of the loop over events */
  fclose(input);

  /* Average over all events of G at the interpolation points */
    for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++) G2mean[k1][k2]/=neve; 

  rmultmean/=neve; /* Average detected particle multiplicity */

  gettimeofday (&Time, 0);
  fin=Time.tv_sec;
  printf("All events read, elapsed time=%d s\n", fin-debut);
  printf("Mean detected multiplicity <M>=%G\n", rmultmean);
  fprintf(output,"  mean detected multiplicity <M>=%G\n\n", rmultmean);

  /* Average over all differential particles of gm (m=1 or 2) */ 
  for(int k1=0; k1<kr; k1++) for(int k2=0; k2<kt; k2++)
    for(int bin=0; bin<Nbins; bin++)
    {
      g22[bin][k1][k2]/=rpmult[bin]; 
      g42[bin][k1][k2]/=rpmult[bin]; 
    }

  /* GENERATING FUNCTIONS OF CUMULANTS
     For integrated flow, compute C = M(<G>^(1/M) - 1), 
     then average over angles (see Eq.(B3)): 
     the result 'Crad' depends on the radius only. 

     For differential flow, compute Dm = <gm> / <G>, multiply by z^m, 
     take the real part, and average over angles (see Eq.(B7)): 
     the result 'Dmrad' depends on the radius only. */
    
  for(int k1=0; k1<kr; k1++)
  {
    for(int k2=0; k2<kt; k2++)
    {
      C2rad[k1]+=rmultmean*(pow(G2mean[k1][k2], 1./rmultmean)-1.);
      for(int bin=0; bin<Nbins; bin++)
      {
        D22rad[bin][k1]+=real(exp(-i*arg[k2])*g22[bin][k1][k2])/G2mean[k1][k2];
        D42rad[bin][k1]+=real(exp(-2.*i*arg[k2])*g42[bin][k1][k2])/
			  G2mean[k1][k2];
      }
    }
    C2rad[k1]/=kt;
    for(int bin=0; bin<Nbins; bin++)
    {
      D22rad[bin][k1] = D22rad[bin][k1]*rz[k1]/kt;
      D42rad[bin][k1] = D42rad[bin][k1]*sqr(rz[k1])/kt;
    }
  }

  /* RECONSTRUCTION OF INTEGRATED FLOW 

     One computes the cumulants, which are related to the
     coefficients in the power series expansion of the generating
     function C = M(G^(1/M) - 1), see Eq.(7). 
     These coefficients are obtained from the interpolation formulas, 
     Eqs.[B5]. */
  
  /* Elliptic flow v2 */
  c2[0]=(3.*C2rad[0]-3.*C2rad[1]/2.+C2rad[2]/3.)/sqr(r0);
  c2[1]=2.*(-5.*C2rad[0]+4.*C2rad[1]-C2rad[2])/pow(r0,4.);
  c2[2]=6.*(3.*C2rad[0]-3.*C2rad[1]+C2rad[2])/pow(r0,6.);

  /* uncorrected v2 using order k=2: Eq.(17a) */ 
  if (c2[0]<0.) V2[0]=0.; 
  else 
  {
    V2[0]=sqrt(c2[0]);
    mflag2[0]=1;
  }
  /* uncorrected v2 using order k=4: Eq.(17b) */
  if (c2[1]>0.) V2[1]=0.; 
  else 
  {
    V2[1]=pow(-c2[1],0.25); 
    mflag2[1]=1;
  }
  /* uncorrected v2 using order k=6: Eq.(17c) */
  if (c2[2]<0.) V2[2]=0.; 
  else 
  {
    V2[2]=pow(c2[2]/4.,1./6.); 
    mflag2[2]=1;
  }

  fprintf(output,"\n\nINTEGRATED ELLIPTIC FLOW v2\n");
  fprintf(output,"\tuncorrected\tcorrected\t+(stat)\t\t-(stat)\n");

  /*  Result of the method to order k=2 */
  fprintf(output, " v2{2}:\t");
  chi2sq=rmultmean*sqr(V2[0]);
  err=sqrt((1.+2.*chi2sq)/neve)/rmultmean;
  fprintf(output, "%.4g\t\t%.4g\t\t+%.4g", V2[0], V2[0],
	  sqrt(sqr(V2[0])+err)-V2[0]);
  if (err<sqr(V2[0]))
    fprintf(output, "\t-%.4g\n", V2[0]-sqrt(sqr(V2[0])-err));
  else
    fprintf(output, "\t-%.4g\n", V2[0]);

  /* Result of the method to order k=4 */
  fprintf(output, " v2{4}:\t");
  chi2sq=rmultmean*sqr(V2[1]);
  err=sqrt(4.*(1.+chi2sq*(4.+chi2sq*(1.+2.*chi2sq)))/neve)/sqr(rmultmean);
  fprintf(output, "%.4g\t\t%.4g\t\t+%.4g", V2[1], V2[1], 
	  pow(pow(V2[1],4.)+err, 0.25)-V2[1]);
  if (err<pow(V2[1],4.))
    fprintf(output, "\t-%.4g\n", V2[1]-pow(pow(V2[1],4.)-err,0.25));
  else
    fprintf(output, "\t-%.4g\n", V2[1]);

  /* Result of the method to order k=6 */
  fprintf(output, " v2{6}:\t");
  chi2sq=rmultmean*sqr(V2[2]);
  err=
    sqrt(12.*(3.+chi2sq*(18.+chi2sq*(9.+chi2sq*(28.+chi2sq*(12.+24.*chi2sq)))))
	   /(4.*neve))/pow(rmultmean,3.);
  fprintf(output, "%.4g\t\t%.4g\t\t+%.4g", V2[2], V2[2],
	  pow(pow(V2[2],6.)+err, 1./6.)-V2[2]);
  if (err<pow(V2[1],6.))
    fprintf(output, "\t-%.4g\n", V2[2]-pow(pow(V2[2],6.)-err, 1./6.));
  else
    fprintf(output, "\t-%.4g\n", V2[2]);

  /* RECONSTRUCTION OF DIFFERENTIAL FLOW
     
     One computes the cumulants, which are related to the coefficients 
     in the power series expansion of the generating function Dm = gm / G, 
     see Eqs.(8) & (B8). 
     These coefficients are obtained by the interpolation formulas
     Eqs.(B9)-(B10). 
   
     The differential flow is then obtained from the cumulants 
     using Eqs.(34)-(35)]. */

  for(int bin=0; bin<Nbins; bin++)
  { /* Loop over differential bins */
    /* Elliptic flow v'2 */
    /* v'2 with respect to v2: Eqs.(B9) */
    d22[bin][0]=(2.*D22rad[bin][0]-D22rad[bin][1]/2.)/sqr(r0);
    d22[bin][1]=(-2.*D22rad[bin][0]+D22rad[bin][1])/pow(r0,4);

    if (mflag2[0]==1) /* if V2[0] does not vanish */
      v22[bin][0]=d22[bin][0]/V2[0]; /* v'2/2{2}, using V2{2}, Eq.(34a) */ 
        /* uncorrected v'2 using highest order, v'2/2{4}, Eq.(34b) */       
    else
      v22[bin][0]=0; 

    chi2sq=rmultmean*sqr(V2[0]);
    err22[bin][0]=sqrt((1.+chi2sq)/(2.*chi2sq*rpmult[bin])); /* Eq.(D35a) */

    if (mflag2[1]==1) /* V2[1] nonvanishing */
      v22[bin][1]=-d22[bin][1]/pow(V2[1],3);  
      /* v'2{4}, using V2{4}, see Eq.(34b) */ 
    else 
        v22[bin][1]=0.;

    chi2sq=rmultmean*sqr(V2[1]);
    err22[bin][1]=sqrt((2.+chi2sq*(6.+chi2sq*(1.+chi2sq)))/
		       (2.*rpmult[bin]*pow(chi2sq, 3.))); /* Eq.(D35c) */

    /* v'4 with respect to v2: Eqs.(B10) */
    d42[bin][0]=(4.*D42rad[bin][0]-D42rad[bin][1]/2.)/pow(r0, 4);
    d42[bin][1]=(-6.*D42rad[bin][0]+1.5*D42rad[bin][1])/pow(r0, 6);

    if (mflag2[1]==1) /* if V2[1] does not vanish */
    {
      v42[bin][0]=d42[bin][0]/sqr(V2[1]); 
      /* v'4/2{3}, using V2{4}, see Eq.(35a) */ 
      v42[bin][1]=-0.5*d42[bin][1]/pow(V2[1], 4);
      /* v'4/2{5}, using V2{4} as well, see Eq.(35b) */       
      }
      else /* if V2[1]=0 */
        v42[bin][0]=v42[bin][1]=0;

      chi2sq=rmultmean*sqr(V2[1]);
      err42[bin][0]=sqrt((2.+chi2sq*(4.+chi2sq))/(2.*rpmult[bin]))/chi2sq; 
      /* Eq.(D36a) */
      err42[bin][1]=sqrt((6.+chi2sq*(24.+chi2sq*(9.+chi2sq*(10.+4.*chi2sq))))/
			 (8.*rpmult[bin]))/sqr(chi2sq); /* Eq.(D36c) */

  } /* End loop over bins */

  vd22av[0]=vd22av[1]=vd42av[0]=vd42av[1]=0.;
  plotvdiff0=fopen("plot_vdiff-cum0.grf","w");
  plotvdiff1=fopen("plot_vdiff-cum1.grf","w");
  fprintf(output, "\n\nDIFFERENTIAL FLOW\n\n");
  for(int bin=0; bin<Nbins; bin++)
  {
    if (rpmult[bin]>0)
    { /* compute the average (over pT) v2 and v4 */
      vd22av[0]+=rpmult[bin]*v22[bin][0]; 
      vd22av[1]+=rpmult[bin]*v22[bin][1]; 
      vd42av[0]+=rpmult[bin]*v42[bin][0];
      vd42av[1]+=rpmult[bin]*v42[bin][1];
    }
    fprintf(output, " . bin %i: average detected multiplicity <M'>=%G\n", bin,
	    ((double) rpmult[bin])/((double) neve));
    fprintf(output, "   v'2/2{2}=%.3g +/- %.3g,", v22[bin][0], err22[bin][0]);
    fprintf(output, "   v'2/2{4}=%.3g +/- %.3g,\n", v22[bin][1], err22[bin][1]);
    fprintf(output, "   v'4/2{3}=%.3g +/- %.3g,", v42[bin][0], err42[bin][0]);
    fprintf(output, "   v'4/2{5}=%.3g +/- %.3g;\n\n", v42[bin][1], 
	    err42[bin][0]);
    fprintf(plotvdiff0, "%i  %.4g  %.4g  %.4g %.4g  %.4g\n", bin, 
            ((double) rpmult[bin])/((double) neve), v22[bin][0], err22[bin][0],
            v42[bin][0], err42[bin][0]);
    fprintf(plotvdiff1, "%i  %.4g  %.4g  %.4g %.4g  %.4g\n", bin, 
            ((double) rpmult[bin])/((double) neve), v22[bin][1], err22[bin][1],
            v42[bin][1], err42[bin][1]);
  }
  vd22av[0]/=(rmultmean*neve);
  vd42av[0]/=(rmultmean*neve);
  vd22av[1]/=(rmultmean*neve);
  vd42av[1]/=(rmultmean*neve);
  fprintf(output, "<v2{2}>=%.3g,  <v2{4}>=%.3g,  <v4{3}>=%.3g,  <v4{5}>=%.3g\n",
          vd22av[0],  vd22av[1], vd42av[0], vd42av[1]);

  gettimeofday (&Time, 0);
  fin=Time.tv_sec;
  fprintf(output, "\nEnd of the analysis, elapsed time=%d s\n", fin-debut);

  fclose(output);
  fclose(plotvdiff0);
  fclose(plotvdiff1);
}