/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2017, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*
* Component: MCPL_input
*
* %I
* Written by: Erik B Knudsen
* Date: Mar 2016
* Origin: DTU Physics
*
* Source-like component that reads neutron state parameters from an mcpl-file.
* %D
* Source-like component that reads neutron state parameters from a binary mcpl-file.
*
* MCPL is short for Monte Carlo Particle List, and is a new format for sharing events
* between e.g. MCNP(X), Geant4 and McStas.
*
* When used with MPI, the --ncount given on the commandline is overwritten by 
* #MPI nodes x #events in the file.
*
* %BUGS
*
* %P
* INPUT PARAMETERS
*
* filename: [str]       Name of neutron mcpl file to read.
* verbose: [ ]          Print debugging information for first 10 particles read.
* polarisationuse: [ ]  If !=0 read polarisation vectors from file.
* Emin: [meV]           Lower energy bound. Particles found in the MCPL-file below the limit are skipped.
* Emax: [meV]           Upper energy bound. Particles found in the MCPL-file above the limit are skipped.
* repeat_count: [1]     Repeat contents of the MCPL file this number of times. NB: When running MPI, repeating is implicit and is taken into account by integer division. Should be combined sith the _smear options!
* E_smear: [1]          When repeating events, make a Gaussian MC choice within E_smear*E around particle energy E
* pos_smear: [m]        When repeating events, make a flat MC choice of position within pos_smear around particle starting position
* dir_smear: [deg]      When repeating events, make a Gaussian MC choice of direction within dir_smear around particle direction
*
* %E
*******************************************************************************/

DEFINE COMPONENT MCPL_input
DEFINITION PARAMETERS ()
SETTING PARAMETERS (string filename=0, polarisationuse=1,verbose=1, Emin=0, Emax=FLT_MAX, int repeat_count=1, E_smear=0, pos_smear=0, dir_smear=0)
OUTPUT PARAMETERS (inputfile,nparticles,read_neutrons,used_neutrons,inactive)
DEPENDENCY "-Wl,-rpath,CMD(mcpl-config --show libdir) -ICMD(mcpl-config --show includedir) -LCMD(mcpl-config --show libdir) -lmcpl"

SHARE
%{
#include <mcpl.h>

%}

DECLARE
%{
mcpl_file_t inputfile;
long long nparticles;
long long read_neutrons;
long long used_neutrons;
int repeat_cnt;
int repeating;
int ismpislave;
int inactive;
%}

INITIALIZE
%{
char line[256];
long long ncount;

if (filename && strlen(filename) && strcmp(filename, "NULL") && strcmp(filename,"0") && repeat_count>0) {
  /* We got a proper filename, do the rest of the work */
  inactive=0;
  if(Emax<Emin){
        fprintf(stderr,"Warning(%s): Nonsensical energy interval: E=[%g,%g]. Aborting.\n",NAME_CURRENT_COMP,Emin,Emax);
        exit(-1);
    }
    /* No need to check if the file opens correctly since mcpl will
     * abort internally if it cannot open the file.*/
    inputfile = mcpl_open_file(filename);

  
    if ( !(nparticles=mcpl_hdr_nparticles(inputfile)) ) {
        fprintf(stderr,"Warning(%s): MCPL-file reports no present particles. Foolishly trying to go on.\n",NAME_CURRENT_COMP);
        nparticles=ncount;
    }else{
        printf("Message(%s): MCPL file (%s) produced with %s.\n",NAME_CURRENT_COMP,filename,mcpl_hdr_srcname(inputfile));
        printf("Message(%s): MCPL file (%s) contains %ld particles.\n",NAME_CURRENT_COMP,filename,nparticles); 
    }
    repeat_cnt = repeat_count;
    ismpislave=0;
#if defined (USE_MPI)
    repeat_cnt = ceil(1.0*repeat_cnt/mpi_node_count);
    ismpislave = mpi_node_rank;
    MPI_MASTER(
#endif
	     fprintf(stdout, "\n\n Warning: You are using MCPL_input with a repeat_count of %i:\n - Minimum neutron count requested is %ld x %i <= %ld",repeat_count,nparticles,repeat_count,repeat_cnt*nparticles); 
#if defined (USE_MPI)
  fprintf(stdout, " x %i MPI nodes = %ld neutrons total\n",mpi_node_count,mpi_node_count*repeat_cnt*nparticles);
     );
  mcset_ncount(mpi_node_count*repeat_cnt*nparticles);
#else
  fprintf(stdout, " neutrons total\n\n");
  mcset_ncount(repeat_cnt*nparticles);
#endif

    ncount=mcget_ncount();
    fprintf(stdout,"Initialize ncount is %ld\n",ncount);
    read_neutrons=0;
    used_neutrons=0;

#if defined (USE_MPI)   
    MPI_MASTER(
#endif

      if (verbose==1) {
        printf("MCPL_input verbose mode - outputting data on the 10 first read neutrons in MCPL units:\n");
      }

#if defined (USE_MPI) 
    );
    
#endif
    repeating = 0;
} else {
  /* Not presented with a reasonable filename, skip further action */
  inactive=1;
  printf("%s: filename is NULL or empty, skipping init and reads\n", NAME_CURRENT_COMP);
}
%}

