#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <ctype.h>

#include "ldstab.h"

long int ntables;
scan_table_single_def *scan_table_all_poi;
corr_table_def corr_table;


//prototypes
double gauss_prob(double t);
double calculate_prob(long int cluster_size,double cluster_len,long int sample_size,double sample_len,scan_table_single_def *scan_table_all_poi,long int ntables,corr_table_def *corr_table,bool nocorr_flag);
long int find_parent(double *t_start_appo_poi,double *t_stop_appo_poi,long int iappo_last,double t_start_appo,double t_stop_appo);
long int find_peaks_and_howmany_sons(long int *parent_id_poi,long int *ipeak_v,long int *nclean_levels_v,long int *n_levels_above_root_v,bool *peak_flag_poi,long int *howmany_sons_poi,long int nele);
void eval_fwhm(double rho_bkg,double rho_ratio_ref_34,long int *level_id_poi,double *elen_poi,long int *nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int *howmany_sons_poi,long int nele,long int *ipeak_v,long int *nclean_levels_v,double *Tp_v,double *fwhm_v,long int npeaks,long int *interpol_word_v);
void eval_fwhm_v(double rho_bkg,double *rho_ratio_ref_v,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int *howmany_sons_poi,long int nele,long int *ipeak_v,long int *nclean_levels_v,double *Tp_v,double *fwhm_v,long int npeaks, long int *interpol_word_v);
void eval_fwhm_nomatter_of_nclean(double rho_bkg,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int nele,long int *ipeak_v,double *Tp_v,double *fwhm_v,long int npeaks, long int *interpol_word_v);
void eval_fwhm_iterate(double rho_bkg,double *rho_ratio_ref_v,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int *howmany_sons_poi,long int nele,long int *ipeak_v,long int *nclean_levels_v,double *Tp_v,double *fwhm_v,long int npeaks, long int *interpol_word_v);
double rho_average(double Delta_Tp,double Delta_T,double rho_bkg,double rho_p,double fwhm);;
double interpol(long int iinterpol_l,long int iinterpol_h,double *t_start_poi,double *t_stop_poi,long int *Nele_poi,double rho_ratio_l,double rho_ratio_h,double rho_ratio_ref);
FILE *init_survived_cluster_file(char *clusterfile,char *toutformat,char *xoutformat,long int N_tol,double prob_thr);
int fprintf_cluster(FILE *cfile,long int level_id,double t_start,double t_stop,long int Nele,double len,double elen,double chance_prob,char *toutformat,char *xoutformat);
void close_survived_cluster_file(FILE *cfile,long int nclusters);
FILE *init_peaklist_file(char *peaklist_file,char *toutformat,char *xoutformat,long int N_tol,double prob_thr,long int Npeaks);
void help(void);
char *stolower(char *str);
void *mymalloc(size_t size);

