Index: /ICOSIM/CPP/trunk/source/Collimator.h
===================================================================
--- /ICOSIM/CPP/trunk/source/Collimator.h	(revision 16)
+++ /ICOSIM/CPP/trunk/source/Collimator.h	(revision 17)
@@ -52,5 +52,5 @@
 
 
-    virtual void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath) {};
+    virtual void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, long int& nhit, string outputpath) {};
     //these methods are specific for a fluka/crystal collimator. They are just defined here.
 #if defined(FLUKA)
Index: /ICOSIM/CPP/trunk/source/CrystalCollimator.cc
===================================================================
--- /ICOSIM/CPP/trunk/source/CrystalCollimator.cc	(revision 16)
+++ /ICOSIM/CPP/trunk/source/CrystalCollimator.cc	(revision 17)
@@ -32,7 +32,7 @@
 
 
-void CrystalCollimator::collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath)
+void CrystalCollimator::collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, long int& nhit, string outputpath)
 {
-
+//pas: number of turns
     //Calculate proton energy
     double Mproton(0.93827231); //rest mass [GeV]
@@ -62,6 +62,8 @@
 
         for (int k(0); k < bunch.size(); ++k) {
-
+            //simulate the particles that are not lost on the previous aperture
+	    if(bunch[k].inabs !=0){   
             sortie  << bunch[k].coordonnees[0][0] << "," << bunch[k].coordonnees[0][2] << "," << bunch[k].coordonnees[0][1] << "," << bunch[k].coordonnees[0][3] << "," << Eproton[k] << endl;
+	    }
         }
 
@@ -93,6 +95,13 @@
    // sim.general(sim, pas, outputpath, NAME, emitx0, emity0,enum, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);
     
-    sim.general(sim, pas, outputpath, Mirror, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);
+    sim.general(sim, pas, nhit, outputpath, Mirror, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);
 
+    
+    //assign the new bunch coordinates after the crystal
+    for (int k(0); k < bunch.size(); ++k) {
+      bunch[k].inabs =0 ;
+    }
+    
+    
     ifstream enter;
     string file_in;
@@ -105,5 +114,5 @@
     } else {
 
-        for (int k(0); k < bunch.size(); ++k) {
+        for (int k(0); k < bunch.size()-nhit; ++k) {
 
             string rest;
@@ -120,7 +129,11 @@
             getline(enter, rest, ',');
             bunch[k].coordonnees[1][3] = atof(rest.c_str());
-
+            
+	    //particle is not lost (interact) with the crystal
+	    bunch[k].inabs = 1; 
+	    
             getline(enter, rest);
             Eproton[k] = atof(rest.c_str());
+	 
         }
 
@@ -128,4 +141,5 @@
     }
 
+    
     vector <double> gamProton;
     vector <double> betProton;
Index: /ICOSIM/CPP/trunk/source/CrystalCollimator.h
===================================================================
--- /ICOSIM/CPP/trunk/source/CrystalCollimator.h	(revision 16)
+++ /ICOSIM/CPP/trunk/source/CrystalCollimator.h	(revision 17)
@@ -44,5 +44,5 @@
     //method specific to the passage through a Crystal collimator
 
-    void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath);
+    void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, long int& nhit, string outputpath);
 
 
Index: /ICOSIM/CPP/trunk/source/PlotRunSummary1.p
===================================================================
--- /ICOSIM/CPP/trunk/source/PlotRunSummary1.p	(revision 16)
+++ /ICOSIM/CPP/trunk/source/PlotRunSummary1.p	(revision 17)
@@ -1,17 +1,11 @@
 #Gnuplot script file for the PlotRunSummary (Plot of the loss map)
 
+set autoscale xmax
+set autoscale y
+set key default
+set xrange [-1:]
+set ylabel "P' (W/m)"
+set xlabel "Distance from starting point (km)"
+set title "Loss map"
 
-   set autoscale xmax
-   set autoscale y
-   set key default
-   set xrange [-1:]
-   set ylabel "P'(W/m)"
-   set xlabel "Distance from starting point (km)"
-   set title "Loss map"
-
-   load 'labelsIP.p'
-   
-   plot "PlotRunSummary1.dat" using 2:1 title 'Losses' with boxes, \
-      "PlotRunSummary2.dat" using 2:1 title 'IP' with impulses
-
-	 load 'labelsIPend.p'
+plot "PlotRunSummary1.dat" using 2:1 title 'Losses' with boxes
Index: /ICOSIM/CPP/trunk/source/PlotRunSummary2.p
===================================================================
--- /ICOSIM/CPP/trunk/source/PlotRunSummary2.p	(revision 16)
+++ /ICOSIM/CPP/trunk/source/PlotRunSummary2.p	(revision 17)
@@ -1,13 +1,14 @@
-#Gnuplot script for the PlotRunSummary (Plot of the collimator load)
+#Gnuplot script file for the PlotRunSummary (Plot of the collimator load)
 
