/* Directed flow (v1) analysis through the generating-function zeros method. */

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

/* The program then reconstructs integrated elliptic flow (V2) 
   and integrated (V1) & differential (v1) directed flows. */

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

/* Throughout the paper, references are made to the paper:
   N. Borghini, J.-Y. Ollitrault, Nucl. Phys. A 742 (2004) 130 
   [nucl-th/0404087]. */

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

/* */

#define Pi (2.*acos(0.))
#define sqr(x) ((x)*(x))
#define rad2deg (180./Pi)
#define rootJ0 2.4048256  /* First zero of the modified Bessel function J0 */
#define J1rootJ0 0.519147 /* J1(rootJ0) */

/* PARAMETERS USED IN THE RECONSTRUCTION */
/* kr is the number of r values at which one computes the generating function 
   in the reconstruction of integrated elliptic flow V2; 
   kt1, kt2 are the numbers of theta1, theta2 values used in the analysis. */
#define kr 80
#define kt1 5
#define kt2 5

#define eps 0.05 /* the arbitrary parameter epsilon: 
                    see Eq.(6) and the discussion below. */

#define v2max 0.101 /* maximal possible integrated elliptic flow value */
#define v2step 0.001
#define v2min 0.021 /* minimal possible integrated flow value:
		       check that v2min + kr*v2step = v2max!! */

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

/* */

FILE *input, *data, *output, *plotVint, *plotvdiff; 