//*********************************************************************************
int main(int argc, char* argv[]){


  bool maxdev_flag=false;
  double nsigma=-1.;
  long int N_tol=-1;  
  char scantabledir[1000];
  bool scantabledir_flag=false;
  char clusterlistfile[1000];  
  char clusterlist_survived_file[1000];  
  char peaklist_file[1000];  

  FILE *cfile_survived;
  FILE *upeaklist;

   long int level_id_appo=0;
   double t_start_appo=0.;
   double t_stop_appo=0.;
   long int Nele_appo=0;
   double len_appo=0.;
   double elen_appo=0.;
   double rho_appo=0.;
   double chance_prob_appo=0.;

   double prob_thr=-1;
   double nsigmas=3.0;

   bool nocorr_flag=false;
   char toutformat[100];
   char xoutformat[100];
   char myarg[1000];
 
   bool peaklist_file_flag=false;


  int iarg;
  for (iarg=0;iarg<argc;iarg++){
    strcpy(myarg,argv[iarg]);
    if(strcmp(stolower(myarg),"help") == 0){
      help();
      return 0;
    }
  }



  strcpy(clusterlistfile,argv[1]);
  strcpy(clusterlist_survived_file,argv[2]);
  strcpy(scantabledir,"./");

  iarg=3;
  while (iarg<argc){
    if(strcmp(stolower(argv[iarg]),"maxdev") == 0){
      maxdev_flag=true;
    }else if(strcmp(stolower(argv[iarg]),"nocorr") == 0){
      nocorr_flag=true;
    }else if(strcmp(stolower(argv[iarg]),"scantabledir") == 0){
      if(iarg+1 < argc){
	scantabledir_flag=true;
	strcpy(scantabledir,argv[iarg+1]);
	printf("scantabledir: %s\n",scantabledir);
	iarg++;
      }
      
    }else if(strcmp(stolower(argv[iarg]),"peaklist_file") == 0){
      if(iarg+1 < argc){
	peaklist_file_flag=true;
	strcpy(peaklist_file,argv[iarg+1]);
	printf("writing to peaklist file: %s\n",peaklist_file);
	iarg++;
      }
      
    }else if(strcmp(stolower(argv[iarg]),"nsigma") == 0){
      if(iarg+1 < argc){
	nsigma=atof(argv[iarg+1]);
	if((nsigma > 3.5) || (nsigma < 2)){
	  fprintf(stderr,"NSIGMA value outside allowed range: [2,3.5]\n");
	  exit(-1);
	}
	iarg++;
      }
    }else if(strcmp(stolower(argv[iarg]),"thr_prob") == 0){
      if(iarg+1 < argc){
	prob_thr=atof(argv[iarg+1]);
	printf("thr prob: %lf\n",prob_thr);
	if((prob_thr < 2.E-4) || (prob_thr > 0.023)){
	  fprintf(stderr,"threshold probability outside allowed range: [2,3.5] std. dev.\n");
	  exit(-1);
	}
	iarg++;
      }
    }else{
      fprintf(stderr,"unknown parameter: %s, bye\n",argv[iarg]);
      exit(-1);
    }
    
    iarg++;
  }
  //printf("config read\n");

  scan_table_all_poi=NULL;
  printf("LOADING SCAN TABLES\n");
  if(scantabledir_flag){
    scan_table_all_poi=load_all_scan_table_fromdir(&ntables,scantabledir);
    printf("LOADING CORR TABLE\n");
     load_corr_table_fromdir(&corr_table,scantabledir);
  }else{
    scan_table_all_poi=load_all_scan_table(&ntables);
    printf("LOADING CORR TABLE\n");
    load_corr_table(&corr_table);
  }
  
  long int iele_ground=-1;
  
  
  if(nsigma > 0) printf("nsigma: %lf\n",nsigma);
  if(maxdev_flag) printf("working in maxdev mode\n");
  
  if(nsigma > 0) nsigmas=nsigma;
  if(prob_thr < 0)prob_thr=gauss_prob(nsigmas);
   
  char dummys[1000];
  printf("reading file with candidate clusters: %s\n",clusterlistfile);
  FILE *reportfile_stream=fopen(clusterlistfile,"r");
  if(reportfile_stream == NULL){
    fprintf(stderr,"unable to open input cluster file, bye\n");
    exit(-1);
  }

  long int nele_appo;      
  fscanf(reportfile_stream,"%s %li\n",dummys,&nele_appo);
  printf("nclusters: %li\n",nele_appo);
  fscanf(reportfile_stream,"%s %li\n",dummys,&N_tol);
  printf("tolerance parameter: %li\n",N_tol);
  fscanf(reportfile_stream,"%s %s\n",dummys,toutformat);
  printf("t format: %s\n",toutformat);
  fscanf(reportfile_stream,"%s %s\n",dummys,xoutformat);
  printf("X format: %s\n",xoutformat);
  fgets(dummys,1000,reportfile_stream);
  fgets(dummys,1000,reportfile_stream);
  fgets(dummys,1000,reportfile_stream);
  fgets(dummys,1000,reportfile_stream);
  fgets(dummys,1000,reportfile_stream);
  fgets(dummys,1000,reportfile_stream);
  fgets(dummys,1000,reportfile_stream);
   //exit(0);

  cfile_survived=init_survived_cluster_file(clusterlist_survived_file,toutformat,xoutformat,N_tol,prob_thr);

  long int *level_id_appo_poi=(long int *)mymalloc(sizeof(long int)*nele_appo);
  double *t_start_appo_poi=(double *)mymalloc(sizeof(double)*nele_appo);
  double *t_stop_appo_poi=(double *)mymalloc(sizeof(double)*nele_appo);
  long int *Nele_appo_poi=(long int *)mymalloc(sizeof(long int)*nele_appo);
  double *len_appo_poi=(double *)mymalloc(sizeof(double)*nele_appo);
  double *elen_appo_poi=(double *)mymalloc(sizeof(double)*nele_appo);
  double *chance_prob_appo_poi=(double *)mymalloc(sizeof(double)*nele_appo);
  long int *parent_id_appo_poi=(long int *)mymalloc(sizeof(long int)*nele_appo);

  //read the root and add it to the list of survived clusters:
  fscanf(reportfile_stream,"%li %lf %lf %li %lf %lf %lf\n",&level_id_appo,&t_start_appo,&t_stop_appo,&Nele_appo,&len_appo,&elen_appo,&rho_appo);
  //printf("%li %lf %lf %li %lf %lf %lf\n",level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,rho_appo);
  *chance_prob_appo_poi=0.;
  *level_id_appo_poi=level_id_appo;
  *t_start_appo_poi=t_start_appo;
  *t_stop_appo_poi=t_stop_appo;
  *Nele_appo_poi=Nele_appo;
  *len_appo_poi=len_appo;
  *elen_appo_poi=elen_appo;
  *parent_id_appo_poi=0;
  char format[1000];
  char format_root[1000];
  char format_cluster[1000];
  char format_maxdev[1000];
  strcpy(format,"%4li %8li ");
  strcat(format,toutformat);
  strcat(format," ");
  strcat(format,toutformat);
  strcat(format," %8li ");
  strcat(format,xoutformat);
  strcat(format," ");
  strcat(format,xoutformat);
  strcat(format," %10.3lE");
  strcat(format," %8.2lE");
   
  strcpy(format_root,"root:    ");
  strcat(format_root,format);
  strcat(format_root,"\n");
  
  strcpy(format_cluster,"survived:");
  strcat(format_cluster,format);
  strcat(format_cluster,"\n");
  
  strcpy(format_maxdev,"maxdev:  ");
  strcat(format_maxdev,format);
  strcat(format_maxdev,"\n");
  
  //printf(format_root,0,level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,rho_appo,chance_prob_appo);
  fprintf_cluster(cfile_survived,level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,chance_prob_appo,toutformat,xoutformat);  // print the root  

  double chance_prob_min=10.;
  long int iele_chance_prob_min=-1;
  long int iele_parent;
  //read all the other elements and add it to the list of survived clusters if they are relevant:
  long int iappo=1;  //the first element was the root
  while(fscanf(reportfile_stream,"%li %lf %lf %li %lf %lf %lf",&level_id_appo,&t_start_appo,&t_stop_appo,&Nele_appo,&len_appo,&elen_appo,&rho_appo) == 7){     
    //printf("%li %lf %lf %li %lf %lf %lf\n",level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,rho_appo);
    if(maxdev_flag){
      iele_parent=0;  //For the maxdev case the comparison is performed with root alsways 
    }else{
      iele_parent=find_parent(t_start_appo_poi,t_stop_appo_poi,iappo,t_start_appo,t_stop_appo);   //find parent
    }
    long int Nele_parent=*(Nele_appo_poi+iele_parent);
     double elen_parent=*(elen_appo_poi+iele_parent); 
     chance_prob_appo=calculate_prob(Nele_appo,len_appo,Nele_parent,elen_parent,scan_table_all_poi,ntables,&corr_table,nocorr_flag);   //eval if cluster is relevant
     if(chance_prob_appo <= prob_thr){
       *(chance_prob_appo_poi+iappo)=chance_prob_appo;
       *(level_id_appo_poi+iappo)=level_id_appo;
       *(t_start_appo_poi+iappo)=t_start_appo;
       *(t_stop_appo_poi+iappo)=t_stop_appo;
       *(Nele_appo_poi+iappo)=Nele_appo;
       *(len_appo_poi+iappo)=len_appo;
       *(elen_appo_poi+iappo)=elen_appo;
       *(parent_id_appo_poi+iappo)=iele_parent;
       //printf(format_cluster,iappo,level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,rho_appo,chance_prob_appo);
       if(!maxdev_flag){
	 fprintf_cluster(cfile_survived,level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,chance_prob_appo,toutformat,xoutformat);  // print the survived cluster
       }
       //printf("parent id:%li\n",iele_parent);
       if(chance_prob_appo < chance_prob_min){
	 chance_prob_min=chance_prob_appo;
	 iele_chance_prob_min=iappo;
       }
       iappo++;
     }
  }
  
  if(!maxdev_flag){
    close_survived_cluster_file(cfile_survived,iappo);
  }else{
    if(iele_chance_prob_min > 0){
      chance_prob_appo=*(chance_prob_appo_poi+iele_chance_prob_min);
      level_id_appo=*(level_id_appo_poi+iele_chance_prob_min);
      t_start_appo=*(t_start_appo_poi+iele_chance_prob_min);
      t_stop_appo=*(t_stop_appo_poi+iele_chance_prob_min);
      Nele_appo=*(Nele_appo_poi+iele_chance_prob_min);
      len_appo=*(len_appo_poi+iele_chance_prob_min);
      elen_appo=*(elen_appo_poi+iele_chance_prob_min);
      rho_appo=(double)Nele_appo/elen_appo;
      //printf(format_maxdev,iappo,level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,rho_appo,chance_prob_appo);
      fprintf_cluster(cfile_survived,level_id_appo,t_start_appo,t_stop_appo,Nele_appo,len_appo,elen_appo,chance_prob_appo,toutformat,xoutformat);  // print the survived cluster
    }else{
      printf("no significant cluster found\n");
    }
    close_survived_cluster_file(cfile_survived,2);
  }
  

  if(maxdev_flag){
    goto end_proc;
  }

  //qui c'è la scrittura dei dati
  long int nele=iappo;
  long int npeaks=nele;   //the number of peaks is lower than nele; but I still do not know
  bool *peak_flag_appo_poi=(bool *)mymalloc(sizeof(bool)*nele);   //it is true if the elemwent is a peak
  long int *howmany_sons_appo_poi=(long int *)mymalloc(sizeof(long int)*nele);                             
  long int *ipeak_v=(long int *)mymalloc(sizeof(long int)*nele);  
  long int *nclean_levels_v=(long int *)mymalloc(sizeof(long int)*nele);
  long int *n_levels_above_root_v=(long int *)mymalloc(sizeof(long int)*nele);
  npeaks=find_peaks_and_howmany_sons(parent_id_appo_poi,ipeak_v,nclean_levels_v,n_levels_above_root_v,peak_flag_appo_poi,howmany_sons_appo_poi,nele);
  double rho_bkg=0.;

  if(peaklist_file_flag){
    upeaklist=init_peaklist_file(peaklist_file,toutformat,xoutformat,N_tol,prob_thr,npeaks);
  }

  if(npeaks >0){
    double *Tp_v=(double *)mymalloc(sizeof(double)*npeaks);
    double *fwhm_v=(double *)mymalloc(sizeof(double)*npeaks);
    long int *interpol_word_v=(long int *)mymalloc(sizeof(long int)*npeaks);
    double rho_ratio_ref_34=0.75;
    eval_fwhm(rho_bkg,rho_ratio_ref_34,level_id_appo_poi,elen_appo_poi,Nele_appo_poi,t_start_appo_poi,t_stop_appo_poi,parent_id_appo_poi,howmany_sons_appo_poi,nele,ipeak_v,nclean_levels_v,Tp_v,fwhm_v,npeaks,interpol_word_v);
    double *Tp78_v=(double *)mymalloc(sizeof(double)*npeaks);
    double *fwhm78_v=(double *)mymalloc(sizeof(double)*npeaks);
    long int *interpol_word78_v=(long int *)mymalloc(sizeof(long int)*npeaks);
    rho_ratio_ref_34=7./8;
    eval_fwhm(rho_bkg,rho_ratio_ref_34,level_id_appo_poi,elen_appo_poi,Nele_appo_poi,t_start_appo_poi,t_stop_appo_poi,parent_id_appo_poi,howmany_sons_appo_poi,nele,ipeak_v,nclean_levels_v,Tp78_v,fwhm78_v,npeaks,interpol_word78_v);

    
    double *Tp_i_v=(double *)mymalloc(sizeof(double)*npeaks);
    double *fwhm_i_v=(double *)mymalloc(sizeof(double)*npeaks);
    long int *interpol_word_i_v=(long int *)mymalloc(sizeof(long int)*npeaks);
    double *rho_ratio_ref_i_v=(double *)mymalloc(sizeof(long int)*npeaks);

    eval_fwhm_iterate(rho_bkg,rho_ratio_ref_i_v,level_id_appo_poi,elen_appo_poi,Nele_appo_poi,t_start_appo_poi,t_stop_appo_poi,parent_id_appo_poi,howmany_sons_appo_poi,nele,ipeak_v,nclean_levels_v,Tp_i_v,fwhm_i_v,npeaks,interpol_word_i_v);


    double *Tp_nomatter_v=(double *)mymalloc(sizeof(double)*npeaks);
    double *fwhm_nomatter_v=(double *)mymalloc(sizeof(double)*npeaks);
    long int *interpol_word_nomatter_v=(long int *)mymalloc(sizeof(long int)*npeaks);
    eval_fwhm_nomatter_of_nclean(rho_bkg,level_id_appo_poi,elen_appo_poi,Nele_appo_poi,t_start_appo_poi,t_stop_appo_poi,parent_id_appo_poi,nele,ipeak_v,Tp_nomatter_v,fwhm_nomatter_v,npeaks,interpol_word_nomatter_v);

    //print peak infos
    long int ipeak;


    char format_peak[1000];
    strcpy(format_peak,"%4li ");
    strcat(format_peak,toutformat);
    strcat(format_peak," %6li ");
    strcat(format_peak,toutformat);
    strcat(format_peak," %10.3lE ");
    strcat(format_peak,xoutformat);
    strcat(format_peak," ");
    strcat(format_peak,toutformat);
    strcat(format_peak," %+6li %6li\n");


    for(ipeak=0;ipeak<npeaks;ipeak++){
      double Tp,Tp_nomatter;
      double fwhm,fwhm78,fwhm_nomatter,fwhm_i;
      long int interpol_word,interpol_word78,interpol_word_i,interpol_word_nomatter;
      long int n_levels_above_root;
      Tp=*(Tp_v+ipeak);
      Tp_nomatter=*(Tp_nomatter_v+ipeak);
      fwhm=*(fwhm_v+ipeak);
      fwhm78=*(fwhm78_v+ipeak);
      fwhm_i=*(fwhm_i_v+ipeak);
      fwhm_nomatter=*(fwhm_nomatter_v+ipeak);
      interpol_word=*(interpol_word_v+ipeak);
      interpol_word78=*(interpol_word78_v+ipeak);
      interpol_word_i=*(interpol_word_i_v+ipeak);
      interpol_word_nomatter=*(interpol_word_nomatter_v+ipeak);
      
      long int index_peak=*(ipeak_v+ipeak);
      long int csize_leaf=*(Nele_appo_poi+index_peak);
      double elen_leaf=*(elen_appo_poi+index_peak);
      double rho_leaf=(double)csize_leaf/elen_leaf;
      double dTp=(*(t_stop_appo_poi+index_peak)-*(t_start_appo_poi+index_peak))*(1.+1./csize_leaf);
      n_levels_above_root=*(n_levels_above_root_v+ipeak);
      //printf("ipeak %4li, Tp %13.7lf, Np %6li, Dtp %8.3lf, rho_p %10.3lE, fwhm %8.3lf %8.3lf, iw %+6li %+6li\n",ipeak,Tp,csize_leaf,dTp,rho_leaf,fwhm,fwhm_nomatter,interpol_word,interpol_word_nomatter);
      //printf("ipeak %4li, Tp %13.7lf, Np %6li, Dtp %8.3lf, rho_p %10.3lE, fwhm %8.3lf %8.3lf %8.3lf %8.3lf, iw %+6li %+6li %+6li %+6li\n",ipeak,Tp,csize_leaf,dTp,rho_leaf,fwhm,fwhm78,fwhm_i,fwhm_nomatter,interpol_word,interpol_word78,interpol_word_i,interpol_word_nomatter);

      if(peaklist_file_flag){
	fprintf(upeaklist,format_peak,index_peak,Tp,csize_leaf,dTp,rho_leaf,elen_leaf,fwhm_i,interpol_word_i,n_levels_above_root);
      }else{
	printf(format_peak,index_peak,Tp,csize_leaf,dTp,rho_leaf,elen_leaf,fwhm_i,interpol_word_i,n_levels_above_root);
      }
    }

  }

  if(peaklist_file_flag){
    fclose(upeaklist);
  }
  

  
 end_proc:
  
  //printf("closing all\n");
  fclose(reportfile_stream);
  free(scan_table_all_poi);
  return 0;
  }


