diff --git a/doc/formats.pdf b/doc/formats.pdf
index d8578e6..6d381d4 100644
Binary files a/doc/formats.pdf and b/doc/formats.pdf differ
diff --git a/doc/formats.tex b/doc/formats.tex
index ac9b569..aac2bbe 100644
--- a/doc/formats.tex
+++ b/doc/formats.tex
@@ -1,9 +1,9 @@
 \documentclass[10pt]{article}
-\usepackage{color}
+\usepackage{xcolor}
 \definecolor{gray}{rgb}{0.7,0.7,0.7}
 \usepackage{framed}
 \usepackage{enumitem}
-\usepackage{longtable}
+\usepackage{adjustbox}
 \usepackage[pdfborder={0 0 0},hyperfootnotes=false]{hyperref}
 
 \addtolength{\textwidth}{3.4cm}
@@ -243,12 +243,67 @@ \section{.rlens.gz}
 Readlength distribution. Distribution is count of alignments of
 specific readlengths. Depending on runmode there might be multiple
 groups as taxid/refs and there will be group specific
-distributions. Distributions for different groups are split by
+distributions.
+\subsection{Old versions (pre april 2026)}
+Distributions for different groups are split by
 newlines. First entry on each line is the identifer
 (chromosomename,taxid). The remaining entries are the number of times
 we have observed an alignment of length (columnnumber). Notice that
 the first 30 column of counts is likely to be zero since reads shorter
 than 30basepairs are normally discarded.
+\subsection{Newer versions: human nice format}
+Each distribution is a newline. Each entry in the distribution is bin:count seperated by tab.
+Example can be seen below. This is the default option that should be enforcable by \colorbox{gray!10}{\texttt{--rlens\_flat\_out 0}}. Output below is clipped at rightside and at bottom. Example shown is from one of the testfiles, and has been generated by \colorbox{gray!10}{\texttt{cut -c1-80 test/output/test\_getdamage\_local.rlens|head }}
+\begin{framed}
+  \small
+\begin{verbatim}
+id	rlen:count
+0	35:33	36:29	37:24	38:38	39:24	40:29	41:36	42:28	43:27	44:36	45:29	46:23	47:32...
+1	35:103	36:106	37:102	38:111	39:109	40:126	41:125	42:164	43:154	44:133	45:152...
+2	35:35	36:44	37:28	38:41	39:33	40:36	41:41	42:33	43:34	44:43	45:41	46:36	47:40...
+3	35:36	36:39	37:25	38:37	39:29	40:34	41:39	42:33	43:34	44:43	45:41	46:33	47:39...
+4	35:100	36:98	37:95	38:112	39:107	40:136	41:127	42:161	43:157	44:138	45:134	46:...
+5	35:92	36:118	37:110	38:125	39:129	40:168	41:158	42:164	43:195	44:180	45:177	46...
+6	35:95	36:130	37:111	38:125	39:129	40:178	41:167	42:148	43:196	44:193	45:185	46...
+7	35:138	36:145	37:131	38:125	39:147	40:179	41:190	42:205	43:196	44:201	45:220	4...
+8	35:95	36:132	37:113	38:137	39:133	40:182	41:162	42:176	43:196	44:190	45:189	46...
+...
+\end{verbatim}
+\end{framed}
+\subsection{Newer versions: human nice format}
+
+
+The parsing of the standard output rlens file can be abit confusing for programs to parse due to the varying columnwidth. Therefore by supplying  \colorbox{gray!10}{\texttt{--rlens\_flat\_out 1}} or \colorbox{gray!10}{\texttt{-z 1}} in relevant commands the output will be flattened in to a simpler but longer format.
+\begin{framed}
+  \small
+\begin{verbatim}
+IFILE=./data/f570b1db7c.dedup.filtered.rname.bam
+../metaDMG-cpp getdamage -r 1 -l 35 -p 5 -o output/test_getdamage_local ${IFILE} -z 1
+gunzip -c output/test_getdamage_local.rlens.gz|head -n20
+##below is the output from the last command above
+ID	rlen	count
+0	35	33
+0	36	29
+0	37	24
+0	38	38
+0	39	24
+0	40	29
+0	41	36
+0	42	28
+0	43	27
+0	44	36
+0	45	29
+0	46	23
+0	47	32
+0	48	23
+0	49	25
+0	50	17
+0	51	16
+0	52	19
+0	53	17
+\end{verbatim}
+\end{framed}
+
 \end{document}
 
 %%% Local Variables:
diff --git a/main_dfit.cpp b/main_dfit.cpp
index 47b050c..c1a8103 100644
--- a/main_dfit.cpp
+++ b/main_dfit.cpp
@@ -25,7 +25,7 @@
 #include "dfit_helppage.h"
 #include "mrand.h"
 
-extern htsFormat *dingding2;
+//extern htsFormat *dingding2;
 
 // the compare function for double values
 static int compare (const void * a, const void * b)