main()
{ /* main() */
  /* Variable declaration */

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

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

  double theta1[kt1], theta2[kt2]; /* theta values */
  double RZ[kr]; /* tabulation points for each theta1, theta2 values */

  /* Generating functions evaluated at the tabulation points: 
     - g[k2][k3] is the generating function g_0^{0,theta2[k2]} evaluated 
       at r=RZ[k3]; 
     - G[][] is its average over events; 
     - modGsq[k2][k3] is the squared modulus of G[k2][k3]. */
  complex<double> g[kt2][kr], G[kt2][kr]; 
  double modGsq[kt2][kr]; 
  double r0[kt2];  /* RZ values which minimize G[k2][] */
  double V2[kt2]; /* Integrated elliptic flow value for theta=theta2[k2]... */
  double V2mean; /* ... and averaged over the angle theta2. */

  /* - geps[k1][k2] is the generating function g_eps^{theta1[k1], theta2[k2]} 
       evaluated at r0[k2]; Geps[][] is its average over events; 
     - g0[k2] is the generating function g_0^{0, theta2[k2]} at r0[k2]; 
     - dlng_dz[k2] is the logarithmic derivative of g_0^{0, theta2[k2]} 
       with respect to z, evaluated at r0[k2]; 
     - dG_dz[] is its average over events of the derivative with respect to z 
       of g_0^{0, theta2[k2]} at r0[k2]; 
     - dlngeps_dw1[k1][k2][bin] is the logarithmic derivative 
     with respect to w1[bin] of g_eps^{theta1[k1], theta2[k2]}, 
     evaluated at r0[k2], for differential particles in bin "bin"; 
     - dGeps_dw1[][][] is the average over events of the derivative of geps
       with respect to w1 at r0[k2]. */
  complex<double> geps[kt1][kt2], Geps[kt1][kt2]; 
  complex<double> g0[kt2]; 
  complex<double> dlng_dz[kt2], dG_dz[kt2]; 
  complex<double> dlngeps_dw1[kt1][kt2][Nbins], dGeps_dw1[kt1][kt2][Nbins]; 

  double V1; /* Reconstructed integrated directed flow */
  double V1sq; /* "Squared" V1 (including the sign of V2) */
  double vd1[Nbins]; /* Differential directed flow v1 in each bin. */
  double v1av; /* Average v1, obtained by averaging the vd1[bin]. */

  double w1[Nbins], w2[Nbins]; /* weights */

  /* Event flow vectors, useful for the computation of statistical errors:
     - Q1/Q2 is the (weighted) flow vector in each event in harmonics n=1/2;
     - Q1mean/Q2mean is its average over events;
     - modQsqmean is the average over event of the squared modulus |Q|^2. */
  complex<double> Q1, Q1mean, Q2, Q2mean; 
  double modQ1sqmean, modQ2sqmean;

  /* - chi1, chi2 are the so-called "resolution parameters", following the 
       definition by J.-Y. Ollitrault in nucl-ex/9711003; 
     - err2[kt2], is the statistical error bar on the integrated flow estimate 
       V2[kt2]; err2mean is the error on V2mean;
     - err1m, err1p are the (asymmetric) error bars on V1;
     - err2mean is the average (over the angle theta2) of the err2[kt2];
     - err11[] is the statistical error bar on the differential flow estimate 
       vd1[] in the corresponding differential bin. */
  double chi2, err2[kt2], err2mean;
  double chi1, err1m, err1p, err11[Nbins]; 

  /* Maximal expected integrated elliptic flow V2
     & step between successive V2 values investigated */
  double V2max, V2step;

  int k; 
  double temp1, temp2, temp3[Nbins]; 

  struct timeval Time;
  int debut, fin;

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

  /* Initialization of angles: 
     for theta1, kt1 angles equally spaced between 0 and 2 Pi (radians) 
     for theta2, kt2 angles equally spaced between 0 and Pi/2 */
  for(int k1=0; k1<kt1; k1++) theta1[k1]=2.*k1*Pi/((float) kt1); 
  for(int k2=0; k2<kt2; k2++) theta2[k2]=((float) k2)*Pi/(2.*kt2); 

  /* The moduli of the tabulation points are chosen so as to correspond to 
     equally spaced (with step v2step) flow values. */
  data=fopen("inputdata.dat","r");
  fscanf(data, "%i %i\n", &k, &rmult); 
  fclose(data);
  V2max=v2max*rmult;
  V2step=v2step*rmult;
  for(int k3=0; k3<kr; k3++) RZ[k3]=rootJ0/(V2max-V2step*k3); 

  /* Initialization of the average generating function 
     used for the determination of V2[k2]. */
  for(int k2=0; k2<kt2; k2++) for(int k3=0; k3<kr; k3++) G[k2][k3]=0.;

  /* In this program, rmultmean is used to estimate statistical errors. */
  rmultmean=0.;
  for(int bin=0; bin<Nbins; bin++) 
  {
    rpmult[bin]=0.;
    w1[bin]=w2[bin]=1.; /* unit weights for the moment... */
  }

  /* Averages of the (complex!) event flow vectors and their squared moduli. */
  Q1mean=Q2mean=0.;
  modQ1sqmean=modQ2sqmean=0.;

  /* Scan of the data */
  printf("Reading the data file to analyse flow\n");
  input=fopen("events.dat", "r");
  fscanf(input, "%i\n", &neve); /* Read the number of events */
  printf("----------\nNumber of events in the run: %i\n", neve);

  output=fopen("zeros2.txt", "w");
  fprintf(output,
	  "FLOW ANALYSIS THROUGH THE ZEROS OF THE GENERATING FUNCTION\n\n");
  fprintf(output, "Parameters in the analysis:\n"); 
  fprintf(output, "%i values of theta1, %i values of theta2, epsilon=%.3f\n", 
	  kt1, kt2, eps); 
  fprintf(output, "\nAvailable statistics:\n  number of events = %i\n", 
	  neve);

  /* RECONSTRUCTION OF INTEGRATED FLOW FROM THE GENERATING FUNCTION */
  printf("\n 1. Reconstruction of integrated elliptic flow\n");
  /* */
  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);

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

    /* Initialization of the generating function g[k2][k3]. */
    for(int k2=0; k2<kt2; k2++) for(int k3=0; k3<kr; k3++) g[k2][k3]=1.;
    /* Initialization of the event flow vector in harmonic n=2. */
    Q2=0.;

    for(long np=0; np<rmult; np++) /* Loop over the event particles */
    {
      fscanf(input,"%f %f", &pT, &phi); /* Read the particle momentum & azimuth;
		                         the former is useless for the moment */
      binpT=pTbin(pT); /* Useless since unit weights... */

      /* For each value of theta2 and r, compute the generating functions g
         using Eq.(1) with epsilon=0 and theta1=0. */
      for(int k2=0; k2<kt2; k2++) /* Loop over theta2 values */
      {
        temp2=w2[binpT]*cos(2.*(phi-theta2[k2])); 
        for(int k3=0; k3<kr; k3++) g[k2][k3]*=(1.+i*RZ[k3]*temp2); 
      }/* End of the loop over theta values */

      /* Compute the event flow vector Q2. */
      Q2+=(w2[binpT]*exp(2.*i*((double) phi)));

    } /* End of the loop over particles */

    /* Compute the average over events of g[][] at each angle theta2 and for 
       each value of r */
    for(int k2=0; k2<kt2; k2++) for(int k3=0; k3<kr; k3++) G[k2][k3]+=g[k2][k3];
    
    Q2mean+=Q2;     /* Compute the average event flow vector... */ 
    modQ2sqmean+=sqr(abs(Q2)); /* and the average of its square modulus */

  } /* End of the loop over events */
  fclose(input);
  rmultmean/=neve; /* Average detected particle multiplicity */

  modQ2sqmean/=neve; /* Average square moduli of the event flow vector */

  gettimeofday(&Time, 0);
  fin=Time.tv_sec;
  printf("All events read, elapsed time=%d s\n", fin-debut);
  printf("Mean detected particle multiplicity <M>=%G\n", rmultmean);
  fprintf(output, "  mean detected particle multiplicity <M>=%G\n", 
	  rmultmean);
  fprintf(output, "  Minimum measurable elliptic flow v2=%.3g\n\n", 
	  rootJ0/sqrt(2.*rmultmean*log(0.5*neve))); 
  /* v = j_{01}/Sqrt(2M ln(Neve/2)) */

  /* RECONSTRUCTION OF INTEGRATED ELLIPTIC FLOW */
  V2mean=0.;

  for(int k2=0; k2<kt2; k2++) /* Loop over theta2 angles */
  { 
    for(int k3=0; k3<kr; k3++) modGsq[k2][k3]=sqr(abs(G[k2][k3])/neve); 
    /* Compute the square modulus of G[][] at each point. */
 
    /* Looking for the first minimum of |G[][k3]|^2 */
    k=0;
    while (modGsq[k2][k]>modGsq[k2][k+1]) k++;
    /* Interpolate the value of V2[k2], see footnote 3 in Nucl. Phys. A 727, 373
       and derive the corresponding position of the first zero of G */
    V2[k2]=V2max-V2step*(k+(modGsq[k2][k-1]-modGsq[k2][k+1])/2./
	   (modGsq[k2][k-1]-2.*modGsq[k2][k]+modGsq[k2][k+1]));
    r0[k2]=rootJ0/V2[k2]; /* cf. Eq.(4) */
    V2mean+=V2[k2]; /* Average V2 */
  } /* End of the loop over theta2 angles */

  V2mean/=kt2; /* Average V2 over angles */

  /* Compute the resolution parameter chi, in harmonic 2,
     using Eqs.(59),(62) of Nucl. Phys. A 727 (2003) 373. */
  chi2=V2mean/
    sqrt(modQ2sqmean-sqr(real(Q2mean)/neve)-sqr(imag(Q2mean)/neve)-sqr(V2mean));

  plotVint=fopen("plot_V2inf.grf", "w");
  temp2=0.;
  for(int k2=0; k2<kt2; k2++) 
  { /* Loop over theta2 angles */
    /* Compute the statistical error bar on the estimates V2^theta2{infty}, 
       using Eq.(90) of NPA 727, 373. */
    err2[k2]=V2[k2]*sqrt((exp(-sqr(rootJ0/chi2)/2.)*BesselJ0(2.*rootJ0)+
			  exp(sqr(rootJ0/chi2)/2.))/2./neve)/rootJ0/J1rootJ0;
    fprintf(plotVint, "%.3g\t%.4g\t%.4g\n", theta2[k2], V2[k2]/rmultmean, 
      err2[k2]/rmultmean); 
    /* Compute the statistical error bar on the average estimate V2{infty}, 
       using Eq.(89) of NPA 727, 373. */
    temp2+=exp(sqr(rootJ0/chi2)*cos(2.*theta2[k2])/2.)*
      BesselJ0(2.*rootJ0*sin(theta2[k2]))+
      exp(-sqr(rootJ0/chi2)*cos(2.*theta2[k2])/2.)*
      BesselJ0(2.*rootJ0*cos(theta2[k2]));    
  }
  fclose(plotVint);
  err2mean=V2mean*sqrt(temp2/2./neve/kt2)/rootJ0/J1rootJ0;

  /* */
  /* RECONSTRUCTION OF DIRECTED FLOW FROM THE GENERATING FUNCTION */
  /* */

  /* Initialization of the average generating functions */
  for(int k2=0; k2<kt2; k2++) 
  {
    dG_dz[k2]=0.; 
    for(int k1=0; k1<kt1; k1++) 
    {
      Geps[k1][k2]=0.;
      for(int bin=0; bin<Nbins; bin++) dGeps_dw1[k1][k2][bin]=0.;
    }
  }

  input=fopen("events.dat", "r");
  fscanf(input, "%i\n", &neve); /* Read the number of events */
  printf("\n 2. Reconstruction of integrated and differential directed flow\n");

  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);

    /* Initializations  */
    for(int k2=0; k2<kt2; k2++) 
    {
      g0[k2]=1.; 
      dlng_dz[k2]=0.; 
      for(int k1=0; k1<kt1; k1++) 
      {
        geps[k1][k2]=1.;
        for(int bin=0; bin<Nbins; bin++) dlngeps_dw1[k1][k2][bin]=0.;
      }
    }
    Q1=0.;

    fscanf(input, "%i", &rmult); /* Read the event multiplicity */
    for(long np=0; np<rmult; np++) /* Loop over the event particles */
    {
      /* Read the particle momentum and azimuth */
      fscanf(input, "%f %f", &pT, &phi); 
      /* Assign a bin to the particle according to its transverse momentum, 
         and increment the number of detected particles in that bin. */
      binpT=pTbin(pT);
      rpmult[binpT]+=1.;

      for(int k2=0; k2<kt2; k2++)
      { /* Compute the derivative dg_0/dz of the generating function g_0^theta2 
           for theta1=0, epsilon=0, at the minimum i*r0^theta2: 
           use the trick dg_0/dz(z) = g_0(z) * d(ln g_0)/dz, since d(ln g_0)/dz
           is easier to compute. 
           The derivative will be useful for the denominators of Eqs.(7) and (9)
           (once it is averaged over events). */
        temp2=w2[binpT]*cos(2.*(phi-theta2[k2]));
        g0[k2]*=(1.+i*r0[k2]*temp2);
        dlng_dz[k2]+=(temp2/(1.+i*r0[k2]*temp2)); 

        for(int k1=0; k1<kt1; k1++)
        { /* Compute the formal derivative dg_eps^{theta1,theta2}/dw_1 
             at the minimum r_0^theta2, for the numerator of Eq.(9). 
             One uses the same trick as above, using the derivative of the 
             logarithm of g_eps^{theta1,theta2}, see Eq.(8): 
             note that while the fraction depends on the bin (only the azimuths
             of particle in the differential bin under study should appear), 
             on the other hand the factor g_eps^{theta1,theta2}(i*r_0^theta2) is
             the same for all bins. */
          temp1=cos(phi-theta1[k1]); 
          geps[k1][k2]*=(1.+i*r0[k2]*(eps*w1[binpT]*temp1+temp2)); 
          dlngeps_dw1[k1][k2][binpT]+=
	    (temp1/(1.+i*r0[k2]*(eps*w1[binpT]*temp1+temp2))); 
        }
      }
      /* Compute the event flow vector in harmonic n=1. */
      Q1+=(w1[binpT]*exp(i*((double) phi)));

    } /* End of the loop over particles */

    for(int k2=0; k2<kt2; k2++) /* loop over theta2 values */
    { 
      /* Add the value of dg_0/dz(i*r_0^theta2) in the present event to the 
         corresponding average, for the denominators of Eqs.(7) and (9). */
      dG_dz[k2]+=dlng_dz[k2]*g0[k2]; 
      for(int k1=0; k1<kt1; k1++) /* loop over theta1 angles */
      {
        /* Compute the average over events of g_eps^{theta1,theta2} at the 
           minimum i*r_0^theta2, for the numerator of Eq.(7). */
        Geps[k1][k2]+=geps[k1][k2]; 
        /* In each differential bin, compute the average over differential 
           particles in that bin of dg_eps^{theta1,theta2}/dw1(i*r_0^theta2): 
           see Eqs.(8) and (9). */
        for(int bin=0; bin<Nbins; bin++)
          dGeps_dw1[k1][k2][bin]+=
	    (i*r0[k2]*eps*dlngeps_dw1[k1][k2][bin]*geps[k1][k2]);
      }
    }
    Q1mean+=Q1;     /* Compute the average event flow vector <Q1>... */ 
    modQ1sqmean+=sqr(abs(Q1)); /* and the average of its square modulus */
  } /* End of the loop over events */
  fclose(input);

  modQ1sqmean/=neve; /* Average square moduli of the event flow vector Q1 */

  /* Average of the generating functions... */
  for(int k2=0; k2<kt2; k2++)
  {
    dG_dz[k2]/=neve; /* ... over events; */
    for(int k1=0; k1<kt1; k1++) 
    {
      Geps[k1][k2]/=neve; /* ... over events: the double average that appears
        in the numerator of Eq.(7) can be rewritten as 
        <...>_{theta1;evts} = <cos(2(theta1-theta2)) * <...>_evts>_theta1, 
        so that one can average g_eps over events now, and perform the average
        over the theta1 angles later; */
      for(int bin=0; bin<Nbins; bin++) dGeps_dw1[k1][k2][bin]/=rpmult[bin]; 
      /* ... over differential particles in the bin!!! 
         Here again, the double average in the numerator of Eq.(9) is performed
         first over differential particles (<...>_psi), while the average over 
         theta1 will come later. */
    }
  }
  for(int bin=0; bin<Nbins; bin++) rpmult[bin]/=neve; 
  /* Compute the mean number of particles per bin */

  /* RECONSTRUCTION OF DIRECTED FLOW */
  V1sq=0.; /* We shall now compute an estimate of the product (V1)^2 sqn(V2), 
              through Eq.(7)... */
  for(int bin=0; bin<Nbins; bin++) vd1[bin]=0.; /* ... as well as an estimate 
              of differential flow v1[] in each differential bin, with the help
              of Eq.(9). */

  for(int k2=0; k2<kt2; k2++)
  {
    /* temp1 will be the sum over theta1 angles of the product of 
       cos(2(theta1-theta2)) times the the real part of the ratio
       <g_eps^{theta1,theta2}(..)>_evts / <dg_0^theta2/dz(..)>_evts, where both
       functions are evaluated at i*r_0^theta2. 
       Here I use the fact that the quantity 
       Re(<cos(2(theta1-theta2))*<...>_evts>_theta1 / <dg0^theta2/dz>_evts) in 
       Eq.(7) can be rewritten 
       <cos(2(theta1-theta2)) * Re(<...>_evts>/<dg0^theta2/dz>_evts)>_theta1. */
    temp1=0.; 
    /* Similarly, the real part appearing in Eq.(9) can be cast in the form
       <cos(2(theta1-theta2)) * Re(<...>_psi>/<dg0^theta2/dz>_evts)>_theta1, 
       and temp3[] will be, for each differential bin, the sum over theta1 of 
       the product of the cosine times the real part of the ratio of generating
       functions. */
    for(int bin=0; bin<Nbins; bin++) temp3[bin]=0.;
    for(int k1=0; k1<kt1; k1++)
    {
      temp1+=real(Geps[k1][k2]/dG_dz[k2])*cos(2.*(theta1[k1]-theta2[k2])); 

      for(int bin=0; bin<Nbins; bin++)
        temp3[bin]+=(real(dGeps_dw1[k1][k2][bin]/dG_dz[k2])*
		     cos(2.*(theta1[k1]-theta2[k2])));
    }
    /* Now we can divide temp1 by the number kt1 of theta1 values, thus 
       performing effectively the average over theta1.
       At once, we divide the (theta2-dependent) result by (r_0^theta2)^3, 
       which is the only other term depending on theta2 in the quantity that is 
       averaged over theta2 in Eq.(7)  */
    V1sq+=(temp1/pow(r0[k2],3.)/kt1); 
    /* Similarly, we now effectively perform the average over theta1 in Eq.(9),
       and compute the sum over theta2 of the theta2-dependent part. */
    for(int bin=0; bin<Nbins; bin++) vd1[bin]+=(temp3[bin]/pow(r0[k2],3.)/kt1);
  }
  
  /* We can now multiply V1sq by the theta2-independent prefactor 8j_{01}/eps^2 
     we had left aside till now, and divide by kt2 — thereby performing the 
     average over theta2. */
  V1sq*=(-8.*rootJ0/sqr(eps)/kt2); 
  /* Next, we compute the estimate of V1 (up to its sign), and correct the sign 
     of our previously-derived estimate of V2. */
  if(V1sq<0.) /* which actually means that V2 is negative */
  {
    V1=sqrt(-V1sq); 
    V2mean=-V2mean; 
  }
  else V1=sqrt(V1sq); 
  /* Eventually, we can compute the estimate of v1 in each bin, using the 
     estimtate V1 we just derived. 
     At the same time, we compute the average v1, average of the differential
     flow estimates over transverse momentum. */
  v1av=0.;
  for(int bin=0; bin<Nbins; bin++)
  { 
    vd1[bin]*=(-4*rootJ0/V1/sqr(eps)/kt2); 
    v1av+=vd1[bin]*rpmult[bin];
  }
  v1av/=rmultmean;


  /* V1 also allows us to compute the rsolution parameter chi1 */
  chi1=V1/
    sqrt(modQ1sqmean-sqr(real(Q1mean)/neve)-sqr(imag(Q1mean)/neve)-sqr(V1));

  temp2=0.;
  for(int k2=0; k2<kt2; k2++) 
    /* Loop over the theta2 angles, to compute the statistical error on the 
       estimates V1 and v1: temp2 will be the sum over the various theta2 
       values of the quantity which is averaged over theta2 in Eq.(28). 
       Note that the latter actually lacks an overall factor cos(2*theta2) 
       multiplying the quantity which is then averaged over theta2! 
       Here this factor is restored. */
    temp2+=(exp(sqr(rootJ0/chi2)*cos(2.*theta2[k2])/2.)*
      BesselJ0(2.*rootJ0*sin(theta2[k2])) -
      exp(-sqr(rootJ0/chi2)*cos(2.*theta2[k2])/2.)*
      BesselJ0(2.*rootJ0*cos(theta2[k2])))*cos(2.*theta2[k2]);
  /* Divide now this sum over theta2 by kt2 and by the factor 4*J1(j_01)^2, 
     see Eq.(28). */
  temp2*=(1/4./sqr(J1rootJ0)/kt2);

  /* Eq.(29) gives, together with Eq.(28), the statistical error bar on (V1)^2.
     This is turn yields asymmetric error bars on V1: the "true" value of V1
     is between sqrt(V1^2-d(V1^2)) and sqrt(V1^2+d(V1^2)). */
  err1m=
    V1-sqrt(sqr(V1)-2.*pow(V1,4.)*(1.+sqr(chi1))*temp2/kt2/pow(chi1,4.)/neve);
  err1p=
    sqrt(sqr(V1)+2.*pow(V1,4.)*(1.+sqr(chi1))*temp2/kt2/pow(chi1,4.)/neve)-V1;
  /* Statistical error bar on v1, given by Eqs.(28) and (30). 
     Remember that rpmult[] is now the average multiplicity per bin, so that 
     the total multiplicity in a bin (N' in the paper) is rpmult[]*neve. */
  for(int bin=0; bin<Nbins; bin++)
    err11[bin]=sqrt((1.+sqr(chi1))*temp2/rpmult[bin]/neve/kt2)/chi1;

  gettimeofday (&Time, 0);
  fin=Time.tv_sec;
  printf("All events read, elapsed time=%d s\n", fin-debut);
  printf("\nThe results will be available in the file zeros2.txt\n");

  /* Elliptic flow v2 */
  fprintf(output, "\nINTEGRATED ELLIPTIC FLOW V2\n");
  fprintf(output, "\taverage\n");
  fprintf(output, " V2/<M>:\t%.3g +/- %.3g\n", V2mean/rmultmean, 
          err2mean/rmultmean);

  /* Directed flow v1 */
  fprintf(output, "\nINTEGRATED DIRECTED FLOW V1\n");
  fprintf(output, " V1/<M>:\t%.3g +%.3g / -%.3g\n", V1/rmultmean, 
          err1p/rmultmean, err1m/rmultmean); 

  plotvdiff=fopen("plot_v1inf.grf", "w");
  fprintf(output, "\nDIFFERENTIAL DIRECTED FLOW\n");
  for(int bin=0; bin<Nbins; bin++)
  {
    fprintf(output, " . bin %i: average detected multiplicity <M'>=%G\n", bin,
	   rpmult[bin]);
    fprintf(output, "   v1=%.3g +/- %.3g\n", vd1[bin], err11[bin]); 
    fprintf(plotvdiff, "%i\t%g\t%.3g\t%.3g\n", bin, rpmult[bin], vd1[bin], 
            err11[bin]);
  }
  fclose(plotvdiff);
  fprintf(output, "\nAVERAGE DIRECTED FLOW <v1>=%.3g\n", v1av);

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

  fclose(output);
} /* main() */