+set autoscale
+set key default
+set style fill solid 0.5
+set ylabel "P_coll (W)"
+set xlabel "Collimators"
+set title "Collimator load distribution"
 
-   set autoscale
-   set key default
-    set title "Collimator load distribution"
-    set ylabel "P_coll (W)"
-    set xlabel "Collimators"
-    set style fill solid 0.5
-    load 'labelscolli.p'
-    plot "PlotRunSummary3.dat" using 1:2 title 'Collimator load' with boxes
+load 'labelscolli.p'
 
-       load 'labelscolliend.p'
+ plot "PlotRunSummary3.dat" using 1:2 title 'Collimator load' with boxes
+
+load 'labelscolliend.p'
Index: /ICOSIM/CPP/trunk/source/PlotRunSummary3.p
===================================================================
--- /ICOSIM/CPP/trunk/source/PlotRunSummary3.p	(revision 16)
+++ /ICOSIM/CPP/trunk/source/PlotRunSummary3.p	(revision 17)
@@ -1,16 +1,13 @@
-#Gnuplot script for PlotRunSummary (Plot of the development after the first impact on collimator)
+#Gnuplot script file for the PlotRunSummary (Plot of the  development after the first impact on collimator)
 
+set autoscale
+set key default
+set ylabel "Number of particles"
+set xlabel "N_turn"
+set title "Time Development after first impact on collimator"
 
-   set autoscale
-   set key default
-    set title "Time Development after first impact on collimator"
-    set xlabel "N_turn"
-    set ylabel "Number of particles"
-    set logscale y
+set logscale y
 
-    plot "PlotRunSummary4.dat" using 1:2 title 'particles in beam' with line, \
-	       "PlotRunSummary4.dat" using 1:3 title 'particles lost on collimators' with line, \
-	       "PlotRunSummary4.dat" using 1:4 title 'particles lost elsewhere' with line
+ plot "PlotRunSummary4.dat" using 1:2 title 'rest particles in beam' with line, "PlotRunSummary4.dat" using 1:3 title 'particles lost on collimators' with line, "PlotRunSummary4.dat" using 1:4 title 'particles lost elsewhere' with line
 
-
-		  unset logscale y
+unset logscale y
Index: /ICOSIM/CPP/trunk/source/PlotTrajectory1.p
===================================================================
--- /ICOSIM/CPP/trunk/source/PlotTrajectory1.p	(revision 16)
+++ /ICOSIM/CPP/trunk/source/PlotTrajectory1.p	(revision 17)
@@ -1,18 +1,15 @@
-#Gnuplot script for PlotTrajectory plot in x
+#Gnuplot script file for the PlotTrajectory in x
 
+set autoscale x
+set yrange [-25:25]
+ set style line 1 lw 1
+set ylabel "x [mm]"
+set xlabel "S-S_IP1 [km]"
+set title "Trajectory in x"
 
-   set autoscale x
-   set yrange [-25:25]
-   set xlabel "S-S_IP1 [km]"
-   set ylabel "x [mm]"
-   set title "Trajectory in x"
-   set style line 1 lw 1
-   set style line 2 lw 1 lt 6
-   set key outside below
+set style line 2 lw 1 lt 6
 
+set key outside below
 
-   load 'plotxco.p'
-
-
-   unset key
-		     
+load 'plotxco.p'
+unset key
Index: /ICOSIM/CPP/trunk/source/PlotTrajectory2.p
===================================================================
--- /ICOSIM/CPP/trunk/source/PlotTrajectory2.p	(revision 16)
+++ /ICOSIM/CPP/trunk/source/PlotTrajectory2.p	(revision 17)
@@ -1,17 +1,15 @@
-#Gnuplot script for PlotTrajectory plot in y
+#Gnuplot script file for the PlotTrajectory in y
 
+set autoscale x
+set yrange [-25:25]
+ set style line 1 lw 1
+set ylabel "y [mm]"
+set xlabel "S-S_IP1 [km]"
+set title "Trajectory in x"
 