//*****************************************
double calculate_prob(long int cluster_size,double cluster_len,long int sample_size,double sample_len,scan_table_single_def *scan_table_all_poi,long int ntables,corr_table_def *corr_table,bool nocorr_flag){


  double mdist_observed=cluster_len/sample_len;

  double *mdist_sigma_poi=(double *)mymalloc(10*sizeof(double));

  interpol_mdist_scan_stat(sample_size,cluster_size,scan_table_all_poi,mdist_sigma_poi, &ntables);    
  int result=interpol_mdist_scan_stat(sample_size,cluster_size,scan_table_all_poi,mdist_sigma_poi, &ntables);
      
      
  double prob_1sigma=1.587E-1;
  double prob_2sigma=2.275E-2;
  double prob_3sigma=1.345E-3;
  double prob_4sigma=3.165E-05;
  double prob_5sigma=2.857E-07;
  
  
//  double prob_0p5sigma=gauss_prob(0.5);
//  double prob_1p5sigma=gauss_prob(1.5);
//  double prob_2p5sigma=gauss_prob(2.5);
//  double prob_3p5sigma=gauss_prob(3.5);
//  double prob_4p5sigma=gauss_prob(4.5);
  
   
  double prob_this;
  double nsigma_this=-1;
  if(mdist_observed <= *(mdist_sigma_poi+4)){
    nsigma_this=5.;
  }else if(mdist_observed <= *(mdist_sigma_poi+9)){
    nsigma_this=4.5+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+9))/(*(mdist_sigma_poi+4)-*(mdist_sigma_poi+9));
  }else if(mdist_observed <= *(mdist_sigma_poi+3)){
    nsigma_this=4.0+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+3))/(*(mdist_sigma_poi+9)-*(mdist_sigma_poi+3));
  }else if(mdist_observed <= *(mdist_sigma_poi+8)){
    nsigma_this=3.5+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+8))/(*(mdist_sigma_poi+3)-*(mdist_sigma_poi+8));
  }else if(mdist_observed <= *(mdist_sigma_poi+2)){
    nsigma_this=3.0+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+2))/(*(mdist_sigma_poi+8)-*(mdist_sigma_poi+2));
  }else if(mdist_observed <= *(mdist_sigma_poi+7)){
    nsigma_this=2.5+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+7))/(*(mdist_sigma_poi+2)-*(mdist_sigma_poi+7));
  }else if(mdist_observed <= *(mdist_sigma_poi+1)){
    nsigma_this=2.0+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+1))/(*(mdist_sigma_poi+7)-*(mdist_sigma_poi+1));
  }else if(mdist_observed <= *(mdist_sigma_poi+6)){
    nsigma_this=1.5+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+6))/(*(mdist_sigma_poi+1)-*(mdist_sigma_poi+6));
  }else if(mdist_observed <= *(mdist_sigma_poi+0)){
    nsigma_this=1.0+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+0))/(*(mdist_sigma_poi+6)-*(mdist_sigma_poi+0));
  }else if(mdist_observed <= *(mdist_sigma_poi+5)){
    nsigma_this=0.5+0.5*((double)(mdist_observed)-*(mdist_sigma_poi+5))/(*(mdist_sigma_poi+0)-*(mdist_sigma_poi+5));
  }

  if(nsigma_this < 0.){
    return 1.;
  }

  if(nocorr_flag){
    prob_this=gauss_prob(nsigma_this);
    //printf("nsigma: %lf\n",nsigma_this);
    return prob_this;
  }
  
  return interpol_correction_table(sample_size,nsigma_this,corr_table);  
}

//*****************************************************************************++
long int find_parent(double *t_start_appo_poi,double *t_stop_appo_poi,long int iappo_last,double t_start_appo,double t_stop_appo){
  //find the parent within the actual list of survived clusters

  //t_start_appo and t_stop_appo are the boundaries of the cluster for which we have to find the parent 
  //iappo_last is the index of the last cluster among the actual list of survived clusters
  
  //bool parent_found=false;
  long int iele=iappo_last;
  //remember that the list of clusters is ordered with respect to threshold
  while(iele >= 0){
    double t_start_parent_candidate=*(t_start_appo_poi+iele);
    double t_stop_parent_candidate=*(t_stop_appo_poi+iele);

    if((t_start_parent_candidate <= t_start_appo) && (t_stop_parent_candidate >= t_start_appo)){
      //parent_found=true;
      return iele;
    }
    iele--;
  }
  fprintf(stderr,"parent not found, why?\n");
  exit(-1);
  
  return 0;
}