@@ -672,6 +672,7 @@ int main_dfit(int argc, char **argv) {
     double CI = 0.95;
     int doCI = 2;
     int rng_type = -1;
+    htsFormat *dingding2 = (htsFormat *)calloc(1, sizeof(htsFormat));
     while (*(++argv)) {
         if (strcasecmp("-h", *argv) == 0)
           HelpPageSimple(stderr);
@@ -947,7 +948,7 @@ int main_dfit(int argc, char **argv) {
       free(outfile_name);
     if(lib_prep)
       free(lib_prep);
-    
+    free(dingding2);
   gettimeofday(&end_time, NULL);
   long seconds = end_time.tv_sec - start_time.tv_sec;
   long microseconds = end_time.tv_usec - start_time.tv_usec;
diff --git a/main_pmd.cpp b/main_pmd.cpp
index d69ffe9..7caef00 100644
--- a/main_pmd.cpp
+++ b/main_pmd.cpp
@@ -15,8 +15,8 @@
 
 #include "profile.h"  // for reconstructRefWithPosHTS
 
-extern htsFormat *dingding2;
-htsFormat *dingding3 = dingding2;
+//extern htsFormat *dingding2;
+//htsFormat *dingding3 = dingding2;
 
 
 int nproc = 0;  // number of reads processed
@@ -159,7 +159,7 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on
     kstr->s = NULL;
     mypair.first = kstr;
     samFile *in = NULL;
-
+    htsFormat *dingding3 = (htsFormat *)calloc(1, sizeof(htsFormat));
     if (refName) {
         char *ref = (char *)malloc(10 + strlen(refName) + 1);
         snprintf(ref,10 + strlen(refName) + 1, "reference=%s", refName);
@@ -241,6 +241,7 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on
     sam_hdr_destroy(hdr);
     sam_close(in);
     free(fname);
+    free(dingding3);
 }
 
 int usage(FILE *fp, int val) {
diff --git a/metaDMG.cpp b/metaDMG.cpp
index 69779c5..e4917ec 100644
--- a/metaDMG.cpp
+++ b/metaDMG.cpp
@@ -29,7 +29,7 @@
 #include "regression.h"  // for main_regression
 #endif
 
-htsFormat *dingding2 = (htsFormat *)calloc(1, sizeof(htsFormat));
+
 typedef std::map<int, char *> int2char;
 int usage_getdamage(FILE *fp) {
     fprintf(fp, "\nUsage: metadamage getdamage [options] <in.bam>|<in.sam>|<in.cram>\n");
@@ -41,6 +41,7 @@ int usage_getdamage(FILE *fp) {
     fprintf(fp, "  -r/--run_mode\t\t 0: global estimate (default)\n\t\t\t 1: damage patterns will be calculated for each chr/scaffold contig\n");
     fprintf(fp, "  -i/--ignore_errors\t continue analyses even if there are errors.\n");
     fprintf(fp, "  -o/--out_prefix\t output prefix (default: meta)\n");
+    fprintf(fp, "  -z/--rlens_flat_out\t make flat output of bins. Nice for computers\n");
     return 1;
 }
 
@@ -88,7 +89,12 @@ std::map<int,mydataD> getval_full_norec(std::map<int, mydataD> &retmap, int2int
   
   //funky modern syntax below
   //loop over all entries in retmap. lizard king 
-  for (const auto &[taxid, data] : retmap) {
+  //  for (const auto &[taxid, data] : retmap) {
+
+  for(  std::map<int, mydataD>::iterator it=retmap.begin();it!=retmap.end();it++){
+    int taxid = it->first;
+    mydataD data = it->second;
+    
     int current = taxid;
     //   fprintf(stderr,"taxid: %d\n",taxid);
      while (true) {
@@ -130,8 +136,8 @@ int mywarn =1;
 
 mydata2 getval_stats_unknownfunctionality(std::map<int, mydata2> &retmap, int2intvec &child, int taxid) {
   if(mywarn>0){
-  fprintf(stderr,"PAS PAA SATAN\n");
-  mywarn--;
+    fprintf(stderr,"PAS PAA SATAN\n");
+    mywarn--;
   }
     // fprintf(stderr,"getval\t%d\t%d\n",taxid,howmany);
     std::map<int, mydata2>::iterator it = retmap.find(taxid);
@@ -193,9 +199,12 @@ int main_getdamage(int argc, char **argv) {
     char *onam = strdup("meta");
     int nthreads = 4;
     int ignore_errors = 0;
+    int rlens_flat_out = 0;
+    htsFormat *dingding2 = (htsFormat *)calloc(1, sizeof(htsFormat));
     // fix thesepro
     static struct option lopts[] = {
         {"threads", required_argument, 0, 'n'},
+	{"rlens_flat_out", required_argument, 0, 'z'},
         {"fasta", required_argument, 0, 'f'},
         {"min_length", required_argument, 0, 'l'},
         {"print_length", required_argument, 0, 'p'},
@@ -207,15 +216,21 @@ int main_getdamage(int argc, char **argv) {
 
     int c;
     while ((c = getopt_long(argc, argv,
-                            "n:f:l:p:r:o:h",
+                            "iz:n:f:l:p:r:o:h",
                             lopts, NULL)) >= 0) {
       switch (c) {
+      case 'i':
+        ignore_errors = 1;
+        break;
       case 'n':
 	nthreads = atoi(optarg);
 	break;
       case 'f':
 	refName = strdup(optarg);
 	break;
+      case 'z':
+	rlens_flat_out = atoi(optarg);
+	break;
       case 'l':
 	minLength = atoi(optarg);
 	break;
@@ -238,7 +253,7 @@ int main_getdamage(int argc, char **argv) {
     }
     if (optind < argc)
         fname = strdup(argv[optind]);
-    fprintf(stderr, "\t-> ./metaDMG-cpp refName: %s min_length: %d print_length: %d run_mode: %d out_prefix: %s nthreads: %d ignore_errors: %d\n", refName, minLength, printLength, runmode, onam, nthreads, ignore_errors);
+    fprintf(stderr, "\t-> ./metaDMG-cpp refName: %s min_length: %d print_length: %d run_mode: %d out_prefix: %s nthreads: %d ignore_errors: %d rlens_flat_out: %d\n", refName, minLength, printLength, runmode, onam, nthreads, ignore_errors,rlens_flat_out);
     if (fname == NULL) {
         usage_getdamage(stderr);
         return 0;
@@ -253,11 +268,15 @@ int main_getdamage(int argc, char **argv) {
 
     if ((fp = sam_open_format(fname, "r", dingding2)) == NULL) {
         fprintf(stderr, "[%s] nonexistant file: %s\n", __FUNCTION__, fname);
-        exit(0);
+        exit(1);
     }
 
     bam1_t *b = bam_init1();
     bam_hdr_t *hdr = sam_hdr_read(fp);
+    if(hdr==NULL){
+      fprintf(stderr,"\t-> Hello Doctor! im afraid there is an error reading the header\n");
+      exit(1);
+    }
     int checkIfSorted(char *str);
     if(checkIfSorted(hdr->text)){
       fprintf(stderr, "Input alignment file is not sorted.");
@@ -292,7 +311,7 @@ int main_getdamage(int argc, char **argv) {
         float mylen = b->core.l_qseq;
 
         int whichref = 0;
-        if (runmode == 1)
+        if (runmode == 1)//<- runmode 0 means local one means global
             whichref = b->core.tid;
         std::map<int, std::vector<float> >::iterator it = gcconts.find(whichref);
         if (it == gcconts.end()) {
@@ -314,7 +333,7 @@ int main_getdamage(int argc, char **argv) {
 
     dmg->printit(stdout, printLength);
     dmg->write(onam, runmode == 1 ? hdr : NULL);
-    dmg->bwrite(onam);
+    dmg->bwrite(onam,rlens_flat_out);
 
     // write stat
     char buf[1024];
@@ -331,8 +350,13 @@ int main_getdamage(int argc, char **argv) {
         assert(it3 != seqlens.end());
         if (0)
             gzprintf(fpstat, "%d\t%lu\t%f\t%f\t%f\t%f\tNA\tNA\n", it->first, it2->second.nreads, mean(it3->second), var(it3->second), mean(it->second), var(it->second));
-        else
-            gzprintf(fpstat, "%s\t%lu\t%f\t%f\t%f\t%f\tNA\tNA\n", sam_hdr_tid2name(hdr, it->first), it2->second.nreads, mean(it3->second), var(it3->second), mean(it->second), var(it->second));
+        else{
+	  if(runmode==1)
+	    gzprintf(fpstat, "%s\t%lu\t%f\t%f\t%f\t%f\tNA\tNA\n", sam_hdr_tid2name(hdr, it->first), it2->second.nreads, mean(it3->second), var(it3->second), mean(it->second), var(it->second));
+	  else
+	    gzprintf(fpstat, "global\t%lu\t%f\t%f\t%f\t%f\tNA\tNA\n", it2->second.nreads, mean(it3->second), var(it3->second), mean(it->second), var(it->second));
+	    
+	}
     }
     gzclose(fpstat);
 
@@ -342,26 +366,13 @@ int main_getdamage(int argc, char **argv) {
     destroy_damage(dmg);
     free(fname);
     free(onam);
-    return 0;
-}
-
-int main_index(int argc, char **argv) {
-    char *infile = argv[1];
-    fprintf(stderr, "infile: %s\n", infile);
-    char onam[strlen(infile) + 20];
-    snprintf(onam,strlen(infile) + 20, "%s.idx", infile);
-    fprintf(stderr, "outfile: %s\n", onam);
-    FILE *fp = NULL;
-    if (((fp = fopen(onam, "wb"))) == NULL) {
-        fprintf(stderr, "Problem opening file\n");
-        return 0;
-    }
-
-    fclose(fp);
+    if(refName) free(refName);
+    free(dingding2);
     return 0;
 }
 
 int main_print(int argc, char **argv) {
+  //  fprintf(stderr,"\nYOYOYOYOYYOOYOOY\n");
     if (argc == 1) {
         fprintf(stderr, "./metaDMG-cpp print file.bdamage.gz [-names names.gz -bam file.bam -ctga -countout -nodes -howmany -r -doOld -nodes]\n");
         return 0;
@@ -416,7 +427,7 @@ int main_print(int argc, char **argv) {
     if (search != -1 && doold == 0) {
         std::map<int, double *> retmap = load_bdamage3(infile, howmany);
         double *dbl = getval(retmap, child, search, howmany);
-        double dbldbl[3 * howmany + 1];  // 3 because ct,ga,other
+        double *dbldbl = new double[3 * howmany + 1];  // 3 because ct,ga,other
         dbldbl[0] = dbl[0];
         for (int i = 0; i < 3 * howmany; i++)
             dbldbl[i + 1] = dbl[1 + i] / dbl[0];
@@ -425,6 +436,7 @@ int main_print(int argc, char **argv) {
         for (int i = 0; i < 3 * howmany; i++)
             fprintf(stdout, "\t%f", dbldbl[1 + i]);
         fprintf(stdout, "\n");
+	delete [] dbldbl;
         return 0;
     }
 
@@ -473,7 +485,7 @@ int main_print(int argc, char **argv) {
         assert(nread == 2 * sizeof(int));
         for (int at = 0; at < printlength; at++) {
             assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16));
-            if (at > howmany)
+            if (at >= howmany)
                 continue;
             if (search == -1 || search == ref_nreads[0]) {
                 if (ctga == 0) {
@@ -527,8 +539,8 @@ int main_print(int argc, char **argv) {
         }
 
         for (int at = 0; at < printlength; at++) {
-            assert(16 * sizeof(int) == bgzf_read(bgfp, data, sizeof(float) * 16));
-            if (at > howmany)
+            assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16));
+            if (at >= howmany)
                 continue;
             if (search == -1 || search == ref_nreads[0]) {
                 if (ctga == 0) {
@@ -654,7 +666,7 @@ int main_print2(int argc, char **argv) {
     if (search != -1 && doold == 0) {
         std::map<int, double *> retmap = load_bdamage3(infile, howmany);
         double *dbl = getval(retmap, child, search, howmany);
-        double dbldbl[3 * howmany + 1];  // 3 because ct,ga,other
+        double *dbldbl = new double [3 * howmany + 1];  // 3 because ct,ga,other
         dbldbl[0] = dbl[0];
         for (int i = 0; i < 3 * howmany; i++)
             dbldbl[i + 1] = dbl[1 + i] / dbl[0];
@@ -663,6 +675,7 @@ int main_print2(int argc, char **argv) {
         for (int i = 0; i < 3 * howmany; i++)
             fprintf(stdout, "\t%f", dbldbl[1 + i]);
         fprintf(stdout, "\n");
+	delete [] dbldbl;
         return 0;
     }
 
@@ -710,7 +723,7 @@ int main_print2(int argc, char **argv) {
         fprintf(stdout, "\n");
     }
 
-    int data[16];
+    float data[16];
 
     while (1) {
         int nread = bgzf_read(bgfp, ref_nreads, 2 * sizeof(int));
@@ -720,7 +733,7 @@ int main_print2(int argc, char **argv) {
         fprintf(stderr, "ref: %d nreads: %d\n", ref_nreads[0], ref_nreads[1]);
         assert(nread == 2 * sizeof(int));
         for (int at = 0; at < printlength; at++) {
-            assert(16 * sizeof(int) == bgzf_read(bgfp, data, sizeof(int) * 16));
+            assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16));
             if ((at + 1) > howmany)
                 continue;
             if (search == -1 || search == ref_nreads[0]) {
@@ -742,7 +755,7 @@ int main_print2(int argc, char **argv) {
                 }
                 if (countout == 1) {
                     for (int i = 0; i < 16; i++)
-                        fprintf(stdout, "\t%d", data[i]);
+                        fprintf(stdout, "\t%f", data[i]);
                     fprintf(stdout, "\n");
                 } else {
                     float flt[16];
@@ -775,7 +788,7 @@ int main_print2(int argc, char **argv) {
         }
 
         for (int at = 0; at < printlength; at++) {
-            assert(16 * sizeof(int) == bgzf_read(bgfp, data, sizeof(int) * 16));
+            assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16));
             if (at + 1 > howmany)
                 continue;
             if (search == -1 || search == ref_nreads[0]) {
@@ -794,7 +807,7 @@ int main_print2(int argc, char **argv) {
                 }
                 if (countout == 1) {
                     for (int i = 0; i < 16; i++)
-                        fprintf(stdout, "\t%d", data[i]);
+                        fprintf(stdout, "\t%f", data[i]);
                     fprintf(stdout, "\n");
                 } else {
                     float flt[16];
@@ -959,23 +972,33 @@ int main_merge(int argc, char **argv) {
 int2int getlcadist(char *fname) {
     //  fprintf(stderr,"fname: %s\n",fname);
     int2int lcadist;
-    int tlen = strlen(fname) + 10;
-    char tmp[tlen];
-    snprintf(tmp, tlen, "%s.stat", fname);
+    if(strlen(fname)>1000){
+      fprintf(stderr,"\t-> Ridiculus long filename: \'%s\'\n",fname);
+      exit(1);
+    }
+    char tmp[1024];
+    snprintf(tmp, 1000, "%s.stat", fname);
     fprintf(stderr, "tmp: %s\n", tmp);
     FILE *fp = NULL;
     fp = fopen(tmp, "rb");
     assert(fp != NULL);
     char buf[4096];
     while (fgets(buf, 4096, fp)) {
-        int key = atoi(strtok(buf, "\t\n "));
-        int val = atoi(strtok(NULL, "\t\n "));
-        lcadist[key] = val;
+      char *tok1 = strtok(buf, "\t\n ");
+      char *tok2 = strtok(NULL, "\t\n ");
+      if(!tok1||!tok2)
+	continue;
+      int key = atoi(tok1);
+      int val = atoi(tok2);
+      lcadist[key] = val;
     }
     fprintf(stderr, "\t-> Done reading: %lu entries from file: \'%s\'\n", lcadist.size(), tmp);
+    fclose(fp);
     return lcadist;
 }
 
+
+
 std::map<int, double *> getcsv(char *fname) {
     std::map<int, double *> ret;
     FILE *fp = NULL;
@@ -1095,20 +1118,18 @@ int main_merge2(int argc, char **argv) {
 }
 
 int main_print_all(int argc, char **argv) {
-    fprintf(stderr, "./metaDMG-cpp print_all file.bdamage.gz -names file.gz [-howmany 5] -nodes trestructure.gz -names fil.gz\n");
+    fprintf(stderr, "./metaDMG-cpp print_all file.bdamage.gz -names file.gz -nodes trestructure.gz -names fil.gz\n");
     if (argc <= 2)
         return 0;
     char *infile_bdamage = NULL;
     char *infile_nodes = NULL;
     char *infile_names = NULL;
-    int howmany = 5;
+
 
     while (*(++argv)) {
         if (strcasecmp("-names", *argv) == 0)
             infile_names = strdup(*(++argv));
-        else if (strcasecmp("-howmany", *argv) == 0)
-            howmany = atoi(*(++argv));
-        else if (strcasecmp("-nodes", *argv) == 0)
+	else if (strcasecmp("-nodes", *argv) == 0)
             infile_nodes = strdup(*(++argv));
         else
             infile_bdamage = strdup(*argv);
@@ -1125,15 +1146,15 @@ int main_print_all(int argc, char **argv) {
 
     if (infile_nodes != NULL)
         parse_nodes(infile_nodes, rank, parent, child, 1);
-
+    int howmany = 5;
     std::map<int, double *> retmap = load_bdamage3(infile_bdamage, howmany);
     fprintf(stderr, "\t-> number of entries in damage pattern file: %lu\n", retmap.size());
     int2char name = parse_names(infile_names);
 
     BGZF *bgfp = NULL;
-    samFile *samfp = NULL;
-
+    
     float presize = retmap.size();
+
     getval(retmap, child, 1, howmany);  // this will do everything
     float postsize = retmap.size();
     fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize);
@@ -1167,6 +1188,7 @@ int main_print_all(int argc, char **argv) {
 }
 
 int main_print_ugly(int argc, char **argv) {
+  htsFormat *dingding2 = (htsFormat *)calloc(1, sizeof(htsFormat));
   if(1)
     fprintf(stderr,"\t-> print_ugly functionality will be removed\n");
   
@@ -1247,7 +1269,7 @@ int main_print_ugly(int argc, char **argv) {
 
     for (std::map<int, mydataD>::iterator it = retmap.begin(); it != retmap.end(); it++) {
         int taxid = it->first;
-        mydataD md = it->second;
+        //mydataD md = it->second;
         if (it->second.nal == 0)
             continue;
         /*
@@ -1358,6 +1380,7 @@ int main_print_ugly(int argc, char **argv) {
       free(infile_lcastat);
     free(kstr->s);
     delete kstr;
+    free(dingding2);
     return 0;
 }
 
@@ -1406,8 +1429,6 @@ int main(int argc, char **argv) {
         return main_pmd(argc, argv);
     if (!strcmp(argv[0], "getdamage"))
         main_getdamage(argc, argv);
-    if (!strcmp(argv[0], "index"))
-        main_index(argc, argv);
     if (!strcmp(argv[0], "print"))
         main_print(argc, argv);
     if (!strcmp(argv[0], "print_all"))
@@ -1428,8 +1449,8 @@ int main(int argc, char **argv) {
         main_lca(argc, argv);
     if (!strcmp(argv[0], "mergedamage"))
       main_mergedamage(argc, argv);
-    free(dingding2);
-    //  fprintf(stderr, "\t[ALL done] cpu-time used =  %.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
-    // fprintf(stderr, "\t[ALL done] walltime used =  %.2f sec\n", (float)(time(NULL) - t2));
+
+    fprintf(stderr, "\t[ALL done] cpu-time used =  %.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
+    fprintf(stderr, "\t[ALL done] walltime used =  %.2f sec\n", (float)(time(NULL) - t2));
     return 0;
 }
diff --git a/ngsLCA.cpp b/ngsLCA.cpp
index 69605ee..1dadda4 100644
--- a/ngsLCA.cpp
+++ b/ngsLCA.cpp
@@ -406,8 +406,8 @@ std::vector<int> purge(std::vector<int> &taxids, std::vector<int> &editdist) {
     return tmpnewvec;
 }
 
-void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hdr, int2char &rank, int2char &name_map, int minmapq, int discard, int editMin, int editMax, double scoreLow, double scoreHigh, int minlength, int lca_rank, char *prefix, int howmany, samFile *fp_usedreads, int skipnorank, int2int &rank2level, int nthreads, int weighttype,long maxreads,samFile *fp_famout) {
-  fprintf(stderr, "[%s] \t-> editMin:%d editmMax:%d scoreLow:%f scoreHigh:%f minlength:%d discard: %d prefix: %s howmany: %d skipnorank: %d weighttype: %d maxreads: %ld\n", __FUNCTION__, editMin, editMax, scoreLow, scoreHigh, minlength, discard, prefix, howmany, skipnorank, weighttype,maxreads);
+void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hdr, int2char &rank, int2char &name_map, int minmapq, int discard, int editMin, int editMax, double scoreLow, double scoreHigh, int minlength, int lca_rank, char *prefix, int howmany, samFile *fp_usedreads, int skipnorank, int2int &rank2level, int nthreads, int weighttype,long maxreads,samFile *fp_famout,int rlens_flat_out) {
+  fprintf(stderr, "[%s] \t-> editMin:%d editmMax:%d scoreLow:%f scoreHigh:%f minlength:%d discard: %d prefix: %s howmany: %d skipnorank: %d weighttype: %d maxreads: %ld rlens_flat_out: %d\n", __FUNCTION__, editMin, editMax, scoreLow, scoreHigh, minlength, discard, prefix, howmany, skipnorank, weighttype,maxreads,rlens_flat_out);
     assert(fp_in != NULL);
     damage *dmg = new damage(howmany, nthreads, 13);
     bam1_t *aln = bam_init1();  // initialize an alignment
@@ -645,7 +645,7 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd
         }
     }
     fprintf(stderr,"maxreads: %ld nreads: %ld\n",maxreads,nreads);
-    dmg->bwrite(prefix);
+    dmg->bwrite(prefix,rlens_flat_out);
 
     specs.clear();
     editdist.clear();
@@ -794,7 +794,7 @@ int main_lca(int argc, char **argv) {
 	fprintf(stderr, "writing headers to %s", p->famout_sam);
     }
     gzprintf(p->fp1,"queryid\tseq\tlen\tnaln\tgc\tlca\ttaxa_path\n");
-    hts(p->fp1, p->hts, *i2i, parent, p->header, rank, name_map, p->minmapq, p->discard, p->editdistMin, p->editdistMax, p->simscoreLow, p->simscoreHigh, p->minlength, lca_rank, p->outnames, p->howmany, usedreads_sam, p->skipnorank, tax2level, p->nthreads, p->weighttype,p->maxreads,famout_sam);
+    hts(p->fp1, p->hts, *i2i, parent, p->header, rank, name_map, p->minmapq, p->discard, p->editdistMin, p->editdistMax, p->simscoreLow, p->simscoreHigh, p->minlength, lca_rank, p->outnames, p->howmany, usedreads_sam, p->skipnorank, tax2level, p->nthreads, p->weighttype,p->maxreads,famout_sam,p->rlens_flat_out);
 
     fprintf(stderr, "\t-> Number of species with reads that map uniquely: %lu\n", specWeight.size());
 
diff --git a/ngsLCA_cli.cpp b/ngsLCA_cli.cpp
index bc24fe8..ac97f32 100644
--- a/ngsLCA_cli.cpp
+++ b/ngsLCA_cli.cpp
@@ -46,6 +46,7 @@ pars *pars_init() {
     p->ignore_errors = 0;
     p->useDump = 0;
     p->maxreads = -1;
+    p->rlens_flat_out = 0;
     return p;
 }
 
@@ -252,6 +253,8 @@ pars *get_pars(int argc, char **argv) {
 	    p->useDump = atoi(val);
 	else if (!strcasecmp("--ignore_errors", key)||!strcasecmp("-i", key))
 	    p->ignore_errors++;
+	else if (!strcasecmp("--rlens_flat_out", key)||!strcasecmp("-rfo", key))
+	  p->rlens_flat_out = atoi(val);
         else if (!strcasecmp("--temp", key)) {
             free(p->tempfolder);
             p->tempfolder = strdup(val);
@@ -295,7 +298,7 @@ pars *get_pars(int argc, char **argv) {
         fprintf(stderr, "\t-> Will output the reads that are used for damage file:\t\'%s\'\n", buf);
         p->usedreads_sam = strdup(buf);
     }
-    if (make_used_reads) {
+    if (make_famout_reads) {
       snprintf(buf, 1024, "%s.famoutreads.bam", p->outnames);
       fprintf(stderr, "\t-> Will output the reads that has lca below family :\t\'%s\'\n", buf);
       p->famout_sam = strdup(buf);
@@ -327,6 +330,7 @@ void print_pars(FILE *fp, pars *p) {
     fprintf(fp, "\t-> --ignore_errors\t%d\n", p->ignore_errors);
     fprintf(fp, "\t-> --temp\t%s\n", p->tempfolder);
     fprintf(fp, "\t-> --filtered_acc2tax\t%s\n", p->filteredAcc2taxfile);
+    fprintf(fp, "\t-> --rlens_flat_out\t%d\n", p->rlens_flat_out);
 }
 
 #ifdef __WITH_MAIN__
diff --git a/ngsLCA_cli.h b/ngsLCA_cli.h
index 6801a63..f8ffbdc 100644
--- a/ngsLCA_cli.h
+++ b/ngsLCA_cli.h
@@ -43,6 +43,7 @@ typedef struct {
     int ignore_errors;
   int useDump;
   long maxreads;
+  int rlens_flat_out;
 } pars;
 
 pars *get_pars(int argc, char **argv);
diff --git a/profile.cpp b/profile.cpp
index 675440e..300eb2e 100644
--- a/profile.cpp
+++ b/profile.cpp
@@ -13,11 +13,26 @@
 #include <cstring>    // for strlen, strtok, memset, strdup
 #include <vector>     // for vector
 
+
+void my_bgzf_write(BGZF *what, const void *the, size_t fuck ){
+  if(bgzf_write(what,the,fuck)!=(ssize_t) fuck){
+    fprintf(stderr,"\t-> Problem writing to file\n");
+    exit(1);
+  }
+}
+
+void my_bgzf_read(BGZF *what, void *the, size_t fuck ){
+  if(bgzf_read(what,the,fuck)!=(ssize_t) fuck){
+    fprintf(stderr,"\t-> Problem writing to file\n");
+    exit(1);
+  }
+}
+
 float **getmatrix(size_t x, size_t y) {
     float **ret = new float *[x];
     for (size_t i = 0; i < x; i++) {
         ret[i] = new float[y];
-        for (int j = 0; j < 16; j++)
+        for (size_t j = 0; j < y; j++)
             ret[i][j] = 0;
     }
     return ret;
@@ -174,8 +189,8 @@ void reconstructRefWithPosHTS(const bam1_t *b, std::pair<kstring_t *, std::vecto
 
     // combine the CIGAR and MD into one single string
     int mdVectorIndex = 0;
-
-    for (unsigned int i = 0; i < strlen(reconstructedTemp); i++) {
+    unsigned int funkydonky = at;
+    for (unsigned int i = 0; i < funkydonky; i++) {
         if (reconstructedTemp[i] == 'M') {  // only look at matches and indels
 
             if (mdVectorIndex < (int)parsedMD.size()) {  // still have mismatches
@@ -212,8 +227,8 @@ void reconstructRefWithPosHTS(const bam1_t *b, std::pair<kstring_t *, std::vecto
         }
     }
 
-    if (strlen(pp.first->s) != b->core.l_qseq) {
-        fprintf(stderr, "Could not recreate the sequence for read: %s pp.first->s: %s strlen():%lu\n", bam_get_qname(b), pp.first->s, strlen(pp.first->s));
+    if (strlen(pp.first->s) !=(size_t) b->core.l_qseq) {
+        fprintf(stderr, "Could not recreate the sequence for read: %s pp.first->s: %s strlen():%zu\n", bam_get_qname(b), pp.first->s, strlen(pp.first->s));
         exit(1);
     }
 
@@ -223,7 +238,7 @@ void reconstructRefWithPosHTS(const bam1_t *b, std::pair<kstring_t *, std::vecto
     }
 }
 
-inline void increaseCounters(const bam1_t *b, const char *reconstructedReference, const std::vector<int> &reconstructedReferencePos, const int &minQualBase, int MAXLENGTH, float **mm5p, float **mm3p, float incval) {
+inline void increaseCounters(const bam1_t *b, const char *reconstructedReference, const int &minQualBase, int MAXLENGTH, float **mm5p, float **mm3p, float incval) {
     const char *alphabetHTSLIB = "NACNGNNNTNNNNNNN";
     char refeBase;
     char readBase;
@@ -261,79 +276,102 @@ inline void increaseCounters(const bam1_t *b, const char *reconstructedReference
             refeBase = readBase;
         }
 
-        refeBase = refToChar[refeBase];
-        readBase = refToChar[readBase];
-
-        if (refeBase != 4 && readBase != 4) {
-            int dist5p = i;
-            int dist3p = b->core.l_qseq - 1 - i;
-
-            if (bam_is_rev(b)) {
-                refeBase = com[refeBase];
-                readBase = com[readBase];
-                // dist5p=int(al.QueryBases.size())-1-i;
-                dist5p = int(b->core.l_qseq) - 1 - i;
-                dist3p = i;
-            }
-
-            if (dist5p < MAXLENGTH)
-                mm5p[dist5p][toIndex[refeBase][readBase]] += incval;
-            if (dist3p < MAXLENGTH)
-                mm3p[dist3p][toIndex[refeBase][readBase]] += incval;
+        refeBase = refToChar[(unsigned char)refeBase];
+        readBase = refToChar[(unsigned char)readBase];
+	if(refeBase<4 && readBase <4){
+	  int dist5p = i;
+	  int dist3p = b->core.l_qseq - 1 - i;
+	  
+	  if (bam_is_rev(b)) {
+	    refeBase = com[(unsigned char)refeBase];
+	    readBase = com[(unsigned char)readBase];
+	    // dist5p=int(al.QueryBases.size())-1-i;
+	    dist5p = int(b->core.l_qseq) - 1 - i;
+	    dist3p = i;
+	  }
+	  
+	  if (dist5p < MAXLENGTH)
+	    mm5p[dist5p][toIndex[(unsigned char)refeBase][(unsigned char)readBase]] += incval;
+	  if (dist3p < MAXLENGTH)
+	    mm3p[dist3p][toIndex[(unsigned char)refeBase][(unsigned char)readBase]] += incval;
         }
     }
 }
 
+//this function will not get a commen
+void mysuperduper_function(triple &t, size_t old_m, size_t new_m) {
+  assert(t.rlens!=NULL&&old_m<new_m);
+  assert((size_t)t.rlens_m == old_m);//just so ser ser double sure
+  size_t *tmp = new size_t[new_m];
+  memset(tmp, 0, new_m * sizeof(size_t));
+  memcpy(tmp, t.rlens, t.rlens_m * sizeof(size_t));
+  delete [] t.rlens;
+  
+  t.rlens = tmp;
+  t.rlens_m = new_m;
+}
+
 int damage::damage_analysis(bam1_t *b, int which, float incval) {
+  size_t rlen = (size_t)b->core.l_qseq;
   //  fprintf(stderr,"\t-> incval: %f\n",incval);
+  //ok we need to do something
+  if (rlen + 10 >= temp_len) {//so lets be sure,
+    temp_len = b->core.l_qseq+10;//and then be really sure
+    kroundup32(temp_len);
+    free(reconstructedTemp);
+    reconstructedTemp = (char *)calloc(temp_len, 1);
+    
+  }
+  memset(reconstructedTemp, 0, temp_len);//maybe bottleneck, should not be needed.
+  
   if (assoc.find(which) == assoc.end()) {
-    triple val = {0, getmatrix(MAXLENGTH, 16), getmatrix(MAXLENGTH, 16),new size_t[500]};
-    for(int i=0;i<500;i++)
-      val.rlens[i] = 0;
+    triple val = {0, getmatrix(MAXLENGTH, 16), getmatrix(MAXLENGTH, 16),new size_t[temp_len],temp_len};
+    memset(val.rlens, 0, temp_len * sizeof(size_t));
     assoc[which] = val;
-    mm5pF = val.mm5pF;
-    mm3pF = val.mm3pF;
+    //mm5pF = val.mm5pF;
+    //mm3pF = val.mm3pF;
     //    fprintf(stderr,"has added which: %d\n",which);
   }
-    std::map<int, triple>::iterator it = assoc.find(which);
-    it->second.nreads++;
-    //    fprintf(stderr,"[%s] it->first:%d it->second.nreads:%d\n",__FUNCTION__,it->first,it->second.nreads);
-    assert(b->core.l_qseq<500);
-    it->second.rlens[b->core.l_qseq] =  it->second.rlens[b->core.l_qseq] + 1; 
-    if (b->core.l_qseq - 10 > temp_len) {
-        temp_len = b->core.l_qseq;
-        kroundup32(temp_len);
-        free(reconstructedTemp);
-        reconstructedTemp = (char *)calloc(temp_len, 1);
-    }
-    memset(reconstructedTemp, 0, temp_len);
-    reconstructRefWithPosHTS(b, reconstructedReference, reconstructedTemp);
-    increaseCounters(b, reconstructedReference.first->s, reconstructedReference.second, minQualBase, MAXLENGTH, it->second.mm5pF, it->second.mm3pF, incval);
-    return 0;
+
+  
+  std::map<int, triple>::iterator it = assoc.find(which);
+  //so lets make sure that 
+  if (rlen +10 >= it->second.rlens_m) {
+    size_t new_m = rlen + 10;
+    kroundup32(new_m);
+    mysuperduper_function(it->second, it->second.rlens_m, new_m);
+  }
+  
+  it->second.nreads++;
+  it->second.rlens[rlen]++;
+  //  fprintf(stderr,"er vi her fprintf%zu rlen: %zu\n",it->second.nreads,rlen);
+  reconstructRefWithPosHTS(b, reconstructedReference, reconstructedTemp);
+  increaseCounters(b, reconstructedReference.first->s, minQualBase, MAXLENGTH, it->second.mm5pF, it->second.mm3pF, incval);
+  return 0;
 }
 
 void damage::write(char *fname, bam_hdr_t *hdr) {
-    // fprintf(stderr,"Dumping asso.size(): %lu\n",assoc.size());
-    char *outname = strdup("metaout");
-    if (fname) {
-        free(outname);
-        outname = fname;
-    }
+    // fprintf(stderr,"Dumping asso.size(): %zu\n",assoc.size());
+  const char *outname = fname ? fname : "metaout";
     kstring_t kstr;
     kstr.l = kstr.m = 0;
     kstr.s = NULL;
     char onam[1024];
     snprintf(onam, 1024, "%s.res.gz", outname);
-    fprintf(stderr, "\t-> Will dump: \'%s\' this contains damage patterns for: %lu items\n", onam, assoc.size());
+    fprintf(stderr, "\t-> Will dump: \'%s\' this contains damage patterns for: %zu items\n", onam, assoc.size());
     BGZF *fp = my_bgzf_open(onam, nthreads);
 
     for (std::map<int, triple>::iterator it = assoc.begin(); it != assoc.end(); it++) {
         if (it->second.nreads == 0)  // should never happen
             continue;
-        if (hdr != NULL)
-	  ksprintf(&kstr, "%s\t%lu", hdr->target_name[it->first], it->second.nreads);
-        else
-            ksprintf(&kstr, "%lu", it->second.nreads);
+        if (hdr != NULL){
+	  if (it->first < 0 || it->first >= hdr->n_targets) {
+	    fprintf(stderr, "Invalid target index: %d\n", it->first);
+	    exit(1);
+	  }
+	  ksprintf(&kstr, "%s\t%zu", hdr->target_name[it->first], it->second.nreads);
+        }else
+            ksprintf(&kstr, "%zu", it->second.nreads);
         for (int l = 0; l < MAXLENGTH; l++) {
             for (int i = 0; i < 16; i++)
                 ksprintf(&kstr, "\t%.0f", it->second.mm5pF[l][i]);
@@ -343,58 +381,76 @@ void damage::write(char *fname, bam_hdr_t *hdr) {
                 ksprintf(&kstr, "\t%.0f", it->second.mm3pF[l][i]);
         }
         ksprintf(&kstr, "\n");
-        assert(bgzf_write(fp, kstr.s, kstr.l) == kstr.l);
+        my_bgzf_write(fp, kstr.s, kstr.l);
         kstr.l = 0;
     }
     bgzf_close(fp);
     free(kstr.s);
 }
 
-void damage::bwrite(char *fname) {
+//default output in rlens is ref bin:count[], FLAT_OUT 1 is ref tab bin tab count newline
+void damage::bwrite(char *fname,int FLAT_OUT = 0) {
 
     char onam[1024];
     snprintf(onam, 1024, "%s.bdamage.gz", fname);
-    fprintf(stderr, "\t-> Will dump: \'%s\' this contains damage patterns for: %lu items\n", onam, assoc.size());
+    fprintf(stderr, "\t-> Will dump: \'%s\' this contains damage patterns for: %zu items\n", onam, assoc.size());
     BGZF *fp = my_bgzf_open(onam, nthreads);
-    assert(bgzf_write(fp, &MAXLENGTH, sizeof(int)) == sizeof(int));
+    my_bgzf_write(fp, &MAXLENGTH, sizeof(int));
 
     for (std::map<int, triple>::iterator it = assoc.begin(); it != assoc.end(); it++) {
         if (it->second.nreads == 0)  // should never happen
             continue;
-        assert(bgzf_write(fp, &it->first, sizeof(int)) == sizeof(int));
-        assert(bgzf_write(fp, &it->second.nreads, sizeof(int)) == sizeof(int));
+        my_bgzf_write(fp, &it->first, sizeof(int));
+	int32_t nreads_out;
+	
+	if (it->second.nreads > std::numeric_limits<int32_t>::max()) {
+	  fprintf(stderr, "nreads overflow\n");
+	  exit(1);
+	}
+	
+	nreads_out = static_cast<int32_t>(it->second.nreads); 
+
+	
+        my_bgzf_write(fp, &nreads_out, sizeof(int));
         for (int l = 0; l < MAXLENGTH; l++)
-            assert(bgzf_write(fp, it->second.mm5pF[l], sizeof(int) * 16) == sizeof(int) * 16);
+            my_bgzf_write(fp, it->second.mm5pF[l], sizeof(float) * 16);
 
         for (int l = 0; l < MAXLENGTH; l++)
-            assert(bgzf_write(fp, it->second.mm3pF[l], sizeof(int) * 16) == sizeof(int) * 16);
+            my_bgzf_write(fp, it->second.mm3pF[l], sizeof(float) * 16);
     }
     bgzf_close(fp);
 
     snprintf(onam, 1024, "%s.rlens.gz", fname);
-    fprintf(stderr, "\t-> Will dump: \'%s\' this contains read length distributions for: %lu items\n", onam, assoc.size());
+    fprintf(stderr, "\t-> Will dump: \'%s\' this contains read length distributions for: %zu items\n", onam, assoc.size());
     fp = my_bgzf_open(onam, nthreads);
     kstring_t kstr3000;
     kstr3000.s = NULL;
     kstr3000.l = kstr3000.m = 0;
-    ksprintf(&kstr3000,"id");
-    for(int i=0;i<500;i++)
-      ksprintf(&kstr3000,"\trlen%d",i);
-    ksprintf(&kstr3000,"\n");
+    if(FLAT_OUT ==0 )
+      ksprintf(&kstr3000,"id\trlen:count\n");
+    else
+      ksprintf(&kstr3000,"ID\trlen\tcount\n");
     for (std::map<int, triple>::iterator it = assoc.begin(); it != assoc.end(); it++) {
         if (it->second.nreads == 0)  // should never happen
             continue;
-	ksprintf(&kstr3000,"%d",it->first);
-	for(int i=0;i<500-1;i++)
-	  ksprintf(&kstr3000,"\t%lu",it->second.rlens[i]);
-	ksprintf(&kstr3000,"\t%lu\n",it->second.rlens[499]);
+	if(FLAT_OUT==0)
+	  ksprintf(&kstr3000,"%d",it->first);
+	for(size_t i=0;i<it->second.rlens_m;i++)
+	  if(it->second.rlens[i]>0){
+	    if(FLAT_OUT == 0)
+	      ksprintf(&kstr3000,"\t%zu:%zu",i,it->second.rlens[i]);
+	    else
+	      ksprintf(&kstr3000,"%d\t%zu\t%zu\n",it->first,i,it->second.rlens[i]);
+	  }
+	if(FLAT_OUT==0)
+	  ksprintf(&kstr3000,"\n");
 	if(kstr3000.l>1000000){
-	  assert(bgzf_write(fp,kstr3000.s,kstr3000.l)==kstr3000.l);
+	  my_bgzf_write(fp,kstr3000.s,kstr3000.l);
 	  kstr3000.l  = 0;
 	}
 	  
     }
-    assert(bgzf_write(fp,kstr3000.s,kstr3000.l)==kstr3000.l);
+    my_bgzf_write(fp,kstr3000.s,kstr3000.l);
     if(kstr3000.l>0)
       free(kstr3000.s);
     kstr3000.l  = 0;
@@ -412,14 +468,17 @@ int printresults_grenaud2(FILE *fp, float **mm5p, int lengthMaxToPrint) {
     for (int l = 0; l < lengthMaxToPrint; l++) {
         fprintf(fp, "%*d\t", int(log10(lengthMaxToPrint)) + 1, l);
         for (int n1 = 0; n1 < 4; n1++) {
-            int totalObs = 0;
+            float totalObs = 0.0;
             for (int n2 = 0; n2 < 4; n2++)
                 totalObs += mm5p[l][4 * n1 + n2];
 
             for (int n2 = 0; n2 < 4; n2++) {
                 if (n1 == n2)
                     continue;
-                fprintf(fp, "%*.*f", 1 + 5 + 1, 5, std::max(0.0, double(mm5p[l][4 * n1 + n2]) / double(totalObs)));
+		if(totalObs>0)
+		  fprintf(fp, "%*.*f", 1 + 5 + 1, 5, std::max(0.0, double(mm5p[l][4 * n1 + n2]) / double(totalObs)));
+		else
+		  fprintf(fp,"0.0");
                 if (!(n1 == 3 && n2 == 2))
                     fprintf(fp, "\t");
             }
@@ -430,10 +489,12 @@ int printresults_grenaud2(FILE *fp, float **mm5p, int lengthMaxToPrint) {
 }
 
 void damage::printit(FILE *fp, int l) {
-    if (mm5pF)
-        printresults_grenaud2(fp, mm5pF, l);
-    if (mm3pF)
-        printresults_grenaud2(fp, mm3pF, l);
+  fprintf(fp,"Printing mismatchmatric for first entry\n");
+  if(assoc.size()==0)
+    return;
+  auto it = assoc.begin();
+  printresults_grenaud2(fp, it->second.mm5pF, l);
+  printresults_grenaud2(fp, it->second.mm3pF, l);
 }
 
 #ifdef __WITH_MAIN__
@@ -519,14 +580,14 @@ int main(int argc, char *argv[]) {
     sam_hdr_destroy(h);
     bam_destroy1(b);
     sam_close(fp);
-    fprintf(stderr, "nreads: %lu\n", dmg->assoc.begin()->second.nreads);
+    fprintf(stderr, "nreads: %zu\n", dmg->assoc.begin()->second.nreads);
     std::map<int,triple>::iterator it = dmg->assoc.begin();
     assert(it!=dmg->assoc.end());
     triple trpl = it->second;
     printresults_grenaud2(stdout, trpl.mm5pF, lengthMaxToPrint);
     printresults_grenaud2(stdout, trpl.mm3pF, lengthMaxToPrint);
     for (int i = 0; 0 & i < 16; i++)
-        fprintf(stdout, "%lu\t", dmg->mm5pF[0][i]);
+        fprintf(stdout, "%zu\t", dmg->mm5pF[0][i]);
     destroy_damage(dmg);
     return 0;
 }
@@ -542,28 +603,31 @@ std::map<int, double *> load_bdamage3(const char *fname, int howmany) {
 
     if (((bgfp = bgzf_open(infile, "r"))) == NULL) {
         fprintf(stderr, "Could not open input BDamage file: %s\n", infile);
-        exit(0);
+        exit(1);
     }
 
     std::map<int, double *> retmap;
 
     int printlength;
-    assert(sizeof(int) == bgzf_read(bgfp, &printlength, sizeof(int)));
+    my_bgzf_read(bgfp, &printlength, sizeof(int));
 
     if (howmany > printlength) {
         fprintf(stderr, "\t-> Problem binary file has data for: %d positions, but you are requesting merge with: %d positions \n", printlength, howmany);
         fprintf(stderr, "\t-> Solutions set -howmany to lower value\n");
-        exit(0);
+        exit(1);
     }
 
     int ref_nreads[2];
 
-    int data[16];
+
     while (1) {
         int nread = bgzf_read(bgfp, ref_nreads, 2 * sizeof(int));
         if (nread == 0)
             break;
-        assert(nread == 2 * sizeof(int));
+        if(nread!=2 * sizeof(int)){
+	  fprintf(stderr,"\t-> Tralala file looks corrupt\n");
+	  exit(1);
+	}
 
         double *formap = new double[1 + 3 * howmany];
         for (int i = 0; i < 1 + 3 * howmany; i++)
@@ -571,9 +635,9 @@ std::map<int, double *> load_bdamage3(const char *fname, int howmany) {
         //    fprintf(stderr,"formap: %p\n",formap);
         int incer = 0;
         formap[incer++] = ref_nreads[1];
-
+	float data[16];
         for (int at = 0; at < printlength; at++) {
-            assert(16 * sizeof(int) == bgzf_read(bgfp, data, sizeof(int) * 16));
+            my_bgzf_read(bgfp, data, sizeof(float) * 16);
             float flt[16];                 // this will contain the float representation of counts
             for (int i = 0; i < 4; i++) {  // loop over A*,C*,G*,T*
                 double tsum = 0;
@@ -597,7 +661,7 @@ std::map<int, double *> load_bdamage3(const char *fname, int howmany) {
             }
         }
         for (int at = 0; at < printlength; at++) {
-            assert(16 * sizeof(int) == bgzf_read(bgfp, data, sizeof(int) * 16));
+            my_bgzf_read(bgfp, data, sizeof(float) * 16);
 
             float flt[16];
             for (int i = 0; i < 4; i++) {
@@ -628,7 +692,7 @@ std::map<int, double *> load_bdamage3(const char *fname, int howmany) {
     //  exit(0);
     if (bgfp)
         bgzf_close(bgfp);
-    fprintf(stderr, "\t-> Done loading binary bdamage.gz file. It contains: %lu\n", retmap.size());
+    fprintf(stderr, "\t-> Done loading binary bdamage.gz file. It contains: %zu\n", retmap.size());
     for (std::map<int, double *>::iterator it = retmap.begin(); 0 && it != retmap.end(); it++)
         fprintf(stderr, "it->second:%p\n", it->second);
 
@@ -644,12 +708,12 @@ std::map<int, mydataD> load_bdamage_full(const char *fname, int &printlength) {
 
     if (((bgfp = bgzf_open(infile, "r"))) == NULL) {
         fprintf(stderr, "Could not open input BDamage file: %s\n", infile);
-        exit(0);
+        exit(1);
     }
 
     std::map<int, mydataD> retmap;
     printlength = 0;
-    assert(sizeof(int) == bgzf_read(bgfp, &printlength, sizeof(int)));
+    my_bgzf_read(bgfp, &printlength, sizeof(int));
 
     int ref_nreads[2];
 
@@ -657,7 +721,10 @@ std::map<int, mydataD> load_bdamage_full(const char *fname, int &printlength) {
         int nread = bgzf_read(bgfp, ref_nreads, 2 * sizeof(int));
         if (nread == 0)
             break;
-        assert(nread == 2 * sizeof(int));
+        if(nread!=2 * sizeof(int)){
+	  fprintf(stderr,"\t-> Tralala looks like a corrupt file\n");
+	  exit(1);
+	}
         mydataD md;
 	md.howmany = printlength;
         md.fwD = new double[16 * printlength];
@@ -666,13 +733,13 @@ std::map<int, mydataD> load_bdamage_full(const char *fname, int &printlength) {
 
         float tmp[16];
         for (int i = 0; i < printlength; i++) {
-            assert(16 * sizeof(float) == bgzf_read(bgfp, tmp, sizeof(float) * 16));
+            my_bgzf_read(bgfp, tmp, sizeof(float) * 16);
             for (int ii = 0; ii < 16; ii++)
                 md.fwD[i * 16 + ii] = tmp[ii];
         }
 
         for (int i = 0; i < printlength; i++) {
-            assert(16 * sizeof(float) == bgzf_read(bgfp, tmp, sizeof(float) * 16));
+            my_bgzf_read(bgfp, tmp, sizeof(float) * 16);
             for (int ii = 0; ii < 16; ii++)
                 md.bwD[i * 16 + ii] = tmp[ii];
         }
@@ -681,7 +748,7 @@ std::map<int, mydataD> load_bdamage_full(const char *fname, int &printlength) {
 
     if (bgfp)
         bgzf_close(bgfp);
-    fprintf(stderr, "\t-> Done loading binary bdamage.gz file. It contains: %lu\n", retmap.size());
+    fprintf(stderr, "\t-> Done loading binary bdamage.gz file. It contains: %zu\n", retmap.size());
 #if 0
   for(std::map<int,mydata>::iterator it = retmap.begin();it!=retmap.end();it++)
     fprintf(stderr,"it->second:%p\n",it->second);
@@ -720,7 +787,7 @@ std::map<int, mydata2> load_lcastat(const char *fname,int skipfirstline) {
     if (fp)
         gzclose(fp);
 
-    fprintf(stderr, "\t-> Done loading lcastat file It contains: %lu\n", retmap.size());
+    fprintf(stderr, "\t-> Done loading lcastat file It contains: %zu\n", retmap.size());
     for (std::map<int, mydata2>::iterator it = retmap.begin(); 0 && it != retmap.end(); it++)
         fprintf(stderr, "%d->(%d,%f,%f,%f,%f)\n", it->first, it->second.nreads, it->second.data[0], it->second.data[1], it->second.data[2], it->second.data[3]);
 
diff --git a/profile.h b/profile.h
index 712d85d..cdf2857 100644
--- a/profile.h
+++ b/profile.h
@@ -21,7 +21,8 @@ typedef struct {
   size_t nreads;
   float **mm5pF;
   float **mm3pF;
-  size_t *rlens;
+  size_t *rlens; //this is the counts. Not hashmap to speed up things
+  size_t rlens_m; // the length of the above. surprise
 } triple;
 
 class damage {
@@ -32,15 +33,14 @@ class damage {
     int MAXLENGTH;
     int minQualBase;  // currently not set; should be set in init_damage
     int nclass;
-    float **mm5pF;  // will point to first. Just a hack
-    float **mm3pF;  // will point to first. Just a hack
     std::pair<kstring_t *, std::vector<int> > reconstructedReference;
     std::map<int, triple> assoc;
     void write(char *prefix, bam_hdr_t *hdr);
-    void bwrite(char *prefix);
+  void bwrite(char *prefix,int FLAT_OUT);
     int damage_analysis(bam1_t *b, int whichclass, float incval);
     void printit(FILE *fp, int l);
-    int temp_len;
+  size_t temp_len;//<- this is the maxlenght of the reconstructed reference etc. This might be updated to a higher value.
+  //also when we initialize a new rlens for a ref/taxid that will be used. 
     damage(int maxlen, int nthd, int minqb) {
         temp_len = 512;
         MAXLENGTH = maxlen;
@@ -51,7 +51,6 @@ class damage {
         kstr->l = kstr->m = 0;
         kstr->s = NULL;
         reconstructedReference.first = kstr;
-        mm5pF = mm3pF = NULL;
     }
     ~damage() { free(reconstructedTemp); }
 };
diff --git a/test/output.md5 b/test/output.md5
index 1fe08d7..fc232a9 100644
--- a/test/output.md5
+++ b/test/output.md5
@@ -5,7 +5,7 @@ e46712d6bc1c9301abbc7467b5f720e6  output/test_aggregate.stat
 783c885b2ee8ac6469ab070b972e73d8  output/test_getdamage_global.bdamage
 3398e353af15d1646228cb8bbaf4a5a8  output/test_getdamage_global.bdamage.tsv
 f0882ca11587cea63aa120d1a64c8d14  output/test_getdamage_global.res
-3fda91168f2a07d82312e4589841df1e  output/test_getdamage_global.stat
+2052a9eb4a55313dafbf0ef80a6b751d  output/test_getdamage_global.stat
 71308047bbc8f703fdd65c798de91cfa  output/test_getdamage_local.bdamage
 93ec5114d069f984ccaaa5079c15db1e  output/test_getdamage_local.bdamage.tsv
 d95e54a3886b9f411bde8d96de4af32e  output/test_getdamage_local.res
@@ -22,7 +22,7 @@ dc5bf1edd27b494e23f72d41d718600d  data/f570b1db7c.dedup.filtered.bam
 #c73c56f8a3e2dc9b13c949ec1404bd50  data/f570b1db7c.dedup.filtered.rname.bam
 #3568c671c0fa09797284746bb0528ddd  data/f570b1db7c.dedup.filtered.rname.bam ##apparantly samtools 1.20 has changed the within qname order when sorting by readname. This is of cause fine but some seqs would be reverse complemented. When switching to the new order then this entry should be used instead of the above
 2ae4149efce77b91786d4a13b5451c3b  output/test_lca.lca.sub
-53943f1ceab5c86c6d815fc55d0363c7  output/test_getdamage_global.rlens
-57013f84420d6460ae5e15aa0092dc1e  output/test_getdamage_local.rlens
-abc7a9ec248d48204d5f417f4c30c59d  output/test_lca.rlens
-9574c8051682dc50c040e5e8c16704a6  output/test_mergedamage.rlens
+ce19e1efdedb24f0353dc52fab2c985b  output/test_getdamage_global.rlens
+8700fac9ccb1414c9004cfcc4a122829  output/test_getdamage_local.rlens
+16e2ed7dbf1959748f6a41d5aaf35d05  output/test_lca.rlens
+1d038fdfdb159ab5efb069648c14b775  output/test_mergedamage.rlens
\ No newline at end of file
diff --git a/test/output.md5.macos b/test/output.md5.macos
index c48b285..949f696 100644
--- a/test/output.md5.macos
+++ b/test/output.md5.macos
@@ -3,7 +3,7 @@ e46712d6bc1c9301abbc7467b5f720e6  output/test_aggregate.stat
 783c885b2ee8ac6469ab070b972e73d8  output/test_getdamage_global.bdamage
 3398e353af15d1646228cb8bbaf4a5a8  output/test_getdamage_global.bdamage.tsv
 f0882ca11587cea63aa120d1a64c8d14  output/test_getdamage_global.res
-3fda91168f2a07d82312e4589841df1e  output/test_getdamage_global.stat
+2052a9eb4a55313dafbf0ef80a6b751d  output/test_getdamage_global.stat
 71308047bbc8f703fdd65c798de91cfa  output/test_getdamage_local.bdamage
 93ec5114d069f984ccaaa5079c15db1e  output/test_getdamage_local.bdamage.tsv
 d95e54a3886b9f411bde8d96de4af32e  output/test_getdamage_local.res
@@ -19,7 +19,7 @@ eaa8bb5f32aa0ef05bcf923d0c2bd711  output/test_lca_taxa.uglyprint.mismatch
 a23c34eaca4a1d1742dd65f98d70efbf  output/test_mergedamage.bdamage
 dc5bf1edd27b494e23f72d41d718600d  data/f570b1db7c.dedup.filtered.bam
 2ae4149efce77b91786d4a13b5451c3b  output/test_lca.lca.sub
-53943f1ceab5c86c6d815fc55d0363c7  output/test_getdamage_global.rlens
-57013f84420d6460ae5e15aa0092dc1e  output/test_getdamage_local.rlens
-abc7a9ec248d48204d5f417f4c30c59d  output/test_lca.rlens
-9574c8051682dc50c040e5e8c16704a6  output/test_mergedamage.rlens
\ No newline at end of file
+ce19e1efdedb24f0353dc52fab2c985b  output/test_getdamage_global.rlens
+8700fac9ccb1414c9004cfcc4a122829  output/test_getdamage_local.rlens
+16e2ed7dbf1959748f6a41d5aaf35d05  output/test_lca.rlens
+1d038fdfdb159ab5efb069648c14b775  output/test_mergedamage.rlens