-   set autoscale x
-   set yrange [-25:25]
-   set xlabel "S-S_IP1 [km]"
-   set ylabel "y [mm]"
-   set title "Trajectory in y"
-   set style line 1 lw 1
-   set style line 2 lw 1 lt 6
-   set key outside below
+set style line 2 lw 1 lt 6
 
+set key outside below
 
-   load 'plotyco.p'
-
-
-   unset key
+load 'plotyco.p'
+unset key
Index: /ICOSIM/CPP/trunk/source/crystal_dan_FINAL_VERSION_CHvsVR.f
===================================================================
--- /ICOSIM/CPP/trunk/source/crystal_dan_FINAL_VERSION_CHvsVR.f	(revision 16)
+++ /ICOSIM/CPP/trunk/source/crystal_dan_FINAL_VERSION_CHvsVR.f	(revision 17)
@@ -453,6 +453,6 @@
           
 c          write(*,*) "y = ", y
-           write(*,*) "before enter VR: x = ", x, "  xp = ",xp
-           write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
+c           write(*,*) "before enter VR: x = ", x, "  xp = ",xp
+c           write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
             
           xp=xp+0.45*(xp/xpcrit+1)*Ang_avr
@@ -465,6 +465,6 @@
      +      xp ,yp,PC)
      
-          write(*,*) "after enter VR: x = ", x, "  xp = ",xp
-          write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
+c          write(*,*) "after enter VR: x = ", x, "  xp = ",xp
+c          write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
           
           x = x + 0.5*s_length * xp
@@ -501,21 +501,21 @@
 !     +          DLRi(IS),xp ,yp,PC)
 
-            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
-            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
+c            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
+c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
            
             CALL CALC_ION_LOSS(IS,PC,s_length-Srefl,DESt)
             
-            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
-            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
-            write(*,*) "IS = ", IS, " NAM = ", NAM
-            write(*,*)  "s_length-Srefl = ",s_length-Srefl
-            write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS)
-            write(*,*) "DLRi(IS) = ", DLRi(IS)
+c            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
+c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
+c            write(*,*) "IS = ", IS, " NAM = ", NAM
+c            write(*,*)  "s_length-Srefl = ",s_length-Srefl
+c            write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS)
+c            write(*,*) "DLRi(IS) = ", DLRi(IS)
             
             CALL MOVE_AM_(IS,NAM,s_length-Srefl,DESt,DLYi(IS),
      +          DLRi(IS),xp ,yp,PC)
      
-            write(*,*) "after enter VR: x = ", x, "  xp = ",xp
-            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
+c            write(*,*) "after enter VR: x = ", x, "  xp = ",xp
+c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
            
             x = x + 0.5 * xp * (s_length - Srefl)
@@ -618,5 +618,5 @@
 c      endif
 c      
-      WRITE(*,*)'Crystal process: ',PROC
+c      WRITE(*,*)'Crystal process: ',PROC
 c      write(*,*)'xp_final :', xp
 c      WRITE(*,*)'Crystal process: ',PROC,'Chann Angle',Ch_angle/1000,   +
@@ -838,10 +838,10 @@
 
      
-      write(*,*) "In Move_AM(): "
+c      write(*,*) "In Move_AM(): "
 
 
       xp=xp*1000
       yp=yp*1000
-      write(*,*)'xp initial:', xp, 'yp initial', yp
+c      write(*,*)'xp initial:', xp, 'yp initial', yp
 C. DEI - dE/dx stoping energy
       PC    = PC - DEI*DZ    ! energy lost because of ionization process[GeV]
@@ -849,11 +849,11 @@
 C. DYA - rms of coloumb scattering
       DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
-      write(*,*)'dya= ',dya
+c      write(*,*)'dya= ',dya
       
       ran_x = RAN_GAUSS(1.)
-      write(*,*) "ran_x = ", ran_x
+c      write(*,*) "ran_x = ", ran_x
       
       ran_y = RAN_GAUSS(1.)
-      write(*,*) "ran_y = ", ran_y
+c      write(*,*) "ran_y = ", ran_y
       
       kxmcs = DYA*ran_x
@@ -867,6 +867,6 @@
       
 c       write(*,*) "ran_x = ", ran_x, "ran_y = ", ran_y
-      write(*,*) "kxmcs = ", kxmcs, "kymcs = ", kymcs
-      write(*,*)'xp + kxmcs:', XP, 'yp + kxmcs', YP
+c      write(*,*) "kxmcs = ", kxmcs, "kymcs = ", kymcs
+c      write(*,*)'xp + kxmcs:', XP, 'yp + kxmcs', YP
       
 cc     XP    = XP+DYA*RAN_GAUSS(1.)