//********************************************
double gauss_prob(double t){
  
  double result=(double)1.-((erf(t/pow(2.,0.5)))/2+0.5);
  return result;
}
 
//*************************************
char *stolower(char *str){
  long int i;
  for(i = 0; str[i]; i++){
    str[i] = tolower(str[i]);
  }
  return str;
}
 
//****************************************
 long int find_peaks_and_howmany_sons(long int *parent_id_poi,long int *ipeak_v,long int *nclean_levels_v,long int *n_levels_above_root_v,bool *peak_flag_poi,long int *howmany_sons_poi,long int nele){
  
  //finds the peaks of the light curve (searches if a generic cluster has a parent. If not, it is a peak)
  
  long int ipeak=0;
  long int npeaks;
  long int iele;
  long int iele2;
  
  long int parent_id;
  
  for(iele=0;iele<nele;iele++){
    *(peak_flag_poi+iele)=true;
    *(howmany_sons_poi+iele)=0;
    for(iele2=iele+1;iele2<nele;iele2++){
       parent_id=*(parent_id_poi+iele2);
       if(parent_id == iele){
	 *(peak_flag_poi+iele)=false;
	 *(howmany_sons_poi+iele)= *(howmany_sons_poi+iele)+1;
       }
    }
    if(*(peak_flag_poi+iele)){
      *(ipeak_v+ipeak)=iele;
      //printf("cluster: %li is a peak\n",iele);
      ipeak++;
    }
  }
  npeaks=ipeak;
  
  long int ipeak_index;
  long int nbrothers;
  for(ipeak=0;ipeak<npeaks;ipeak++){
    *(nclean_levels_v+ipeak)=1;
    ipeak_index=*(ipeak_v+ipeak);
    parent_id=*(parent_id_poi+ipeak_index);
    nbrothers=*(howmany_sons_poi+parent_id);
    while((nbrothers == 1) && (parent_id != 0)){    //while it is the only son, and the parent is not the root 
      *(nclean_levels_v+ipeak)=*(nclean_levels_v+ipeak)+1;
      parent_id=*(parent_id_poi+parent_id);
      nbrothers=*(howmany_sons_poi+parent_id);
    }
    //printf("flare with peak cluster %li has %li clean levels\n",ipeak_index,*(nclean_levels_v+ipeak));
  }

  for(ipeak=0;ipeak<npeaks;ipeak++){
    ipeak_index=*(ipeak_v+ipeak);
    parent_id=*(parent_id_poi+ipeak_index);
    *(n_levels_above_root_v+ipeak)=1;
    while(parent_id != 0){    //while it is the only son, and the parent is not the root 
      *(n_levels_above_root_v+ipeak)=*(n_levels_above_root_v+ipeak)+1;
      parent_id=*(parent_id_poi+parent_id);
    }
  }


  
  return npeaks;
}

