*$ CREATE COMSCW.FOR
*COPY COMSCW
*
*=== comscw  ===========================================================*
*
      DOUBLE PRECISION FUNCTION COMSCW ( IJ    , XA    , YA    , ZA    ,
     &                                   MREG  , RULL  , LLO   , ICALL )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1989-2005      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     New comscw for FLUKA9x-200x                                      *
*                                                                      *
*     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
*     !!! This is a completely dummy routine for Fluka9x/200x. !!!     *
*     !!! The  name has been kept the same as for older  Fluka !!!     *
*     !!! versions for back-compatibility, even though  Comscw !!!     *
*     !!! is applied only to estimators which didn't exist be- !!!     *
*     !!! fore Fluka89.                                        !!!     *
*     !!! User  developed versions  can be used for  weighting !!!     *
*     !!! density-like quantities at runtime                   !!!     *
*     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
*                                                                      *
*     Input variables:                                                 *
*                                                                      *
*           Ij = (generalized) particle code                           *
*     Xa,Ya,Za = position                                              *
*         Mreg = region number                                         *
*         Rull = amount to be deposited                                *
*          Llo = particle generation                                   *
*        Icall = call id                                               *
*                                                                      *
*     Output variables:                                                *
*                                                                      *
*       Comscw = factor the scored amount will be multiplied by        *
*       Lsczer = logical flag, if true no amount will be scored        *
*                regardless of Comscw                                  *
*                                                                      *
*     Useful variables (common SCOHLP):                                *
*                                                                      *
*     Energy/Star binnings/scorings (Comscw):                          *
*          ISCRNG = 1 --> Energy density  binning                      *
*          ISCRNG = 2 --> Star   density  binning                      *
*          ISCRNG = 3 --> Residual nuclei scoring                      *
*          ISCRNG = 4 --> Momentum transfer density binning            *
*          ISCRNG = 5 --> Activity density binning                     *
*          JSCRNG = # of the binning                                   *
*                                                                      *
*     Useful variables (common SOUEVT):                                *
*                                                                      *
*          X,Y,Zsoevt(i) = position    of the i_th source particle     *
*          TX,Y,Zsoev(i) = direction   of the i_th source particle     *
*              Wtsoev(i) = weight      of the i_th source particle     *
*              Pmsoev(i) = momentum    of the i_th source particle     *
*              Tksoev(i) = kin. energy of the i_th source particle     *
*              Agsoev(i) = age         of the i_th source particle     *
*              Aksoev(i) = Kaon ampl.  of the i_th source particle     *
*              Ussoev(i) = user var.   of the i_th source particle     *
*              Ijsoev(i) = identity    of the i_th source particle     *
*              Nrsoev(i) = region      of the i_th source particle     *
*              Nlsoev(i) = lattice     of the i_th source particle     *
*                Npsoev  = number of the source particles              *
*                                                                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(FLKMAT)'
      INCLUDE '(SCOHLP)'
      INCLUDE '(USRBIN)'
      INCLUDE '(SOUEVT)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(TRACKR)'
      INCLUDE '(PAPROP)'
*     

       LOGICAL LFIRST
       SAVE LFIRST
       DATA LFIRST /.TRUE./

       
       REAL X,Y
       REAL X_A, X_Z
       REAL ENERGY
       REAL X_AMU
       REAL X_P

*       X_AMU = 938.271996
       X_AMU = 931.494088

       X = 0;
       Y = 0;

* ../../flutil/ldpm3qmd comscw.f

* return message from first call
       IF (LFIRST) THEN
          WRITE(LUNOUT,*) 'COMSCW TEST_NB'
          LFIRST = .FALSE.
       END IF
       COMSCW=ZERZER

*	WRITE(LUNOUT,*) 'COMSCW TEST_NB *****', TITUSB(JSCRNG), '****'


* All ions:
       IF (TITUSB(JSCRNG) .EQ. 'ENERGY  ') THEN
          COMSCW=ONEONE
*          WRITE(LUNOUT,*) 'COMSCW TEST_NB ALL'
       END IF
*************************** FUN STARTS HERE ********
* Translate into resonse
       IF (TITUSB(JSCRNG) .EQ. 'ALANINE ') THEN
*     debugging only

*	  WRITE(LUNOUT,*) 'COMSCW TEST_NB _INSIDEAL__', TITUSB(JSCRNG)
          
          X_P = Ptrack * 1000.0
          
*          ENERGY = SQRT(X_P*X_P + X_AMU*4*X_AMU*4) 
*     $         - X_AMU*4
*
*          WRITE(LUNOUT,*) '_NB ',TITUSB(JSCRNG),JTRACK,
*     $         X_P, TKSOEV(NPSOEV)*1000, ENERGY
*          WRITE(LUNOUT,*) '_NB ',TITUSB(JSCRNG),JTRACK,
*     $         X_P
* deposit all dose as default with modifications below

          COMSCW=ONEONE

*     Heavy ions with Z > 2:
          IF (JTRACK .EQ. -2) THEN
             CALL USRDCI(JTRACK,IONA,IONZ,IONM)
             
*     translate to log(MeV/nucleon)

             ENERGY = (SQRT(X_P*X_P + X_AMU*IONA*X_AMU*IONA) 
     $            - X_AMU*IONA)/IONA             
             IF (ENERGY .GT. 0) THEN
                X = LOG10(ENERGY)
             END IF

* y = -0.41578 + 0.3092 * x - 0.074314 * x^2 + 0.14397 * x^3 - 0.1506 * x^4 + 0.055563 * x^5 - 0.0069034 * x^6

             IF (X .LT. 3) THEN
                IF (X .GT. -1.03) THEN
                   IF (IONZ .EQ. 3) THEN
                      Y = -0.41578 + 0.3092 * X - 0.074314 * X**2 + 
     $                     0.14397 * X**3 - 0.1506 * X**4 + 
     $                     0.055563 * x**5 - 0.0069034 * x**6
                      COMSCW=10**Y
                   END IF   