@@ -938,6 +938,6 @@
       endif
 
-      write(*,*) "ichoix = ", ichoix
-       write(*,*)'xp final:', xp, 'yp final', yp
+c      write(*,*) "ichoix = ", ichoix
+c       write(*,*)'xp final:', xp, 'yp final', yp
        
 !---------- calculate the related kick -----------
@@ -980,5 +980,5 @@
 
 
-      write(*,*) "tx = ", tx, " tz = ", tz
+c      write(*,*) "tx = ", tx, " tz = ", tz
       
       XP = XP + tx
@@ -990,5 +990,5 @@
       yp=yp/1000
 
-      write(*,*)'xp final:', xp, 'yp final', yp
+c      write(*,*)'xp final:', xp, 'yp final', yp
 !-----------------------------------------------
 ! print out the energy loss and process experienced
Index: /ICOSIM/CPP/trunk/source/crystalprotonfuns.f
===================================================================
--- /ICOSIM/CPP/trunk/source/crystalprotonfuns.f	(revision 16)
+++ /ICOSIM/CPP/trunk/source/crystalprotonfuns.f	(revision 17)
@@ -28,5 +28,5 @@
       save
  
-      write(*,*) "In RAN_GAUSS: "
+c      write(*,*) "In RAN_GAUSS: "
 c      write(*,*) "flag = ",flag
 c      write(*,*) "U1 = ", U1,"U2 = ", U2
@@ -53,5 +53,5 @@
           RAN_GAUSS = x
 
-          write(*,*) "ran_gauss = ", RAN_GAUSS
+c          write(*,*) "ran_gauss = ", RAN_GAUSS
 
         return