//****************************************************
void eval_fwhm(double rho_bkg,double rho_ratio_ref_34,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int *howmany_sons_poi,long int nele,long int *ipeak_v,long int *nclean_levels_v,double *Tp_v,double *fwhm_v,long int npeaks, long int *interpol_word_v){

  //fwhm is evaluated assuming flares with trapezoidal shape: the flare peak defines the higher base of the trapezoid,
  //the length of the higher base is = to elen_peak (the effective length of the cluster defining the peak).
  //for extremely large FWHM, the mean density within the FWHM is 3/4 the event density at the peak,
  //for FWHM = 2*elen_peak  , the mean density within the FWHM is 7/8 the event density at the peak,
  //for an hat shaped flare FWHM= ele_peak , the mean density within the FWHM is = the event density at the peak.
  //two methods for evaluate the FWHM:
  //1)assume the FWHM always large and eval the effective length of an hypothetical cluster with
  //  event density == 3/4 event density of the peak cluster.
  //2)interpolate FWHM at the value of the density which equal the foreseen density for a flare with FWHM=FWHM_try
  //If there is only one cluster above the root.

  long int ipeak;
  long int ipeak_index;
  long int iparent_index;
  long int iactual_level_index;
  long int nclean_levels;
  long int iclean_level;

  long int parent_level_id;
  long int parent_multiplicity;

  bool good_quality_on_levels_flag;


  long int interpol_word34;
  double fwhm_time;
  double fwhm_time34;
  double Tp_time;
  double Delta_Tp_time;
  double Delta_time_l;
  double Delta_Tp;
  double Delta_T;
  double fwhm_try;
  double rho_p,rho_estimate,rho_level,rho_ratio,rho_ratio_l,rho_ratio_l2,rho_ratio_h,rho_ratio_h2;
  double rho_ratio_withp,rho_ratio_withp_l,rho_ratio_withp_l2,rho_ratio_withp_h,rho_ratio_withp_h2,rho_ratio_old;
  long int iinterpol_l,iinterpol_l2,iinterpol_h,iinterpol_h2;
  long int iinterpol34_l,iinterpol34_l2,iinterpol34_h,iinterpol34_h2;
  double rho_ratio_ref=1.;
  double rho_ratio_ref_34_old=0.75;
  //double rho_ratio_ref_34=7./8;   //when fwhm=2*higher_base
  for(ipeak=0;ipeak<npeaks;ipeak++){
    interpol_word34=-100;
    good_quality_on_levels_flag=false;
    nclean_levels=*(nclean_levels_v+ipeak);
    fwhm_time=-1.E9;
    ipeak_index=*(ipeak_v+ipeak);
    Delta_Tp=*(elen_poi+ipeak_index);
    Tp_time=(*(t_start_poi+ipeak_index)+*(t_stop_poi+ipeak_index))/2;
    Delta_Tp_time=(*(t_stop_poi+ipeak_index)-*(t_start_poi+ipeak_index))*(1.+1./(*(Nele_poi+ipeak_index)));
    //printf("\n\nprocessing peak %li, index %li, nclean levels: %li, Tp %lf, Delta_Tp %lf\n",ipeak,ipeak_index,nclean_levels,Tp_time,Delta_Tp_time);
    rho_p=(double)(*(Nele_poi+ipeak_index))/(*(elen_poi+ipeak_index))-rho_bkg;
    //printf("rho_p %lE, elen %lE, size %li\n",rho_p,*(elen_poi+ipeak_index),*(Nele_poi+ipeak_index));
    if(nclean_levels >= 2){   //nclean_levels counts the peak also
      good_quality_on_levels_flag=true;
      //rho_p=(*elen_poi+ipeak_index)/(*(Nele_poi+ipeak_index))-rho_bkg;
      iactual_level_index=*(parent_id_poi+ipeak_index);
      rho_ratio_old=-1.;
      iinterpol_l=-1;
      iinterpol_l2=-1;
      iinterpol_h=-1;
      iinterpol_h2=-1;
      iinterpol34_l=-1;
      iinterpol34_l2=-1;
      iinterpol34_h=ipeak_index;
      rho_ratio_withp_h=1.;
      iinterpol34_h2=ipeak_index;
      rho_ratio_withp_h2=1.;
      for(iclean_level=0;iclean_level<nclean_levels-1;iclean_level++){   //eval on no more than nclean_levels-1 (nclean_levels counts the peak also)
	//rho_level=(*elen_poi+iactual_level_index)/(*(Nele_poi+iactual_level_index));
	rho_level=(double)(*(Nele_poi+iactual_level_index))/(*(elen_poi+iactual_level_index));
	Delta_T=*(elen_poi+iactual_level_index);
	fwhm_try=Delta_T;
	rho_estimate=rho_average(Delta_Tp,Delta_T,rho_bkg,rho_p,fwhm_try);
	rho_ratio=rho_estimate/rho_level;
	rho_ratio_withp=rho_level/rho_p;
	//printf("      rho_p %lE, rho_bkg %lE, rho_level %lE, rho est %lE\n",rho_p,rho_bkg,rho_level,rho_estimate);
	//printf("      elen p %lE, elen level %lE\n",Delta_Tp,Delta_T);
	//printf("      iclean level %li, level %li, rho_ratio: %lf, rho_ratio2 %lf\n",iclean_level,iactual_level_index,rho_ratio,rho_ratio_withp);

	//case interpol
	if(rho_ratio <1){
	  iinterpol_h2=iinterpol_h;
	  iinterpol_h=iactual_level_index;
	  rho_ratio_h2=rho_ratio_h;
	  rho_ratio_h=rho_ratio;
	}
	if(rho_ratio >= 1){
	  if(rho_ratio_old < 1){
	    iinterpol_l=iactual_level_index;
	    rho_ratio_l=rho_ratio;
	  //printf("    preparing for interpolation %li %li\n",iinterpol_l,iinterpol_h);
	  }else{
	    if(iinterpol_l2 < 0){
	      iinterpol_l2=iactual_level_index;
	      rho_ratio_l2=rho_ratio;
	    }
	  }
	}

	//case interpol3/4
	if((rho_ratio_withp >= rho_ratio_ref_34)){
	  iinterpol34_h2=iinterpol34_h;
	  iinterpol34_h=iactual_level_index;
	  rho_ratio_withp_h=rho_ratio_withp;
	}
	if(rho_ratio_withp <= rho_ratio_ref_34){
	  if(iinterpol34_l < 0){
	    iinterpol34_l=iactual_level_index;
	    rho_ratio_withp_l=rho_ratio_withp;
	    //printf("    preparing for interpolation 3/4 %li %li\n",iinterpol34_l,iinterpol34_h);
	  }else{
	    if(iinterpol34_l2 < 0){
	      iinterpol34_l2=iactual_level_index;
	      rho_ratio_withp_l2=rho_ratio_withp;
	    }
	  }
	}

	iactual_level_index=*(parent_id_poi+iactual_level_index);
	rho_ratio_old=rho_ratio;
      }  //end for cicle
      if((iinterpol_h >= 0) && (iinterpol_l >= 0)){
	//printf("    preparing for interpolation %li %li\n",iinterpol_l,iinterpol_h);
	fwhm_time=interpol(iinterpol_l,iinterpol_h,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_l,rho_ratio_h,rho_ratio_ref);
	//printf("    >>>fwhm_time: %lf\n",fwhm_time);
	//printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time,rho_ratio_l,rho_ratio_h);
      }else{
	if((iinterpol_l2 >= 0) && (iinterpol_l >= 0)){
	  //printf("    preparing for extrapolating back %li %li %li\n",iinterpol_l,iinterpol_l2,iinterpol_h);
	  fwhm_time=interpol(iinterpol_l,iinterpol_l2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_l,rho_ratio_l2,rho_ratio_ref);
	  Delta_time_l=*(t_stop_poi+iinterpol_l)-*(t_start_poi+iinterpol_l);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time);
	  //printf("    >>>fwhm_time: %lf %lf %lf %lf\n",fwhm_time,rho_ratio_l,rho_ratio_l2,Delta_time_l);
	}else if((iinterpol_h >= 0) && (iinterpol_h2 >= 0)){
	  //printf("    preparing for extrapolating forward %li %li %li\n",iinterpol_l,iinterpol_l2,iinterpol_h,iinterpol_h2);
	  fwhm_time=interpol(iinterpol_h,iinterpol_h2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_h,rho_ratio_h2,rho_ratio_ref);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time,rho_ratio_h,rho_ratio_h2);	  
	}else{
	  Delta_time_l=*(t_stop_poi+iinterpol_l)-*(t_start_poi+iinterpol_l);
	  //printf("    not enough points (iinterpol_l2 < 0) , rho_ratio_at_l %lf %lf\n",rho_ratio_l,Delta_time_l);
	}
      }
      if((iinterpol34_h >= 0) && (iinterpol34_l >= 0)){
	interpol_word34=nclean_levels;
	//printf("    preparing for interpolation 3/4 %li %li %li\n",iinterpol34_l,iinterpol34_h,interpol_word34);
	fwhm_time34=interpol(iinterpol34_l,iinterpol34_h,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_l,rho_ratio_withp_h,rho_ratio_ref_34);
	//printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	//printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_l,rho_ratio_withp_h);
      }else{
	interpol_word34=-nclean_levels;
	//printf("    preparing for extrapolation 3/4 forward %li %li %li %li %li\n",iinterpol34_l,iinterpol34_l2,iinterpol34_h,iinterpol34_h2,interpol_word34);
	  fwhm_time34=interpol(iinterpol34_h,iinterpol34_h2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_h,rho_ratio_withp_h2,rho_ratio_ref_34);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_h,rho_ratio_withp_h2);
      }


    }else{
      iparent_index=*(parent_id_poi+ipeak_index);
      parent_level_id=*(level_id_poi+iparent_index);
      parent_multiplicity=*(howmany_sons_poi+iparent_index);
      //printf("  no clean levels, proceed with the parent\n");
      //printf("    parent level id: %li, parent multiplicity %li\n",parent_level_id,parent_multiplicity);
      if(parent_level_id != -1){
	rho_level=(double)(*(Nele_poi+iparent_index))/(*(elen_poi+iparent_index));
	//Delta_T=*(elen_poi+iparent_index)/2;              //there are multiple sons, I assume the flare integral size is 1/2 of elen
	Delta_T=*(elen_poi+iparent_index);              //there are multiple sons, I assume the flare integral size is = of elen
	if(Delta_T < Delta_Tp){
	  //printf("    warning: Delta_T < Delta_Tp, assuming full size\n");	  
	  Delta_T=*(elen_poi+iparent_index);
	}

	fwhm_try=Delta_T;
	rho_estimate=rho_average(Delta_Tp,Delta_T,rho_bkg,rho_p,fwhm_try);
	rho_ratio=rho_estimate/rho_level;
	rho_ratio_withp=rho_level/rho_p;
	//printf("      rho_p %lE, rho_bkg %lE, rho_level %lE, rho est %lE\n",rho_p,rho_bkg,rho_level,rho_estimate);
	//printf("      elen p %lE, elen level %lE\n",Delta_Tp,Delta_T);
	//printf("      rho_ratio: %lf, rho_ratio2 %lf\n",rho_ratio,rho_ratio_withp);
	iinterpol34_h=ipeak_index;
	iinterpol34_l=-1;
	if(rho_ratio_withp <= rho_ratio_ref_34){
	interpol_word34=1;
	  iinterpol34_h=ipeak_index;
	  iinterpol34_l=iparent_index;
	  rho_ratio_withp_h=1.;
	  rho_ratio_withp_l=rho_ratio_withp;
	  //printf("    preparing for interpolation 3/4 %li %li\n",iinterpol34_l,iinterpol34_h);
	  fwhm_time34=interpol(iinterpol34_l,iinterpol34_h,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_l,rho_ratio_withp_h,rho_ratio_ref_34);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_l,rho_ratio_withp_h);
	}else{
	  interpol_word34=-1;
	  iinterpol34_h2=ipeak_index;
	  iinterpol34_h=iparent_index;
	  rho_ratio_withp_h2=1.;
	  rho_ratio_withp_h=rho_ratio_withp;
	  //printf("    preparing for extrapolation 3/4 %li %li (including the peak also)\n",iinterpol34_l,iinterpol34_h);
	  fwhm_time34=interpol(iinterpol34_h,iinterpol34_h2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_h,rho_ratio_withp_h2,rho_ratio_ref_34);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_h,rho_ratio_withp_h2);
	}

      }else{    //if the parent is the root
	interpol_word34=0;
	fwhm_time34=Delta_Tp_time;
      }


    }

    *(Tp_v+ipeak)=Tp_time;
    *(fwhm_v+ipeak)=fwhm_time34;
    *(interpol_word_v+ipeak)=interpol_word34;
    //printf("ipeak %li, tp %lf, fwhm %lf, iw %li\n",ipeak,Tp_time,fwhm_time34,interpol_word34);
  }
  
   return;
}