* y = -0.45724 + 0.28374 * x - 0.15943 * x^2 + 0.23361 * x^3 - 0.13681 * x^4 + 0.028907 * x^5 - 0.0016492 * x^6

                   IF (IONZ .EQ. 4) THEN
                      Y = -0.45724 + 0.28374 * X - 0.15943 * X**2 + 
     $                     0.23361 * X**3 - 0.13681 * X**4 + 
     $                     0.028907 * x**5 - 0.0016492 * x**6
                      COMSCW=10**Y
                   END IF   

* y = -0.48919 + 0.27817 * x - 0.21775 * x^2 + 0.27529 * x^3 - 0.12609 * x^4 + 0.017352 * x^5 + 0.00025007 * x^6

                   IF (IONZ .EQ. 5) THEN
                      Y = -0.48919 + 0.27817 * X - 0.21775 * X**2 + 
     $                     0.27529 * X**3 - 0.12609 * X**4 + 
     $                     0.017352 * x**5 + 0.00025007 * x**6
                      COMSCW=10**Y
                   END IF   
                   
* y = -0.51196 + 0.28205 * x - 0.25975 * x^2 + 0.295 * x^3 - 0.12248 * x^4 + 0.015699 * x^5 + 0.00019709 * x^6
                   IF (IONZ .EQ. 6) THEN 
                      Y = -0.51196 + 0.28205 * X - 0.25975 * X**2 + 
     $                     0.295 * X**3 - 0.12248 * X**4 + 
     $                     0.015699 * x**5 + 0.00019709 * x**6
                      COMSCW=10**Y
                   END IF
                END IF   

***                IF (Y .GT. 0) THEN                   
***                   WRITE(LUNOUT,*) '__NB ETA ',TITUSB(JSCRNG),JTRACK,
***     $                  ENERGY, X, 10**Y
***                END IF
             END IF
             IF (X .GE. 3) THEN
                COMSCW = ONEONE                          
             END IF         
             
*     estimate for low energies for Z > 2:
             IF (X .LE. -1.03) THEN
                COMSCW = 0.3
             END IF         
          END IF
*         end of Z > 2 

*     estimate for Fission fragments (208)
          IF (JTRACK .EQ. 208) THEN
             COMSCW = 0.3
          END IF
          
          X_A = 0
          X_Z = 0
*     indicidual scoring: 
*     antiprotons and protons
          IF ((JTRACK .EQ. 2) .OR. (JTRACK .EQ. 1)) THEN
             X_A = 1
             X_Z = 1
          END IF
* Deuterons 
          IF (JTRACK .EQ. -3) THEN
             X_A = 2
             X_Z = 1
          END IF
* Tritons
          IF (JTRACK .EQ. -4) THEN
             X_A = 3
             X_Z = 1
          END IF
* He3
          IF (JTRACK .EQ. -5) THEN
             X_A = 3
             X_Z = 2
          END IF

* He4
          IF (JTRACK .EQ. -6) THEN
             X_A = 4
             X_Z = 2
          END IF
          
          IF (X_A .NE. 0) THEN 

             ENERGY = (SQRT(X_P*X_P + X_AMU*X_A*X_AMU*X_A) 
     $            - X_AMU*X_A)/X_A
             
             IF (ENERGY .GT. 0) THEN
                X = LOG10(ENERGY)
             END IF

             IF (X .GT. -1.03) THEN
* Z = 1         
* y = -0.11157 + 0.20489 * x - 0.17707 * x^2 + 0.10375 * x^3 - 0.041892 * x^4 + 0.0099014 * x^5 - 0.0009912 * x^6
    
                IF (X_Z .EQ. 1) THEN
                   Y = -0.11157 + 0.20489 * X - 0.17707 * X**2 + 
     $                  0.10375 * X**3 - 0.041892 * X**4 + 
     $                  0.0099014 * x**5 - 0.0009912 * x**6
                   COMSCW=10**Y   
                END IF

*     Z = 2             
* y = -0.30799 + 0.34219 * x - 0.13186 * x^2 + 0.098529 * x^3 
* - 0.1059 * x^4 + 0.045992 * x^5 - 0.0065743 * x^6
* y = -0.3078 + 0.34174 * x - 0.13337 * x^2 + 0.10006 * x^3 - 0.1047 * x^4 + 0.044674 * x^5 - 0.0063064 * x^6
                IF (X_Z .EQ. 2) THEN
                   Y = -0.3078 + 0.34219 * X - 0.13186 * X**2 + 
     $                  0.098529 * X**3 - 0.1059 * X**4 + 
     $                  0.045992 * x**5 - 0.0065743 * x**6
                   COMSCW=10**Y
                END IF                     
             END IF

***             IF (Y .GT. 0) THEN                   
***                WRITE(LUNOUT,*) '__NB ETA ',TITUSB(JSCRNG),JTRACK,
***     $               ENERGY, X, 10**Y
***             END IF


*             WRITE(LUNOUT,*) '__NB ',TITUSB(JSCRNG),JTRACK,X_Z,
*     $            ENERGY, X, 10**Y

             IF (X .GE. 3) THEN
                COMSCW = ONEONE                          
             END IF         
             
*     estimate for low energies for Z > 2:
             IF (X .LE. -1.03) THEN
                COMSCW = 0.3
             END IF 
          END IF
*     END OF Z = 1 and Z = 2
       END IF 
*     END OF ALANINE SECTION
       LSCZER = .FALSE.
       RETURN
*===  End of function Comscw.f ===========================================*
       END