TRACE
%{
  long long ncount;

  double nrm;
  const mcpl_particle_t *particle;// = (mcpl_particle_t *) calloc(sizeof(mcpl_particle_t),1);
  if(!inactive) {
    particle = mcpl_read(inputfile);

    ncount=mcget_ncount();
    // fprintf(stdout,"Trace ncount is %ld\n",ncount);      
    if (!particle) {
      if(repeat_cnt>1) {
	/* Trigger rewind of the file and ABSORB to get the first neutron "again" */
	repeating++;
	mcpl_rewind(inputfile);
	particle = mcpl_read(inputfile);
#if defined (USE_MPI)   
	MPI_MASTER(
#endif
       printf("MCPL inputfile %s rewound %i time(s)\n",filename,repeating);
#if defined (USE_MPI) 
    );
#endif
      } else
	ABSORB;       
    }
    if (particle->pdgcode!=2112) {
        /*Either no particle read, particle is not a neutron, or it has invalid energy - terminate to trigger next ray*/
       ABSORB;
    }
    read_neutrons++;
    
    /* check energy range*/
    if ( particle->ekin<Emin*1e-9 || particle->ekin>Emax*1e-9 ) {
        /*Particle energy out of range - terminate to trigger next ray*/
        ABSORB;
    }
    used_neutrons++;
    SCATTER;
#if defined (USE_MPI)   
  MPI_MASTER(
#endif

    if (verbose && used_neutrons<11) {
      printf("id=%"PRIu64" pdg=2112\tekin=%g MeV\tx=%g cm\ty=%g cm\tz=%g cm\tux=%g\tuy=%g\tuz=%g\tt=%g ms\tweight=%g\tpolx=%g\tpoly=%g\tpolz=%g\n",
	     read_neutrons, particle->ekin, particle->position[0], particle->position[1], particle->position[2],
	     particle->direction[0], particle->direction[1], particle->direction[2], particle->time, particle->weight,
	     particle->polarisation[0], particle->polarisation[1], particle->polarisation[2]);
     
    }

#if defined (USE_MPI)   
	    );
#endif    

    /*positions are in cm*/
    x=particle->position[0]/100;
    y=particle->position[1]/100;
    z=particle->position[2]/100;
    
    if (ismpislave || repeating) {
      double tmpx,tmpy,tmpz;
      // Position-MC:
      randvec_target_circle(&tmpx, &tmpy, &tmpz, NULL, 0, 0, 1, 0);
      NORM(tmpx,tmpy,tmpz);
      tmpx *= pos_smear*rand01(); tmpy *= pos_smear*rand01(); tmpz *= pos_smear*rand01();
      x+=tmpx; y+=tmpy; z+=tmpz;
    }
    
    if(polarisationuse){
        sx=particle->polarisation[0];
        sy=particle->polarisation[1];
        sz=particle->polarisation[2];
    }else{
        sx=sy=sz=0;
    }

    nrm = particle->ekin *1e9/VS2E;
    nrm = sqrt(nrm);
    if (ismpislave || repeating) {
      // Energy-MC:
      double tmp=(1.0+E_smear*randpm1());
      //printf("Adjusting energy from %g to",nrm);
      nrm *= (1+E_smear*randpm1());
      //printf(" to %g\n",nrm);
    }
    double d0=particle->direction[0],d1=particle->direction[1],d2=particle->direction[2];
    
    if (ismpislave || repeating) {
      // Direction-MC:
      double tmpx,tmpy,tmpz;
      // Position-MC, only in case of non-zero dir_smear):
      if (dir_smear) {
	randvec_target_circle(&d0, &d1, &d2, NULL, particle->direction[0], particle->direction[1], particle->direction[2], sin(dir_smear*DEG2RAD));
	NORM(d0,d1,d2);
      }
    }
    
    vx=d0*nrm;
    vy=d1*nrm;
    vz=d2*nrm;   

    /*time in ms:*/
    t=particle->time*1e-3;
    /*weight in unspecified units:*/
    p=particle->weight;

    /* Correct for repetition, by repeat_count and/or MPI */
    p /= repeat_cnt;
#if defined (USE_MPI)   
    p /= mpi_node_count;
#endif    
    SCATTER;
  }
%}

SAVE
%{
  if(!inactive) mcpl_close_file(inputfile);  
%}

FINALLY
%{
  if(!inactive) {
    long long ncount;

    ncount=mcget_ncount();
    
    if (used_neutrons!=read_neutrons){
        fprintf(stdout,"Message(%s): You have used %ld of %ld neutrons available in the MCPL file.\n",NAME_CURRENT_COMP,used_neutrons,read_neutrons);
    }
    if (ncount != used_neutrons){
        fprintf(stderr,"Warning (%s): You requested %ld neutrons from a file which contains %ld particles in general, of which only %ld are neutrons (within the wanted energy interval).\n"
                "Please examine the recorded intensities carefully.\n",NAME_CURRENT_COMP,ncount,nparticles,used_neutrons);
    }
  }
%}

MCDISPLAY
%{
    multiline(5, 0.2,0.2,0.0, -0.2,0.2,0.0, -0.2,-0.2,0.0, 0.2,-0.2,0.0, 0.2,0.2,0.0);
    /*M*/
    multiline(5,-0.085,-0.085,0.0, -0.085,0.085,0.0, -0.045,-0.085,0.0, -0.005,0.085,0.0, -0.005,-0.085,0.0);
    /*I*/
    line(0.045,-0.085,0, 0.045, 0.085,0);    
    line(0.005, 0.085,0, 0.085, 0.085,0);    
    line(0.005,-0.085,0, 0.085,-0.085,0);    
%}

END