//****************************************************
void eval_fwhm_v(double rho_bkg,double *rho_ratio_ref_v,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int *howmany_sons_poi,long int nele,long int *ipeak_v,long int *nclean_levels_v,double *Tp_v,double *fwhm_v,long int npeaks, long int *interpol_word_v){

  //same as eval_fwhm, but rho_ratio_ref is a vector, now


  long int ipeak;
  long int ipeak_index;
  long int iparent_index;
  long int iactual_level_index;
  long int nclean_levels;
  long int iclean_level;

  long int parent_level_id;
  long int parent_multiplicity;

  long int interpol_word34;
  double fwhm_time34;
  double Tp_time;
  double Delta_Tp_time;
  double Delta_time_l;
  double Delta_Tp;
  double Delta_T;
  double fwhm_try;
  double rho_p,rho_level;
  double rho_ratio_withp,rho_ratio_withp_l,rho_ratio_withp_l2,rho_ratio_withp_h,rho_ratio_withp_h2;
  long int iinterpol34_l,iinterpol34_l2,iinterpol34_h,iinterpol34_h2;

  double rho_ratio_ref_34;
  for(ipeak=0;ipeak<npeaks;ipeak++){
    interpol_word34=-100;
    nclean_levels=*(nclean_levels_v+ipeak);
    ipeak_index=*(ipeak_v+ipeak);
    Delta_Tp=*(elen_poi+ipeak_index);
    Tp_time=(*(t_start_poi+ipeak_index)+*(t_stop_poi+ipeak_index))/2;
    Delta_Tp_time=(*(t_stop_poi+ipeak_index)-*(t_start_poi+ipeak_index))*(1.+1./(*(Nele_poi+ipeak_index)));
    //printf("\n\nprocessing peak %li, index %li, nclean levels: %li, Tp %lf, Delta_Tp %lf\n",ipeak,ipeak_index,nclean_levels,Tp_time,Delta_Tp_time);
    rho_p=(double)(*(Nele_poi+ipeak_index))/(*(elen_poi+ipeak_index))-rho_bkg;
    rho_ratio_ref_34=*(rho_ratio_ref_v+ipeak);
    //printf("rho_p %lE, elen %lE, size %li\n",rho_p,*(elen_poi+ipeak_index),*(Nele_poi+ipeak_index));
    if(nclean_levels >= 2){   //nclean_levels counts the peak also
      //rho_p=(*elen_poi+ipeak_index)/(*(Nele_poi+ipeak_index))-rho_bkg;
      iactual_level_index=*(parent_id_poi+ipeak_index);
      iinterpol34_l=-1;
      iinterpol34_l2=-1;
      iinterpol34_h=ipeak_index;
      rho_ratio_withp_h=1.;
      iinterpol34_h2=ipeak_index;
      rho_ratio_withp_h2=1.;
      for(iclean_level=0;iclean_level<nclean_levels-1;iclean_level++){   //eval on no more than nclean_levels-1 (nclean_levels counts the peak also)
	//rho_level=(*elen_poi+iactual_level_index)/(*(Nele_poi+iactual_level_index));
	rho_level=(double)(*(Nele_poi+iactual_level_index))/(*(elen_poi+iactual_level_index));
	Delta_T=*(elen_poi+iactual_level_index);
	fwhm_try=Delta_T;
	rho_ratio_withp=rho_level/rho_p;

	//case interpol3/4
	if((rho_ratio_withp >= rho_ratio_ref_34)){
	  iinterpol34_h2=iinterpol34_h;
	  iinterpol34_h=iactual_level_index;
	  rho_ratio_withp_h=rho_ratio_withp;
	}
	if(rho_ratio_withp <= rho_ratio_ref_34){
	  if(iinterpol34_l < 0){
	    iinterpol34_l=iactual_level_index;
	    rho_ratio_withp_l=rho_ratio_withp;
	    //printf("    preparing for interpolation 3/4 %li %li\n",iinterpol34_l,iinterpol34_h);
	  }else{
	    if(iinterpol34_l2 < 0){
	      iinterpol34_l2=iactual_level_index;
	      rho_ratio_withp_l2=rho_ratio_withp;
	    }
	  }
	}

	iactual_level_index=*(parent_id_poi+iactual_level_index);
      }  //end for cicle

      if((iinterpol34_h >= 0) && (iinterpol34_l >= 0)){
	interpol_word34=nclean_levels;
	//printf("    preparing for interpolation 3/4 %li %li %li\n",iinterpol34_l,iinterpol34_h,interpol_word34);
	fwhm_time34=interpol(iinterpol34_l,iinterpol34_h,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_l,rho_ratio_withp_h,rho_ratio_ref_34);
	//printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	//printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_l,rho_ratio_withp_h);
      }else{
	interpol_word34=-nclean_levels;
	//printf("    preparing for extrapolation 3/4 forward %li %li %li %li %li\n",iinterpol34_l,iinterpol34_l2,iinterpol34_h,iinterpol34_h2,interpol_word34);
	  fwhm_time34=interpol(iinterpol34_h,iinterpol34_h2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_h,rho_ratio_withp_h2,rho_ratio_ref_34);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_h,rho_ratio_withp_h2);
      }


    }else{
      iparent_index=*(parent_id_poi+ipeak_index);
      parent_level_id=*(level_id_poi+iparent_index);
      parent_multiplicity=*(howmany_sons_poi+iparent_index);
      //printf("  no clean levels, proceed with the parent\n");
      //printf("    parent level id: %li, parent multiplicity %li\n",parent_level_id,parent_multiplicity);
      if(parent_level_id != -1){
	rho_level=(double)(*(Nele_poi+iparent_index))/(*(elen_poi+iparent_index));
	//Delta_T=*(elen_poi+iparent_index)/2;              //there are multiple sons, I assume the flare integral size is 1/2 of elen
	Delta_T=*(elen_poi+iparent_index);              //there are multiple sons, I assume the flare integral size is = of elen
	if(Delta_T < Delta_Tp){
	  //printf("    warning: Delta_T < Delta_Tp, assuming full size\n");	  
	  Delta_T=*(elen_poi+iparent_index);
	}

	fwhm_try=Delta_T;
	rho_ratio_withp=rho_level/rho_p;
	//printf("      rho_p %lE, rho_bkg %lE, rho_level %lE, rho est %lE\n",rho_p,rho_bkg,rho_level,rho_estimate);
	//printf("      elen p %lE, elen level %lE\n",Delta_Tp,Delta_T);
	//printf("      rho_ratio: %lf, rho_ratio2 %lf\n",rho_ratio,rho_ratio_withp);
	iinterpol34_h=ipeak_index;
	iinterpol34_l=-1;
	if(rho_ratio_withp <= rho_ratio_ref_34){
	interpol_word34=1;
	  iinterpol34_h=ipeak_index;
	  iinterpol34_l=iparent_index;
	  rho_ratio_withp_h=1.;
	  rho_ratio_withp_l=rho_ratio_withp;
	  //printf("    preparing for interpolation 3/4 %li %li\n",iinterpol34_l,iinterpol34_h);
	  fwhm_time34=interpol(iinterpol34_l,iinterpol34_h,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_l,rho_ratio_withp_h,rho_ratio_ref_34);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_l,rho_ratio_withp_h);
	}else{
	  interpol_word34=-1;
	  iinterpol34_h2=ipeak_index;
	  iinterpol34_h=iparent_index;
	  rho_ratio_withp_h2=1.;
	  rho_ratio_withp_h=rho_ratio_withp;
	  //printf("    preparing for extrapolation 3/4 %li %li (including the peak also)\n",iinterpol34_l,iinterpol34_h);
	  fwhm_time34=interpol(iinterpol34_h,iinterpol34_h2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_h,rho_ratio_withp_h2,rho_ratio_ref_34);
	  //printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	  //printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_h,rho_ratio_withp_h2);
	}

      }else{    //if the parent is the root
	interpol_word34=0;
	fwhm_time34=Delta_Tp_time;
      }


    }

    *(Tp_v+ipeak)=Tp_time;
    *(fwhm_v+ipeak)=fwhm_time34;
    *(interpol_word_v+ipeak)=interpol_word34;
    //printf("ipeak %li, tp %lf, fwhm %lf, iw %li\n",ipeak,Tp_time,fwhm_time34,interpol_word34);
  }
  
   return;
}



//***********************************************************
double rho_average(double Delta_Tp,double Delta_T,double rho_bkg,double rho_p,double FWHM){

  //compute the event density for the cluster of effective lenght Delta_t  assuming a trapezoid shape
  //trapezoid is defined with a FWHM and an higher base of width Delta_Tp and event density at the higher
  //base = rho_p
  //the average background density is rho_bkg

  double delta=(Delta_T-Delta_Tp)/2;

  double AA=rho_bkg;
  double BB=rho_p*Delta_Tp/Delta_T;

  double CC2=2.*delta/Delta_Tp;
  double CC3=pow((2.*delta),2.)/(4.*(FWHM-Delta_Tp)*Delta_Tp);

  double CC=1.+CC2-CC3;

   return AA+BB*CC;
}

//***********************************************************
double interpol(long int iinterpol_l,long int iinterpol_h,double *t_start_poi,double *t_stop_poi,long int *Nele_poi,double rho_ratio_l,double rho_ratio_h,double rho_ratio_ref){

  double Delta_t_l=(*(t_stop_poi+iinterpol_l)-*(t_start_poi+iinterpol_l))*(1.+1./(*(Nele_poi+iinterpol_l)));
  double Delta_t_h=(*(t_stop_poi+iinterpol_h)-*(t_start_poi+iinterpol_h))*(1.+1./(*(Nele_poi+iinterpol_l)));;
  double fwhm_time=Delta_t_l+(Delta_t_h-Delta_t_l)*(rho_ratio_ref-rho_ratio_l)/(rho_ratio_h-rho_ratio_l);


  return fwhm_time;
}