Index: /ICOSIM/CPP/trunk/source/lattice.cc
===================================================================
--- /ICOSIM/CPP/trunk/source/lattice.cc	(revision 16)
+++ /ICOSIM/CPP/trunk/source/lattice.cc	(revision 17)
@@ -441,4 +441,7 @@
 {
 
+  cout << "bunch size is: " << bunch.size() << endl;
+  
+  
     if (idpart >= 0) {
         if (idpart > bunch.size()) {
@@ -602,9 +605,11 @@
           read(bunch);
           }*/
-
+	
+	
         for (int p(0); p < bunch.size(); ++p) { //loop througt the particles
 
             double dpopeff;
 
+	    // .inabs = 0, particle lost, .inabs = 1, particle not lost 
             if (bunch[p].inabs == 1) {
 
@@ -743,9 +748,14 @@
 
                             int lost(bunch.size());
-
-                            resColli[icop]->collipassCrystal(bunch, betgam, turn, outputpath);
-
-                            nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
-
+                             
+			    long int nhit = 0;
+			    
+                            resColli[icop]->collipassCrystal(bunch, betgam, turn, nhit, outputpath);
+
+                            //nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
+                            
+                             nhitcolli[icop] = nhit; 
+
+                         //    cout << "after the crystal: bunch size = " << bunch.size()<< endl;
                         }
 
@@ -924,4 +934,8 @@
                         if (inside == false) {
                             cout <<"The particle is lost at element " << reseau[i]->NAME << endl;
+// 			     cout << "Particle ID is: " << p << endl;
+// 			     cout << "bunch size is: " << bunch.size() << endl;
+// 			     cout << "bunch status is: " << bunch[p].inabs << endl;
+// 			     
                             xh0 = bunch[p].coordonnees[0][0];
                             bunch[p].inabs = 0;
@@ -1100,4 +1114,5 @@
             }
         }
+    //    cout << "bunch size is ...." << bunch.size() << endl;
 
         if (plotflag == "Yes") {
@@ -1132,5 +1147,5 @@
     vector <Particle> tempinabs;
 
-    //we contine just with particles that are not lost
+    //we contine to track for the next turn; just with particles that are not lost 
     for (int w(0); w < bunch.size(); ++w) {
         if (bunch[w].inabs != 0) {
@@ -1139,4 +1154,6 @@
     }
 
+ //   cout << "bunch size before the turn: "<< bunch.size() << endl; 
+    
     bunch.clear();
 
@@ -1147,3 +1164,4 @@
     tempinabs.clear();
 
-}
+   // cout << "bunch size after the turn " <<bunch.size() << endl;
+}
Index: /ICOSIM/CPP/trunk/source/partcrys.h
===================================================================
--- /ICOSIM/CPP/trunk/source/partcrys.h	(revision 16)
+++ /ICOSIM/CPP/trunk/source/partcrys.h	(revision 17)
@@ -44,6 +44,9 @@
     int nvcapt;
     int nvrefl;
+    int nabsorbed;
     int nout;          //number of particles that not hit the crystal collimator
-    int effet;        //Number that defini the effect :1 --> Out ;;; 2 --> Amorphous ;;; 3 --> Channeling ;;; 4 --> Dechanneling ;;; 5 --> V. Reflection ;;; 6 --> V. Capture
+    int effet;        //Number that defini the effect :1 --> Out ;;; 2 --> Amorphous ;;; 3 --> Channeling ;;; 4 --> Dechanneling ;;; 5 --> V. Reflection ;;; 6 --> V. Capture, 7 --> absorbed
+   long int nhit; //number of particles that hit the crystal collimator 
+    
 
 };
Index: /ICOSIM/CPP/trunk/source/simcrys.cc
===================================================================
--- /ICOSIM/CPP/trunk/source/simcrys.cc	(revision 16)
+++ /ICOSIM/CPP/trunk/source/simcrys.cc	(revision 17)
@@ -576,5 +576,5 @@
 * original crystal routine
 * *******************************************************************************/
-void SimCrys::general(SimCrys& sim, const int& pas, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt)
+void SimCrys::general(SimCrys& sim, const int& pas, long int& nhit, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt)
 {
     srand ( time(NULL) );
@@ -659,5 +659,7 @@
 
 
-
+    
+     
+     
     //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS-------------------------------------------------------------------------
     c_length = crys.Cry_length;
@@ -692,5 +694,17 @@
        * xmin=x-, xmax=x+, etc. from the input file
        */
-     long int    nhit    = 0;
+   
+
+      //initialize particle variables
+      part.n_part = 0;
+      part.nabsorbed = 0;
+      part.namor = 0;
+      part.nchann = 0;
+      part.ndechann = 0;
+      part.nout =0;
+      part.nvcapt = 0;
+      part.nvrefl = 0; 
+      part.nhit = 0;
+      nhit    = 0;
     double fracab  = 0;
   //  int   n_chan  = 0;          // valentina :initialize to zero the counters for crystal effects 
@@ -748,4 +762,8 @@
       while (!entree.eof()) {
 	count = count + 1;
+	
+	//particle hit the crystal or not
+	bool hitflag = false;
+	  
 	
 	getline(entree, rest, ',');
@@ -927,5 +945,5 @@
 	      s_impact = s_in0;       //!(for the first impact)
               
-	      if(1){
+	      if(0){
 		cout<<"hit the cry entrance face"<<endl;
 		cout<<"impact at s = "<< s_impact<< ",  x = "<<x_in0<<endl;
@@ -1062,5 +1080,9 @@
 	      //cout<<"----------------------1111111111111111111111111-----------------"<<endl;
 	      
-	
+	   if(layerflag == 0){
+	      part.namor = part.namor + 1;     /////////////COUNT/////////////
+              part.effet = 2;
+	      }
+	      
 	      if (strncmp(Proc,"out",3)==0) 
 		crossing  = false;
@@ -1069,5 +1091,7 @@
 	      
 	      if (crossing == true){
-		nhit = nhit + 1;
+ 		hitflag = true;
+ 		nhit = nhit + 1;
+ 		part.nhit = nhit; 
 	//	lhit = 100000000*ie + ITURN;
 		impact = x_in0;            
@@ -1261,6 +1285,7 @@
 		      xp_in0 = xp_shift;
 		    }                 ////////////////////////// fin IF NUMERO 17//////////////////
-		    
+		    hitflag = true;
 		    nhit = nhit + 1;
+		    part.nhit = nhit; 
 		    //lhit = 100000000*ie + ITURN;
 		    impact = x_in0;            
@@ -1339,5 +1364,5 @@
 
 
-	    nabs = 0; 
+	    nabs = 0;   //particle lost or not lost.....
 	    part.effet = nabs;
 	    //cout<< "debug - S Coll RF 2" ,  s_rot <<endl;
@@ -1375,4 +1400,19 @@
 	     part.nvcapt = part.nvcapt + 1;     /////////////COUNT/////////////
             part.effet = 6;}
+             if(strncmp(Proc,"absorbed",8==0)){
+ 	     part.nabsorbed = part.nabsorbed +1;
+ 	      part.effet = 7;
+ 	    }
+	     /*           
+ if (PROC(1:3).eq.'pne')
+                 PROC_dan=7
+             if (PROC(1:3).eq.'ppe')
+                 PROC_dan=8
+               if (PROC(1:4).eq.'diff')
+                 PROC_dan=9
+               if (PROC(1:4).eq.'ruth')
+                 PROC_dan=10
+	     */
+    
             
               //  ?????????????????????????????????????
@@ -1471,11 +1511,15 @@
 	  av_pc1 = av_pc1 + part.PC;
 
-
+      if(0){
 	  cout<<"  X apres :"<<part.x<<"  XP apres :"<<part.xp<<"  Y apres :"<<part.y<<"  YP apres :"<<part.yp<<"  PC apres :"<<part.PC<<endl;
 	  //cout<<endl;
+      }
+      //only record the particles that hit the collimator.
+      if(hitflag != true){
 	  //sortieIco.setf(ios::scientific);
-	  sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl;
+	  sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC  << endl;
 	  //}//end of for....
-	  
+      }
+      
 	}
       } //End of while
@@ -1525,4 +1569,9 @@
     file_out(pas, outputpath);
     
+    if(1){
+      cout << "number of particles before the crystal: " << count<< endl;
+      cout << "number of particle hit the crystal: " << part.nhit << endl;
+      cout << "number of particles after the crystal: " << count - part.nhit << endl;
+    }
     //cout<<endl;
     //cout<<endl;
@@ -1667,5 +1716,6 @@
     cout << "N.OUT :" << part.nout << endl;
     cout << "N.VREFLECTION :" << part.nvrefl << endl;
-    cout << "N.VCAPTURE :" << part.nvcapt << endl;
+    cout << "N.VCAPTURE:" << part.nvcapt << endl;
+    cout << "N.ABSORBED: "<< part.nabsorbed << endl;
 
     //cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!";
@@ -1688,5 +1738,5 @@
         out << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
         out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
-        out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
+        out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << 1e6*crys.Cry_length / crys.Rcurv << endl;
         out << "Critical Radius :" << crys.Rcrit << endl;
         out << "Ratio :" << crys.ratio << endl;
@@ -1702,6 +1752,10 @@
         out << "N.OUT :" << part.nout << endl;
         out << "N.VREFLECTION :" << part.nvrefl << endl;
-        out << "N.VCAPTURE :" << part.nvcapt << endl << endl << endl;
-
+        out << "N.VCAPTURE :" << part.nvcapt << endl ;
+	out << "N.ABSORBED: "<< part.nabsorbed << endl;
+	out << "N. of particles that hit the crystal collimator: " << part.nhit << endl;
+	cout << endl;
+	cout << endl;
+	
         out.close();
     }
Index: /ICOSIM/CPP/trunk/source/simcrys.h
===================================================================
--- /ICOSIM/CPP/trunk/source/simcrys.h	(revision 16)
+++ /ICOSIM/CPP/trunk/source/simcrys.h	(revision 17)
@@ -53,5 +53,5 @@
     //THE GENERAL FUNCTION OF THE TEST THAT REGROUP ALL THE PREVIOUS FUNCTIONS!!!!
 
-    void general(SimCrys& sim, const int& pas, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt);
+    void general(SimCrys& sim, const int& pas, long int& nhit, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt);
 
 
Index: /ICOSIM/CPP/trunk/source/simulation.cc
===================================================================
--- /ICOSIM/CPP/trunk/source/simulation.cc	(revision 16)
+++ /ICOSIM/CPP/trunk/source/simulation.cc	(revision 17)
@@ -78,4 +78,6 @@
     particlefile = inputfile + "initial.dat";
 
+    cout << "READ THE BUNCH OF PARTICLES FROM THE FILE: " << particlefile <<endl;
+    
     //if the name of the file is not 'initial.dat', uncomment the following two lines
 
@@ -663,5 +665,5 @@
       cout <<"********************************"<<endl;
         cout << "RunningFlag = " <<  RunningFlag << endl;
-	cout << "READ THE BUNCH OF PARTICLES FROM THE FILE: " <<inputfile <<endl;
+	
         readParticle(Apr, Zpr, inputfile);
 
@@ -676,5 +678,5 @@
       cout << "RunningFlag = " <<  RunningFlag << endl;
       cout << "GENERATION OF THE BUNCH OF PARTICLES." << endl;
-        genbunch(i0, nparti, partdistr, r1r2skin, nsigi * nsigi * emix, nsigi * nsigi * emiy, sigdpp, lat.reseau[i0]->BETX, lat.reseau[i0]->ALFX, lat.reseau[i0]->DX, lat.reseau[i0]->DPX, lat.reseau[i0]->BETY, lat.reseau[i0]->ALFY, lat.reseau[i0]->DY, lat.reseau[i0]->DPY, nsigi);
+        genbunch(i0, nparti, partdistr, r1r2skin, nsigi * nsigi * emix, nsigi * nsigi * emiy, sigdpp, lat.reseau[i0]->BETX, lat.reseau[i0]->ALFX, lat.reseau[i0]->DX, lat.reseau[i0]->DPX, lat.reseau[i0]->BETY, lat.reseau[i0]->ALFY, lat.reseau[i0]->DY, lat.reseau[i0]->DPY, nsigi, inputfile);
 
         //we set the id of the particles
@@ -742,5 +744,4 @@
       cout <<"********************************"<<endl;
 	cout << "RunningFlag = " <<  RunningFlag << endl;
-	cout << "READ THE BUNCH OF PARTICLES FROM THE FILE: " <<inputfile <<endl;
         readParticle(Apr, Zpr, inputfile);
     } else if (RunningFlag == 3) {
@@ -749,5 +750,5 @@
       cout << "RunningFlag = " <<  RunningFlag << endl;
       cout << "GENERATION OF THE BUNCH OF PARTICLES." << endl;
-        genbunch(i0, nparti, partdistr, r1r2skin, nsigi * nsigi * emix, nsigi * nsigi * emiy, sigdpp, lat.reseau[i0]->BETX, lat.reseau[i0]->ALFX, lat.reseau[i0]->DX, lat.reseau[i0]->DPX, lat.reseau[i0]->BETY, lat.reseau[i0]->ALFY, lat.reseau[i0]->DY, lat.reseau[i0]->DPY, nsigi);
+        genbunch(i0, nparti, partdistr, r1r2skin, nsigi * nsigi * emix, nsigi * nsigi * emiy, sigdpp, lat.reseau[i0]->BETX, lat.reseau[i0]->ALFX, lat.reseau[i0]->DX, lat.reseau[i0]->DPX, lat.reseau[i0]->BETY, lat.reseau[i0]->ALFY, lat.reseau[i0]->DY, lat.reseau[i0]->DPY, nsigi, inputfile);
     }
 
