Changeset 17 in ZHANGProjects


Ignore:
Timestamp:
Sep 30, 2014, 11:05:33 AM (10 years ago)
Author:
zhangj
Message:

Added comments...

Location:
ICOSIM/CPP/trunk/source
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • ICOSIM/CPP/trunk/source/Collimator.h

    r9 r17  
    5252
    5353
    54     virtual void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath) {};
     54    virtual void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, long int& nhit, string outputpath) {};
    5555    //these methods are specific for a fluka/crystal collimator. They are just defined here.
    5656#if defined(FLUKA)

  • ICOSIM/CPP/trunk/source/CrystalCollimator.cc

    r15 r17  
    3232
    3333
    34 void CrystalCollimator::collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath)
     34void CrystalCollimator::collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, long int& nhit, string outputpath)
    3535{
    36 
     36//pas: number of turns
    3737    //Calculate proton energy
    3838    double Mproton(0.93827231); //rest mass [GeV]
    … …  
    6262
    6363        for (int k(0); k < bunch.size(); ++k) {
    64 
     64            //simulate the particles that are not lost on the previous aperture
     65            if(bunch[k].inabs !=0){   
    6566            sortie  << bunch[k].coordonnees[0][0] << "," << bunch[k].coordonnees[0][2] << "," << bunch[k].coordonnees[0][1] << "," << bunch[k].coordonnees[0][3] << "," << Eproton[k] << endl;
     67            }
    6668        }
    6769
    … …  
    9395   // sim.general(sim, pas, outputpath, NAME, emitx0, emity0,enum, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);
    9496   
    95     sim.general(sim, pas, outputpath, Mirror, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);
     97    sim.general(sim, pas, nhit, outputpath, Mirror, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);
    9698
     99   
     100    //assign the new bunch coordinates after the crystal
     101    for (int k(0); k < bunch.size(); ++k) {
     102      bunch[k].inabs =0 ;
     103    }
     104   
     105   
    97106    ifstream enter;
    98107    string file_in;
    … …  
    105114    } else {
    106115
    107         for (int k(0); k < bunch.size(); ++k) {
     116        for (int k(0); k < bunch.size()-nhit; ++k) {
    108117
    109118            string rest;
    … …  
    120129            getline(enter, rest, ',');
    121130            bunch[k].coordonnees[1][3] = atof(rest.c_str());
    122 
     131           
     132            //particle is not lost (interact) with the crystal
     133            bunch[k].inabs = 1;
     134           
    123135            getline(enter, rest);
    124136            Eproton[k] = atof(rest.c_str());
     137         
    125138        }
    126139
    … …  
    128141    }
    129142
     143   
    130144    vector <double> gamProton;
    131145    vector <double> betProton;

  • ICOSIM/CPP/trunk/source/CrystalCollimator.h

    r5 r17  
    4444    //method specific to the passage through a Crystal collimator
    4545
    46     void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath);
     46    void collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, long int& nhit, string outputpath);
    4747
    4848

  • ICOSIM/CPP/trunk/source/PlotRunSummary1.p

    r5 r17  
    11#Gnuplot script file for the PlotRunSummary (Plot of the loss map)
    22
     3set autoscale xmax
     4set autoscale y
     5set key default
     6set xrange [-1:]
     7set ylabel "P' (W/m)"
     8set xlabel "Distance from starting point (km)"
     9set title "Loss map"
    310
    4    set autoscale xmax
    5    set autoscale y
    6    set key default
    7    set xrange [-1:]
    8    set ylabel "P'(W/m)"
    9    set xlabel "Distance from starting point (km)"
    10    set title "Loss map"
    11 
    12    load 'labelsIP.p'
    13    
    14    plot "PlotRunSummary1.dat" using 2:1 title 'Losses' with boxes, \
    15       "PlotRunSummary2.dat" using 2:1 title 'IP' with impulses
    16 
    17          load 'labelsIPend.p'
     11plot "PlotRunSummary1.dat" using 2:1 title 'Losses' with boxes

  • ICOSIM/CPP/trunk/source/PlotRunSummary2.p

    r5 r17  
    1 #Gnuplot script for the PlotRunSummary (Plot of the collimator load)
     1#Gnuplot script file for the PlotRunSummary (Plot of the collimator load)
    22
     3set autoscale
     4set key default
     5set style fill solid 0.5
     6set ylabel "P_coll (W)"
     7set xlabel "Collimators"
     8set title "Collimator load distribution"
    39
    4    set autoscale
    5    set key default
    6     set title "Collimator load distribution"
    7     set ylabel "P_coll (W)"
    8     set xlabel "Collimators"
    9     set style fill solid 0.5
    10     load 'labelscolli.p'
    11     plot "PlotRunSummary3.dat" using 1:2 title 'Collimator load' with boxes
     10load 'labelscolli.p'
    1211
    13        load 'labelscolliend.p'
     12 plot "PlotRunSummary3.dat" using 1:2 title 'Collimator load' with boxes
     13
     14load 'labelscolliend.p'

  • ICOSIM/CPP/trunk/source/PlotRunSummary3.p

    r5 r17  
    1 #Gnuplot script for PlotRunSummary (Plot of the development after the first impact on collimator)
     1#Gnuplot script file for the PlotRunSummary (Plot of the development after the first impact on collimator)
    22
     3set autoscale
     4set key default
     5set ylabel "Number of particles"
     6set xlabel "N_turn"
     7set title "Time Development after first impact on collimator"
    38
    4    set autoscale
    5    set key default
    6     set title "Time Development after first impact on collimator"
    7     set xlabel "N_turn"
    8     set ylabel "Number of particles"
    9     set logscale y
     9set logscale y
    1010
    11     plot "PlotRunSummary4.dat" using 1:2 title 'particles in beam' with line, \
    12                "PlotRunSummary4.dat" using 1:3 title 'particles lost on collimators' with line, \
    13                "PlotRunSummary4.dat" using 1:4 title 'particles lost elsewhere' with line
     11 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
    1412
    15 
    16                   unset logscale y
     13unset logscale y

  • ICOSIM/CPP/trunk/source/PlotTrajectory1.p

    r5 r17  
    1 #Gnuplot script for PlotTrajectory plot in x
     1#Gnuplot script file for the PlotTrajectory in x
    22
     3set autoscale x
     4set yrange [-25:25]
     5 set style line 1 lw 1
     6set ylabel "x [mm]"
     7set xlabel "S-S_IP1 [km]"
     8set title "Trajectory in x"
    39
    4    set autoscale x
    5    set yrange [-25:25]
    6    set xlabel "S-S_IP1 [km]"
    7    set ylabel "x [mm]"
    8    set title "Trajectory in x"
    9    set style line 1 lw 1
    10    set style line 2 lw 1 lt 6
    11    set key outside below
     10set style line 2 lw 1 lt 6
    1211
     12set key outside below
    1313
    14    load 'plotxco.p'
    15 
    16 
    17    unset key
    18                      
     14load 'plotxco.p'
     15unset key

  • ICOSIM/CPP/trunk/source/PlotTrajectory2.p

    r5 r17  
    1 #Gnuplot script for PlotTrajectory plot in y
     1#Gnuplot script file for the PlotTrajectory in y
    22
     3set autoscale x
     4set yrange [-25:25]
     5 set style line 1 lw 1
     6set ylabel "y [mm]"
     7set xlabel "S-S_IP1 [km]"
     8set title "Trajectory in x"
    39
    4    set autoscale x
    5    set yrange [-25:25]
    6    set xlabel "S-S_IP1 [km]"
    7    set ylabel "y [mm]"
    8    set title "Trajectory in y"
    9    set style line 1 lw 1
    10    set style line 2 lw 1 lt 6
    11    set key outside below
     10set style line 2 lw 1 lt 6
    1211
     12set key outside below
    1313
    14    load 'plotyco.p'
    15 
    16 
    17    unset key
     14load 'plotyco.p'
     15unset key

  • ICOSIM/CPP/trunk/source/crystal_dan_FINAL_VERSION_CHvsVR.f

    r16 r17  
    453453         
    454454c          write(*,*) "y = ", y
    455            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
    456            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
     455c           write(*,*) "before enter VR: x = ", x, "  xp = ",xp
     456c           write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
    457457           
    458458          xp=xp+0.45*(xp/xpcrit+1)*Ang_avr
    … …  
    465465     +      xp ,yp,PC)
    466466     
    467           write(*,*) "after enter VR: x = ", x, "  xp = ",xp
    468           write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
     467c          write(*,*) "after enter VR: x = ", x, "  xp = ",xp
     468c          write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
    469469         
    470470          x = x + 0.5*s_length * xp
    … …  
    501501!     +          DLRi(IS),xp ,yp,PC)
    502502
    503             write(*,*) "before enter VR: x = ", x, "  xp = ",xp
    504             write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
     503c            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
     504c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
    505505           
    506506            CALL CALC_ION_LOSS(IS,PC,s_length-Srefl,DESt)
    507507           
    508             write(*,*) "before enter VR: x = ", x, "  xp = ",xp
    509             write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
    510             write(*,*) "IS = ", IS, " NAM = ", NAM
    511             write(*,*)  "s_length-Srefl = ",s_length-Srefl
    512             write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS)
    513             write(*,*) "DLRi(IS) = ", DLRi(IS)
     508c            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
     509c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
     510c            write(*,*) "IS = ", IS, " NAM = ", NAM
     511c            write(*,*)  "s_length-Srefl = ",s_length-Srefl
     512c            write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS)
     513c            write(*,*) "DLRi(IS) = ", DLRi(IS)
    514514           
    515515            CALL MOVE_AM_(IS,NAM,s_length-Srefl,DESt,DLYi(IS),
    516516     +          DLRi(IS),xp ,yp,PC)
    517517     
    518             write(*,*) "after enter VR: x = ", x, "  xp = ",xp
    519             write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
     518c            write(*,*) "after enter VR: x = ", x, "  xp = ",xp
     519c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
    520520           
    521521            x = x + 0.5 * xp * (s_length - Srefl)
    … …  
    618618c      endif
    619619c     
    620       WRITE(*,*)'Crystal process: ',PROC
     620c      WRITE(*,*)'Crystal process: ',PROC
    621621c      write(*,*)'xp_final :', xp
    622622c      WRITE(*,*)'Crystal process: ',PROC,'Chann Angle',Ch_angle/1000,   +
    … …  
    838838
    839839     
    840       write(*,*) "In Move_AM(): "
     840c      write(*,*) "In Move_AM(): "
    841841
    842842
    843843      xp=xp*1000
    844844      yp=yp*1000
    845       write(*,*)'xp initial:', xp, 'yp initial', yp
     845c      write(*,*)'xp initial:', xp, 'yp initial', yp
    846846C. DEI - dE/dx stoping energy
    847847      PC    = PC - DEI*DZ    ! energy lost because of ionization process[GeV]
    … …  
    849849C. DYA - rms of coloumb scattering
    850850      DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
    851       write(*,*)'dya= ',dya
     851c      write(*,*)'dya= ',dya
    852852     
    853853      ran_x = RAN_GAUSS(1.)
    854       write(*,*) "ran_x = ", ran_x
     854c      write(*,*) "ran_x = ", ran_x
    855855     
    856856      ran_y = RAN_GAUSS(1.)
    857       write(*,*) "ran_y = ", ran_y
     857c      write(*,*) "ran_y = ", ran_y
    858858     
    859859      kxmcs = DYA*ran_x
    … …  
    867867     
    868868c       write(*,*) "ran_x = ", ran_x, "ran_y = ", ran_y
    869       write(*,*) "kxmcs = ", kxmcs, "kymcs = ", kymcs
    870       write(*,*)'xp + kxmcs:', XP, 'yp + kxmcs', YP
     869c      write(*,*) "kxmcs = ", kxmcs, "kymcs = ", kymcs
     870c      write(*,*)'xp + kxmcs:', XP, 'yp + kxmcs', YP
    871871     
    872872cc     XP    = XP+DYA*RAN_GAUSS(1.)
    … …  
    938938      endif
    939939
    940       write(*,*) "ichoix = ", ichoix
    941        write(*,*)'xp final:', xp, 'yp final', yp
     940c      write(*,*) "ichoix = ", ichoix
     941c       write(*,*)'xp final:', xp, 'yp final', yp
    942942       
    943943!---------- calculate the related kick -----------
    … …  
    980980
    981981
    982       write(*,*) "tx = ", tx, " tz = ", tz
     982c      write(*,*) "tx = ", tx, " tz = ", tz
    983983     
    984984      XP = XP + tx
    … …  
    990990      yp=yp/1000
    991991
    992       write(*,*)'xp final:', xp, 'yp final', yp
     992c      write(*,*)'xp final:', xp, 'yp final', yp
    993993!-----------------------------------------------
    994994! print out the energy loss and process experienced

  • ICOSIM/CPP/trunk/source/crystalprotonfuns.f

    r15 r17  
    2828      save
    2929 
    30       write(*,*) "In RAN_GAUSS: "
     30c      write(*,*) "In RAN_GAUSS: "
    3131c      write(*,*) "flag = ",flag
    3232c      write(*,*) "U1 = ", U1,"U2 = ", U2
    … …  
    5353          RAN_GAUSS = x
    5454
    55           write(*,*) "ran_gauss = ", RAN_GAUSS
     55c          write(*,*) "ran_gauss = ", RAN_GAUSS
    5656
    5757        return

  • ICOSIM/CPP/trunk/source/lattice.cc

    r15 r17  
    441441{
    442442
     443  cout << "bunch size is: " << bunch.size() << endl;
     444 
     445 
    443446    if (idpart >= 0) {
    444447        if (idpart > bunch.size()) {
    … …  
    602605          read(bunch);
    603606          }*/
    604 
     607       
     608       
    605609        for (int p(0); p < bunch.size(); ++p) { //loop througt the particles
    606610
    607611            double dpopeff;
    608612
     613            // .inabs = 0, particle lost, .inabs = 1, particle not lost
    609614            if (bunch[p].inabs == 1) {
    610615
    … …  
    743748
    744749                            int lost(bunch.size());
    745 
    746                             resColli[icop]->collipassCrystal(bunch, betgam, turn, outputpath);
    747 
    748                             nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
    749 
     750                             
     751                            long int nhit = 0;
     752                           
     753                            resColli[icop]->collipassCrystal(bunch, betgam, turn, nhit, outputpath);
     754
     755                            //nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
     756                           
     757                             nhitcolli[icop] = nhit;
     758
     759                         //    cout << "after the crystal: bunch size = " << bunch.size()<< endl;
    750760                        }
    751761
    … …  
    924934                        if (inside == false) {
    925935                            cout <<"The particle is lost at element " << reseau[i]->NAME << endl;
     936//                           cout << "Particle ID is: " << p << endl;
     937//                           cout << "bunch size is: " << bunch.size() << endl;
     938//                           cout << "bunch status is: " << bunch[p].inabs << endl;
     939//                           
    926940                            xh0 = bunch[p].coordonnees[0][0];
    927941                            bunch[p].inabs = 0;
    … …  
    11001114            }
    11011115        }
     1116    //    cout << "bunch size is ...." << bunch.size() << endl;
    11021117
    11031118        if (plotflag == "Yes") {
    … …  
    11321147    vector <Particle> tempinabs;
    11331148
    1134     //we contine just with particles that are not lost
     1149    //we contine to track for the next turn; just with particles that are not lost
    11351150    for (int w(0); w < bunch.size(); ++w) {
    11361151        if (bunch[w].inabs != 0) {
    … …  
    11391154    }
    11401155
     1156 //   cout << "bunch size before the turn: "<< bunch.size() << endl;
     1157   
    11411158    bunch.clear();
    11421159
    … …  
    11471164    tempinabs.clear();
    11481165
    1149 }
     1166   // cout << "bunch size after the turn " <<bunch.size() << endl;
     1167}

  • ICOSIM/CPP/trunk/source/partcrys.h

    r15 r17  
    4444    int nvcapt;
    4545    int nvrefl;
     46    int nabsorbed;
    4647    int nout;          //number of particles that not hit the crystal collimator
    47     int effet;        //Number that defini the effect :1 --> Out ;;; 2 --> Amorphous ;;; 3 --> Channeling ;;; 4 --> Dechanneling ;;; 5 --> V. Reflection ;;; 6 --> V. Capture
     48    int effet;        //Number that defini the effect :1 --> Out ;;; 2 --> Amorphous ;;; 3 --> Channeling ;;; 4 --> Dechanneling ;;; 5 --> V. Reflection ;;; 6 --> V. Capture, 7 --> absorbed
     49   long int nhit; //number of particles that hit the crystal collimator
     50   
    4851
    4952};

  • ICOSIM/CPP/trunk/source/simcrys.cc

    r16 r17  
    576576* original crystal routine
    577577* *******************************************************************************/
    578 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)
     578void 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)
    579579{
    580580    srand ( time(NULL) );
    … …  
    659659
    660660
    661 
     661   
     662     
     663     
    662664    //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS-------------------------------------------------------------------------
    663665    c_length = crys.Cry_length;
    … …  
    692694       * xmin=x-, xmax=x+, etc. from the input file
    693695       */
    694      long int    nhit    = 0;
     696   
     697
     698      //initialize particle variables
     699      part.n_part = 0;
     700      part.nabsorbed = 0;
     701      part.namor = 0;
     702      part.nchann = 0;
     703      part.ndechann = 0;
     704      part.nout =0;
     705      part.nvcapt = 0;
     706      part.nvrefl = 0;
     707      part.nhit = 0;
     708      nhit    = 0;
    695709    double fracab  = 0;
    696710  //  int   n_chan  = 0;          // valentina :initialize to zero the counters for crystal effects
    … …  
    748762      while (!entree.eof()) {
    749763        count = count + 1;
     764       
     765        //particle hit the crystal or not
     766        bool hitflag = false;
     767         
    750768       
    751769        getline(entree, rest, ',');
    … …  
    927945              s_impact = s_in0;       //!(for the first impact)
    928946             
    929               if(1){
     947              if(0){
    930948                cout<<"hit the cry entrance face"<<endl;
    931949                cout<<"impact at s = "<< s_impact<< ",  x = "<<x_in0<<endl;
    … …  
    10621080              //cout<<"----------------------1111111111111111111111111-----------------"<<endl;
    10631081             
    1064        
     1082           if(layerflag == 0){
     1083              part.namor = part.namor + 1;     /////////////COUNT/////////////
     1084              part.effet = 2;
     1085              }
     1086             
    10651087              if (strncmp(Proc,"out",3)==0)
    10661088                crossing  = false;
    … …  
    10691091             
    10701092              if (crossing == true){
    1071                 nhit = nhit + 1;
     1093                hitflag = true;
     1094                nhit = nhit + 1;
     1095                part.nhit = nhit;
    10721096        //      lhit = 100000000*ie + ITURN;
    10731097                impact = x_in0;           
    … …  
    12611285                      xp_in0 = xp_shift;
    12621286                    }                 ////////////////////////// fin IF NUMERO 17//////////////////
    1263                    
     1287                    hitflag = true;
    12641288                    nhit = nhit + 1;
     1289                    part.nhit = nhit;
    12651290                    //lhit = 100000000*ie + ITURN;
    12661291                    impact = x_in0;           
    … …  
    13391364
    13401365
    1341             nabs = 0;
     1366            nabs = 0;   //particle lost or not lost.....
    13421367            part.effet = nabs;
    13431368            //cout<< "debug - S Coll RF 2" ,  s_rot <<endl;
    … …  
    13751400             part.nvcapt = part.nvcapt + 1;     /////////////COUNT/////////////
    13761401            part.effet = 6;}
     1402             if(strncmp(Proc,"absorbed",8==0)){
     1403             part.nabsorbed = part.nabsorbed +1;
     1404              part.effet = 7;
     1405            }
     1406             /*           
     1407 if (PROC(1:3).eq.'pne')
     1408                 PROC_dan=7
     1409             if (PROC(1:3).eq.'ppe')
     1410                 PROC_dan=8
     1411               if (PROC(1:4).eq.'diff')
     1412                 PROC_dan=9
     1413               if (PROC(1:4).eq.'ruth')
     1414                 PROC_dan=10
     1415             */
     1416   
    13771417           
    13781418              //  ?????????????????????????????????????
    … …  
    14711511          av_pc1 = av_pc1 + part.PC;
    14721512
    1473 
     1513      if(0){
    14741514          cout<<"  X apres :"<<part.x<<"  XP apres :"<<part.xp<<"  Y apres :"<<part.y<<"  YP apres :"<<part.yp<<"  PC apres :"<<part.PC<<endl;
    14751515          //cout<<endl;
     1516      }
     1517      //only record the particles that hit the collimator.
     1518      if(hitflag != true){
    14761519          //sortieIco.setf(ios::scientific);
    1477           sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl;
     1520          sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC  << endl;
    14781521          //}//end of for....
    1479          
     1522      }
     1523     
    14801524        }
    14811525      } //End of while
    … …  
    15251569    file_out(pas, outputpath);
    15261570   
     1571    if(1){
     1572      cout << "number of particles before the crystal: " << count<< endl;
     1573      cout << "number of particle hit the crystal: " << part.nhit << endl;
     1574      cout << "number of particles after the crystal: " << count - part.nhit << endl;
     1575    }
    15271576    //cout<<endl;
    15281577    //cout<<endl;
    … …  
    16671716    cout << "N.OUT :" << part.nout << endl;
    16681717    cout << "N.VREFLECTION :" << part.nvrefl << endl;
    1669     cout << "N.VCAPTURE :" << part.nvcapt << endl;
     1718    cout << "N.VCAPTURE:" << part.nvcapt << endl;
     1719    cout << "N.ABSORBED: "<< part.nabsorbed << endl;
    16701720
    16711721    //cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!";
    … …  
    16881738        out << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
    16891739        out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
    1690         out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
     1740        out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << 1e6*crys.Cry_length / crys.Rcurv << endl;
    16911741        out << "Critical Radius :" << crys.Rcrit << endl;
    16921742        out << "Ratio :" << crys.ratio << endl;
    … …  
    17021752        out << "N.OUT :" << part.nout << endl;
    17031753        out << "N.VREFLECTION :" << part.nvrefl << endl;
    1704         out << "N.VCAPTURE :" << part.nvcapt << endl << endl << endl;
    1705 
     1754        out << "N.VCAPTURE :" << part.nvcapt << endl ;
     1755        out << "N.ABSORBED: "<< part.nabsorbed << endl;
     1756        out << "N. of particles that hit the crystal collimator: " << part.nhit << endl;
     1757        cout << endl;
     1758        cout << endl;
     1759       
    17061760        out.close();
    17071761    }

  • ICOSIM/CPP/trunk/source/simcrys.h

    r9 r17  
    5353    //THE GENERAL FUNCTION OF THE TEST THAT REGROUP ALL THE PREVIOUS FUNCTIONS!!!!
    5454
    55     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);
     55    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);
    5656
    5757

  • ICOSIM/CPP/trunk/source/simulation.cc

    r15 r17  
    7878    particlefile = inputfile + "initial.dat";
    7979
     80    cout << "READ THE BUNCH OF PARTICLES FROM THE FILE: " << particlefile <<endl;
     81   
    8082    //if the name of the file is not 'initial.dat', uncomment the following two lines
    8183
    … …  
    663665      cout <<"********************************"<<endl;
    664666        cout << "RunningFlag = " <<  RunningFlag << endl;
    665         cout << "READ THE BUNCH OF PARTICLES FROM THE FILE: " <<inputfile <<endl;
     667       
    666668        readParticle(Apr, Zpr, inputfile);
    667669
    … …  
    676678      cout << "RunningFlag = " <<  RunningFlag << endl;
    677679      cout << "GENERATION OF THE BUNCH OF PARTICLES." << endl;
    678         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);
     680        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);
    679681
    680682        //we set the id of the particles
    … …  
    742744      cout <<"********************************"<<endl;
    743745        cout << "RunningFlag = " <<  RunningFlag << endl;
    744         cout << "READ THE BUNCH OF PARTICLES FROM THE FILE: " <<inputfile <<endl;
    745746        readParticle(Apr, Zpr, inputfile);
    746747    } else if (RunningFlag == 3) {
    … …  
    749750      cout << "RunningFlag = " <<  RunningFlag << endl;
    750751      cout << "GENERATION OF THE BUNCH OF PARTICLES." << endl;
    751         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);
     752        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);
    752753    }
    753754
    … …  
    827828
    828829
    829     //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
     830    //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
    830831    cout <<" "<<endl;
    831832      cout <<"********************************"<<endl;
    … …  
    11011102
    11021103        for (int m(0); m < lat.nhitcolli.size(); ++m) {
     1104            cout << "m = " << m << "nhitcollio[m] is: " << nhitcollio[m] << " lat.nhitcolli[m] is:  " << lat.nhitcolli[m] << endl;
    11031105            nhitcollio[m] = nhitcollio[m] + lat.nhitcolli[m];
    11041106        }
    … …  
    11061108        for (int m(0); m < lat.nhitspoiler.size(); ++m) {
    11071109            nhitspoilero[m] = nhitspoilero[m] + lat.nhitspoiler[m];
     1110           
    11081111        }
    11091112
    … …  
    18711874
    18721875
    1873 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)
     1876void 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)
    18741877{
    18751878
    18761879  //save bunch distribution
    18771880 
    1878   //string path1;
    1879   //path1 = outputpath + "initial.dat";
    1880   ofstream sortie("initial.dat");
     1881  string outfile;
     1882  outfile = inputfile + "initial.dat";
     1883  ofstream sortie(outfile.c_str());
    18811884 
    18821885  if (sortie.fail()) {
    … …  
    19992002        sortie << endl;
    20002003
    2001         sortie << "Number of particles that have hit another element during a preceding turn (asumhits):";
     2004        sortie << "(Total) Number of particles that have hit another element during the preceding turns (asumhits):";
    20022005
    20032006        for (int i(0); i < asumhits.size(); ++i) {
    … …  
    20072010        sortie << endl;
    20082011
    2009         sortie << "Number of particles that have hit the collimators, for each collimator (nhitcolli):";
     2012        sortie << "(Total) Number of particles that have hit the collimators, for each collimator (nhitcolli):";
    20102013
    20112014        for (int i(0); i < nhitcollio.size(); ++i) {
    … …  
    21022105
    21032106        for (int i(0); i < plosslocal.size(); ++i) {
    2104             outplot1 << setw(10) <<  plosslocal[i]*PlossPb / NLostTotal << setw(10) << zlosslocal[i] / 1000 << endl;
     2107            outplot1<< setw(10) <<  plosslocal[i]*PlossPb / NLostTotal << setw(10) << zlosslocal[i] / 1000 << endl;
    21052108        }
    21062109
    … …  
    21902193    vector <int> newnhitcolli, ipsnew;
    21912194
     2195    cout << " size of nhitcollio.size():  "<<nhitcollio.size()<<endl;
    21922196    for (int i(0); i < nhitcollio.size(); ++i) {
    21932197        if (removedCollimators[i] == 0) {
    … …  
    22092213
    22102214        for (int i(0); i < newnhitcolli.size(); ++i) {
    2211             outplot3 << setw(10) << i + 1 << setw(10) <<  newnhitcolli[i]*PlossPb / NLostTotal << setw(10) << ipsnew[i] << setw(20) << lat.reseau[ipsnew[i]]->NAME << endl;
     2215          cout << newnhitcolli[i]*PlossPb / NLostTotal << endl;
     2216            outplot3 << setw(10) << i + 1 << setw(20) <<  newnhitcolli[i]*PlossPb / NLostTotal << setw(10) << ipsnew[i] << setw(20) << lat.reseau[ipsnew[i]]->NAME << endl;
    22122217        }
    22132218        outplot3.close();
    … …  
    22622267
    22632268        for (int i(0); i < asumrem.size(); ++i) {
    2264             outplot4 << setw(10) << i << setw(10) <<  asumrem[i] << setw(10) << asumhitcolli[i] << setw(10) << asumhits[i] << endl;
     2269            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;
    22652270        }
    22662271        outplot4.close();
    … …  
    22732278    //Writing PlotRunSummary1.p
    22742279
     2280    //flag to switch to generate GUNPLOT script or not, default is false
     2281    bool gnuplotflag = true;
     2282    if(gnuplotflag)
     2283    {
    22752284    ofstream PlotRun1;
    22762285    string plotrunname1;
    … …  
    22882297        PlotRun1 << "set key default" << endl;
    22892298        PlotRun1 << "set xrange [-1:]" << endl;
    2290         PlotRun1 << "set ylabel \"P\'(W/m)\"" << endl;
     2299        PlotRun1 << "set ylabel \"P\' (W/m)\"" << endl;
    22912300        PlotRun1 << "set xlabel \"Distance from starting point (km)\"" << endl;
    22922301        PlotRun1 << "set title \"Loss map\"" << endl << endl;
    … …  
    23522361    }
    23532362
     2363    }
     2364   
    23542365}
    23552366

  • ICOSIM/CPP/trunk/source/simulation.h

    r8 r17  
    8383    //method generating an initial bunch of particles
    8484
    85     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);
     85    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);
    8686
    8787
Note: See TracChangeset for help on using the changeset viewer.