//********************************************************************************************************
void eval_fwhm_nomatter_of_nclean(double rho_bkg,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int nele,long int *ipeak_v,double *Tp_v,double *fwhm_v,long int npeaks, long int *interpol_word_v){


  //here I proceed exactly as in eval_fwhm but I do not care about parent multiplicity, I always stop before the root


 
  long int  n_nomatter_levels;
  long int ipeak,ipeak_index;
  double rho_level;
  double rho_ratio_ref_34=0.75;
  //double rho_ratio_ref_34=7./8;
  long int interpol_word34;
  double fwhm_time34;
  double Delta_Tp,Delta_Tp_time;
  double Delta_T;
  double Tp_time,rho_p;
  long int iactual_level_index;
  long int iinterpol34_l,iinterpol34_l2,iinterpol34_h,iinterpol34_h2;
  double rho_ratio_withp,rho_ratio_withp_l,rho_ratio_withp_l2,rho_ratio_withp_h,rho_ratio_withp_h2;

  for(ipeak=0;ipeak<npeaks;ipeak++){
    n_nomatter_levels=1;
    interpol_word34=-100;
    fwhm_time34=-1.E9;
    ipeak_index=*(ipeak_v+ipeak);
    Delta_Tp=*(elen_poi+ipeak_index);
    Tp_time=(*(t_start_poi+ipeak_index)+*(t_stop_poi+ipeak_index))/2;
    Delta_Tp_time=(*(t_stop_poi+ipeak_index)-*(t_start_poi+ipeak_index))*(1.+1./(*(Nele_poi+ipeak_index)));
    //printf("\n\nprocessing peak %li, index %li, Tp %lf, Delta_Tp %lf\n",ipeak,ipeak_index,Tp_time,Delta_Tp_time);
    rho_p=(double)(*(Nele_poi+ipeak_index))/(*(elen_poi+ipeak_index))-rho_bkg;
    //printf("rho_p %lE, elen %lE, size %li\n",rho_p,*(elen_poi+ipeak_index),*(Nele_poi+ipeak_index));


    iinterpol34_l=-1;
    iinterpol34_l2=-1;
    iinterpol34_h=ipeak_index;
    rho_ratio_withp_h=1.;
    iinterpol34_h2=ipeak_index;
    rho_ratio_withp_h2=1.;


    iactual_level_index=*(parent_id_poi+ipeak_index);
    while(iactual_level_index > 0){    //find interpol up to the root (excluded)
      n_nomatter_levels++;

      rho_level=(double)(*(Nele_poi+iactual_level_index))/(*(elen_poi+iactual_level_index));
      Delta_T=*(elen_poi+iactual_level_index);
      rho_ratio_withp=rho_level/rho_p;
      //printf("iactual index %li, ipeak index %li, ipeak %li, rho ratio %lf\n",iactual_level_index,ipeak_index,ipeak,rho_ratio_withp);

      //case interpol3/4
      if((rho_ratio_withp >= rho_ratio_ref_34)){
	iinterpol34_h2=iinterpol34_h;
	iinterpol34_h=iactual_level_index;
	rho_ratio_withp_h=rho_ratio_withp;
      }
      if(rho_ratio_withp <= rho_ratio_ref_34){
	if(iinterpol34_l < 0){
	  iinterpol34_l=iactual_level_index;
	  rho_ratio_withp_l=rho_ratio_withp;
	  //printf("    preparing for interpolation 3/4 %li %li\n",iinterpol34_l,iinterpol34_h);
	}else{
	  if(iinterpol34_l2 < 0){
	    iinterpol34_l2=iactual_level_index;
	    rho_ratio_withp_l2=rho_ratio_withp;
	  }
	}
      }
      

      iactual_level_index=*(parent_id_poi+iactual_level_index);    //next level
    }   //end while cicle


    if(n_nomatter_levels > 1){   //if I have more than one level above root
      if((iinterpol34_h >= 0) && (iinterpol34_l >= 0)){
	interpol_word34=n_nomatter_levels;
	//printf("    preparing for interpolation 3/4 %li %li\n",iinterpol34_l,iinterpol34_h);
	fwhm_time34=interpol(iinterpol34_l,iinterpol34_h,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_l,rho_ratio_withp_h,rho_ratio_ref_34);
	//printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	//printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_l,rho_ratio_withp_h);
      }else{
	interpol_word34=-n_nomatter_levels;
	//printf("    preparing for extrapolation 3/4 forward %li %li %li %li\n",iinterpol34_l,iinterpol34_l2,iinterpol34_h,iinterpol34_h2);
	fwhm_time34=interpol(iinterpol34_h,iinterpol34_h2,t_start_poi,t_stop_poi,Nele_poi,rho_ratio_withp_h,rho_ratio_withp_h2,rho_ratio_ref_34);
	//printf("    >>>fwhm_time: %lf\n",fwhm_time34);
	//printf("    >>>fwhm_time: %lf %lf %lf\n",fwhm_time34,rho_ratio_withp_h,rho_ratio_withp_h2);
      }
    }else{
      interpol_word34=0;
      fwhm_time34=Delta_Tp_time;
    }

    *(Tp_v+ipeak)=Tp_time;
    *(fwhm_v+ipeak)=fwhm_time34;
    *(interpol_word_v+ipeak)=interpol_word34;

  }  //end for loop on peaks
  
  return;
}

//**********************************************************
void eval_fwhm_iterate(double rho_bkg,double *rho_ratio_ref_i_v,long int *level_id_poi,double *elen_poi,long int *Nele_poi,double *t_start_poi,double *t_stop_poi,long int *parent_id_poi,long int *howmany_sons_poi,long int nele,long int *ipeak_v,long int *nclean_levels_v,double *Tp_v,double *fwhm_i_v,long int npeaks, long int *interpol_word_i_v){


    double *Tp_l_v=(double *)mymalloc(sizeof(double)*npeaks);
    double *fwhm_l_v=(double *)mymalloc(sizeof(double)*npeaks);
    long int *interpol_word_l_v=(long int *)mymalloc(sizeof(long int)*npeaks);
    double *rho_ratio_ref_l_v=(double *)mymalloc(sizeof(long int)*npeaks);
    long int ipeak;
    for(ipeak=0;ipeak<npeaks;ipeak++){
      *(rho_ratio_ref_l_v+ipeak)=0.75;
    }
    eval_fwhm_v(rho_bkg,rho_ratio_ref_l_v,level_id_poi,elen_poi,Nele_poi,t_start_poi,t_stop_poi,parent_id_poi,howmany_sons_poi,nele,ipeak_v,nclean_levels_v,Tp_l_v,fwhm_l_v,npeaks,interpol_word_l_v);

    double *Tp_h_v=(double *)mymalloc(sizeof(double)*npeaks);
    double *fwhm_h_v=(double *)mymalloc(sizeof(double)*npeaks);
    long int *interpol_word_h_v=(long int *)mymalloc(sizeof(long int)*npeaks);
    double *rho_ratio_ref_h_v=(double *)mymalloc(sizeof(long int)*npeaks);
    for(ipeak=0;ipeak<npeaks;ipeak++){
      *(rho_ratio_ref_h_v+ipeak)=7./8;
    }
    eval_fwhm_v(rho_bkg,rho_ratio_ref_h_v,level_id_poi,elen_poi,Nele_poi,t_start_poi,t_stop_poi,parent_id_poi,howmany_sons_poi,nele,ipeak_v,nclean_levels_v,Tp_h_v,fwhm_h_v,npeaks,interpol_word_h_v);



    for(ipeak=0;ipeak<npeaks;ipeak++){
      long int index_peak=*(ipeak_v+ipeak);
      double dTp=*(t_stop_poi+index_peak)-*(t_start_poi+index_peak);
      double fwhm_l=*(fwhm_l_v+ipeak);
      double fwhm_h=*(fwhm_h_v+ipeak);
      double elen_peak_2_fwhm_l_ratio=dTp/fwhm_l;
      double elen_peak_2_fwhm_h_ratio=dTp/fwhm_h;
      double elen_peak_2_fwhm_i_ratio=0.5*(elen_peak_2_fwhm_l_ratio+elen_peak_2_fwhm_h_ratio);

      if(elen_peak_2_fwhm_l_ratio > 0.5){
	*(rho_ratio_ref_i_v+ipeak)=0.75+1./4*elen_peak_2_fwhm_l_ratio;
      }else if(elen_peak_2_fwhm_h_ratio < 0.1){
	*(rho_ratio_ref_i_v+ipeak)=0.75+1./4*elen_peak_2_fwhm_h_ratio;	
      }else{
	*(rho_ratio_ref_i_v+ipeak)=0.75+1./4*elen_peak_2_fwhm_i_ratio;	
      }
     
    }

    double *Tp_i_v=(double *)mymalloc(sizeof(double)*npeaks);
    eval_fwhm_v(rho_bkg,rho_ratio_ref_i_v,level_id_poi,elen_poi,Nele_poi,t_start_poi,t_stop_poi,parent_id_poi,howmany_sons_poi,nele,ipeak_v,nclean_levels_v,Tp_i_v,fwhm_i_v,npeaks,interpol_word_i_v);


//    long int iloop;
//    long int nloop=100;
//    double *fwhm_i_old_v=(double *)mymalloc(sizeof(double)*npeaks);
//    for(iloop=0;iloop<nloop;iloop++){
//      for(ipeak=0;ipeak<npeaks;ipeak++){
//	long int index_peak=*(ipeak_v+ipeak);
//	double dTp=*(t_stop_poi+index_peak)-*(t_start_poi+index_peak);
//	double fwhm_i=*(fwhm_i_v+ipeak);
//	double elen_peak_2_fwhm_i_ratio=dTp/fwhm_i;
//	*(rho_ratio_ref_i_v+ipeak)=0.75+1./4*elen_peak_2_fwhm_i_ratio;     
//	*(fwhm_i_old_v+ipeak)=fwhm_i; 
//	printf("peak: %li, fwhm: %lE\n",ipeak,fwhm_i);
//     }
//      eval_fwhm_v(rho_bkg,rho_ratio_ref_i_v,level_id_poi,elen_poi,Nele_poi,t_start_poi,t_stop_poi,parent_id_poi,howmany_sons_poi,nele,ipeak_v,nclean_levels_v,Tp_i_v,fwhm_i_v,npeaks,interpol_word_i_v);
//    }
//    double dfwhm_rel;
//    for(ipeak=0;ipeak<npeaks;ipeak++){
//      dfwhm_rel=abs(*(fwhm_i_v+ipeak)/(*(fwhm_i_old_v+ipeak))-1);
//      printf("peak: %li, rel error on fwhm: %lE\n",ipeak,dfwhm_rel);
//    }

    return;
}