@@ -827,5 +828,5 @@
 
 
-    //We take the particles hitting collimators during the previous tracking (we take their coordinates at the beginning of the last turn before the hit), and do a more precise tracking
+    //We take the particles hitting collimators during the previous tracking (we take their coordinates at the beginning of the last turn before the hit), and do a  preciser tracking
     cout <<" "<<endl;
       cout <<"********************************"<<endl;
@@ -1101,4 +1102,5 @@
 
         for (int m(0); m < lat.nhitcolli.size(); ++m) {
+	    cout << "m = " << m << "nhitcollio[m] is: " << nhitcollio[m] << " lat.nhitcolli[m] is:  " << lat.nhitcolli[m] << endl; 
             nhitcollio[m] = nhitcollio[m] + lat.nhitcolli[m];
         }
@@ -1106,4 +1108,5 @@
         for (int m(0); m < lat.nhitspoiler.size(); ++m) {
             nhitspoilero[m] = nhitspoilero[m] + lat.nhitspoiler[m];
+	    
         }
 
@@ -1871,12 +1874,12 @@
 
 
-void Simulation::genbunch(const int& i0, const int& npart, const string& partdistr, const double& r1r2skin, const double& emx, const double& emy, const double& sigdpp, const double& bx, const double& ax, const double& dx, const double& dpx, const double& by, const double& ay, const double& dy, const double& dpy, const double& nsigi)
+void Simulation::genbunch(const int& i0, const int& npart, const string& partdistr, const double& r1r2skin, const double& emx, const double& emy, const double& sigdpp, const double& bx, const double& ax, const double& dx, const double& dpx, const double& by, const double& ay, const double& dy, const double& dpy, const double& nsigi, const string& inputfile)
 {
 
   //save bunch distribution
   
-  //string path1;
-  //path1 = outputpath + "initial.dat";
-  ofstream sortie("initial.dat");
+  string outfile;
+  outfile = inputfile + "initial.dat";
+  ofstream sortie(outfile.c_str());
   
   if (sortie.fail()) {
@@ -1999,5 +2002,5 @@
         sortie << endl;
 
-        sortie << "Number of particles that have hit another element during a preceding turn (asumhits):";
+        sortie << "(Total) Number of particles that have hit another element during the preceding turns (asumhits):";
 
         for (int i(0); i < asumhits.size(); ++i) {
@@ -2007,5 +2010,5 @@
         sortie << endl;
 
-        sortie << "Number of particles that have hit the collimators, for each collimator (nhitcolli):";
+        sortie << "(Total) Number of particles that have hit the collimators, for each collimator (nhitcolli):";
 
         for (int i(0); i < nhitcollio.size(); ++i) {
@@ -2102,5 +2105,5 @@
 
         for (int i(0); i < plosslocal.size(); ++i) {
-            outplot1 << setw(10) <<  plosslocal[i]*PlossPb / NLostTotal << setw(10) << zlosslocal[i] / 1000 << endl;
+            outplot1<< setw(10) <<  plosslocal[i]*PlossPb / NLostTotal << setw(10) << zlosslocal[i] / 1000 << endl;
         }
 
@@ -2190,4 +2193,5 @@
     vector <int> newnhitcolli, ipsnew;
 
+    cout << " size of nhitcollio.size():  "<<nhitcollio.size()<<endl;
     for (int i(0); i < nhitcollio.size(); ++i) {
         if (removedCollimators[i] == 0) {
@@ -2209,5 +2213,6 @@
 
         for (int i(0); i < newnhitcolli.size(); ++i) {
-            outplot3 << setw(10) << i + 1 << setw(10) <<  newnhitcolli[i]*PlossPb / NLostTotal << setw(10) << ipsnew[i] << setw(20) << lat.reseau[ipsnew[i]]->NAME << endl;
+	  cout << newnhitcolli[i]*PlossPb / NLostTotal << endl;
+            outplot3 << setw(10) << i + 1 << setw(20) <<  newnhitcolli[i]*PlossPb / NLostTotal << setw(10) << ipsnew[i] << setw(20) << lat.reseau[ipsnew[i]]->NAME << endl;
         }
         outplot3.close();
@@ -2262,5 +2267,5 @@
 
         for (int i(0); i < asumrem.size(); ++i) {
-            outplot4 << setw(10) << i << setw(10) <<  asumrem[i] << setw(10) << asumhitcolli[i] << setw(10) << asumhits[i] << endl;
+            outplot4 << setw(10) << i << setw(10) <<  asumrem[i] << setw(10) << asumhitcolli[i] << setw(10) << asumhits[i] << setw(10) << (asumrem[i]+asumhitcolli[i]+asumhits[i])<<endl;
         }
         outplot4.close();
@@ -2273,4 +2278,8 @@
     //Writing PlotRunSummary1.p
 
+    //flag to switch to generate GUNPLOT script or not, default is false
+    bool gnuplotflag = true;
+    if(gnuplotflag)
+    {
     ofstream PlotRun1;
     string plotrunname1;
@@ -2288,5 +2297,5 @@
         PlotRun1 << "set key default" << endl;
         PlotRun1 << "set xrange [-1:]" << endl;
-        PlotRun1 << "set ylabel \"P\'(W/m)\"" << endl;
+        PlotRun1 << "set ylabel \"P\' (W/m)\"" << endl;
         PlotRun1 << "set xlabel \"Distance from starting point (km)\"" << endl;
         PlotRun1 << "set title \"Loss map\"" << endl << endl;
@@ -2352,4 +2361,6 @@
     }
 
+    }
+    
 }
 
Index: /ICOSIM/CPP/trunk/source/simulation.h
===================================================================
--- /ICOSIM/CPP/trunk/source/simulation.h	(revision 16)
+++ /ICOSIM/CPP/trunk/source/simulation.h	(revision 17)
@@ -83,5 +83,5 @@
     //method generating an initial bunch of particles
 
-    void genbunch(const int& i0, const int& npart, const string& partdistr, const double& r1r2skin, const double& emx, const double& emy, const double& sigdpp, const double& bx, const double& ax, const double& dx, const double& dpx, const double& by, const double& ay, const double& dy, const double& dpy, const double& nsigi);
+    void genbunch(const int& i0, const int& npart, const string& partdistr, const double& r1r2skin, const double& emx, const double& emy, const double& sigdpp, const double& bx, const double& ax, const double& dx, const double& dpx, const double& by, const double& ay, const double& dy, const double& dpy, const double& nsigi, const string& inputfile);
 
 
