Changeset 17 in ZHANGProjects for ICOSIM/CPP/trunk/source/simcrys.cc

Timestamp:

Sep 30, 2014, 11:05:33 AM (10 years ago)

Author:

zhangj

Message:

Added comments...

File:

• ICOSIM/CPP/trunk/source/simcrys.cc (modified) (15 diffs)

ICOSIM/CPP/trunk/source/simcrys.cc

@@ -576 +576 @@
 * 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 +659 @@
 
 
-
+    
+     
+     
     //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS-------------------------------------------------------------------------
     c_length = crys.Cry_length;
@@ -692 +694 @@
        * 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 +762 @@
       while (!entree.eof()) {
         count = count + 1;
+        
+        //particle hit the crystal or not
+        bool hitflag = false;
+          
         
         getline(entree, rest, ',');
@@ -927 +945 @@
               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 +1080 @@
               //cout<<"----------------------1111111111111111111111111-----------------"<<endl;
               
-        
+           if(layerflag == 0){
+              part.namor = part.namor + 1;     /////////////COUNT/////////////
+              part.effet = 2;
+              }
+              
               if (strncmp(Proc,"out",3)==0) 
                 crossing  = false;
@@ -1069 +1091 @@
               
               if (crossing == true){
-                nhit = nhit + 1;
+                hitflag = true;
+                nhit = nhit + 1;
+                part.nhit = nhit; 
         //      lhit = 100000000*ie + ITURN;
                 impact = x_in0;            
@@ -1261 +1285 @@
                       xp_in0 = xp_shift;
                     }                 ////////////////////////// fin IF NUMERO 17//////////////////
-                    
+                    hitflag = true;
                     nhit = nhit + 1;
+                    part.nhit = nhit; 
                     //lhit = 100000000*ie + ITURN;
                     impact = x_in0;            
@@ -1339 +1364 @@
 
 
-            nabs = 0; 
+            nabs = 0;   //particle lost or not lost.....
             part.effet = nabs;
             //cout<< "debug - S Coll RF 2" ,  s_rot <<endl;
@@ -1375 +1400 @@
              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 +1511 @@
           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 +1569 @@
     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 +1716 @@
     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 +1738 @@
         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 +1752 @@
         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();
     }