//****************************************************************************
FILE *init_survived_cluster_file(char *clusterfile,char *toutformat,char *xoutformat,long int N_tol,double prob_thr){

  //open cluster file
  
  printf("opening survived clusters list file: %s\n",clusterfile);

  FILE *cfile=fopen(clusterfile,"w");
  if(cfile == NULL){
    fprintf(stderr,"unable to open file for survived clusters, bye\n");
    exit(-1);
  }

  long int nclusters_dummy=-1;
  fprintf(cfile,"Nclusters %8li\n",nclusters_dummy);
  fprintf(cfile,"N_tol %8li\n",N_tol);
  fprintf(cfile,"P thr %10.3lE\n",prob_thr);
  fprintf(cfile,"t_format %s\n",toutformat);
  fprintf(cfile,"x_format %s\n",xoutformat);
  //file columns:
  //level id of \Delta_{thr}   (-1 for root)
  //  t_start           
  //  t_stop
  //  cluster size
  //  cluster len
  //  cluster effective length    
  fprintf(cfile,"col1: \"level id\" (-1 is for the root)\n");
  fprintf(cfile,"col2: \"t start\"\n");
  fprintf(cfile,"col3: \"t stop\"\n");
  fprintf(cfile,"col4: \"cluster size\"\n");
  fprintf(cfile,"col5: \"cluster length\" (in X domain if \"twocolumns\" mode has been selected)\n");
  fprintf(cfile,"col6: \"effective cluster length\" (in X domain if \"twocolumns\" mode has been selected)\n");  
  fprintf(cfile,"col7: \"event density\" (in X domain if \"twocolumns\" mode has been selected)\n");  
  fprintf(cfile,"col8: \"cluster coincidence probability in the null hypothesis that the events within its own parent are uniformly distributed\"\n");  

  return cfile;
}

//************************************
int fprintf_cluster(FILE *cfile,long int level_id,double t_start,double t_stop,long int Nele,double len,double elen,double chance_prob,char *toutformat,char *xoutformat){
  char format[200];
  double rho=(double)Nele/elen;

  strcpy(format,"%8li ");
  strcat(format,toutformat);
  strcat(format," ");
  strcat(format,toutformat);
  strcat(format," %8li ");
  strcat(format,xoutformat);
  strcat(format," ");
  strcat(format,xoutformat);
  strcat(format," %10.3lE %8.2lE\n");

  int res=fprintf(cfile,format,level_id,t_start,t_stop,Nele,len,elen,rho);
  return res;
}

//*************************************
void close_survived_cluster_file(FILE *cfile,long int nclusters){
  
  //write ther true number of clusters at the head of the cluster file
  fseek(cfile,0,SEEK_SET);
  fprintf(cfile,"Nclusters %8li\n",nclusters);
  fseek(cfile,0,SEEK_END);
  fclose(cfile);  
  return;
}

//****************************************************************************
FILE *init_peaklist_file(char *peaklist_file,char *toutformat,char *xoutformat,long int N_tol,double prob_thr,long int Npeaks){

  //open peaklist file

  printf("opening peaklist file: %s\n",peaklist_file);
  FILE *pfile=fopen(peaklist_file,"w");
  if(pfile == NULL){
    fprintf(stderr,"unable to open file for flare list, bye\n");
    exit(-1);
  }

  fprintf(pfile,"Npeaks %8li\n",Npeaks);
  fprintf(pfile,"N_tol %8li\n",N_tol);
  fprintf(pfile,"P thr %10.3lE\n",prob_thr);
  fprintf(pfile,"t_format %s\n",toutformat);
  fprintf(pfile,"x_format %s\n",xoutformat);
  fprintf(pfile,"col1: \"index of the peak cluster within the list of survived clusters\"\n");
  fprintf(pfile,"col2: \"t peak\"\n");
  fprintf(pfile,"col3: \"peak cluster size\"\n");
  fprintf(pfile,"col4: \"effective peak cluster length\" (in t domain)\n");  
  fprintf(pfile,"col5: \"event density at peak\" (in X domain if \"twocolumns\" mode has been selected)\n");  
  fprintf(pfile,"col6: \"effective peak cluster length\" (in X domain if \"twocolumns\" mode has been selected)\n");
  fprintf(pfile,"col7: \"flare FWHM\" (in t domain)\n");  
  fprintf(pfile,"col8: \"FWHM_interpol_word\" (>0: interpolation, <0: extrapolation, 0: the peak is just over the root)\n");
  fprintf(pfile,"       the absolute value of FWHM_interpol_word is the number of clusters encountered descending the tree,\n");
  fprintf(pfile,"       starting from the peak up to the first cluster with more than one son (and excluding this cluster)\n");
  fprintf(pfile,"col9: \"Peak cluster height\", e.g., the number of levels from root to the peak (exluding the root)\n");

  return pfile;
}

//*****************************************
void *mymalloc(size_t size){
  void *mypoi=malloc(size);
  if(mypoi == NULL){
    fprintf(stderr,"unable to get memory access\n");
    exit(-1);
  }

  return mypoi;
}

//*******************************
void help(void){

  printf("rmrndm usage:\n");
  printf("./rmrndm <input file with candidate cluster list> <output filename for the list of survived clusters>\n");
  printf("\n");
  printf("optional arguments: \n");
  printf("                 MAXDEV\n");
  printf("                 NOCORR\n");
  printf("                 THR_PROB <threshold probability>\n");
  printf("                 NSIGMA <nsigma>\n");
  printf("                 SCANTABLEDIR <path to the stab.txt and ctab.txt files>\n");
  printf("                 PEAKLIST_FILE <output filename with the list of flares>\n");
  printf("Each cluster is put in the list of survived clusters when the probability\n");
  printf("for it to form (in the null hypothesis that the events of its parent are uniformly distributed) is\n");
  printf("lower than a threshold probability (threshold probability = 1 - c.l.).\n");
  printf("NB: threshold probability can be set directly (THR_PROB parameter)\n");
  printf("    or in standard gaussian units (NSIGMA parameter, with default value= 3)\n");
  printf("\n");
  printf("When rmrndm works in MAXDEV mode, only the most relevant cluster (with respect to root) is reported\n");
  printf("When NOCORR is set, probabilities are evaluated using simple scan-statistics\n");
  printf("NB: NOCORR and MAXDEV must be used and must be set together to prepare a new correction table\n");
  return;
}