[no commit message]

.hgignore

@@ -0,0 +1,5 @@
+\.orig$
+\.orig\..*$
+\.chg\..*$
+\.rej$
+\.conflict\~$

nbactions.xml

@@ -0,0 +1,37 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<actions>
+        <action>
+            <actionName>run</actionName>
+            <goals>
+                <goal>process-classes</goal>
+                <goal>org.codehaus.mojo:exec-maven-plugin:1.2.1:exec</goal>
+            </goals>
+            <properties>
+                <exec.args>-ea -classpath %classpath hep.dataforge.trapping.Trapping D:\\PlayGround\\Trapping.res</exec.args>
+                <exec.executable>java</exec.executable>
+            </properties>
+        </action>
+        <action>
+            <actionName>debug</actionName>
+            <goals>
+                <goal>process-classes</goal>
+                <goal>org.codehaus.mojo:exec-maven-plugin:1.2.1:exec</goal>
+            </goals>
+            <properties>
+                <exec.args>-Xdebug -Xrunjdwp:transport=dt_socket,server=n,address=${jpda.address} -ea -classpath %classpath ${packageClassName} D:\\PlayGround\\Trapping.res</exec.args>
+                <exec.executable>java</exec.executable>
+                <jpda.listen>true</jpda.listen>
+            </properties>
+        </action>
+        <action>
+            <actionName>profile</actionName>
+            <goals>
+                <goal>process-classes</goal>
+                <goal>org.codehaus.mojo:exec-maven-plugin:1.2.1:exec</goal>
+            </goals>
+            <properties>
+                <exec.args>${profiler.args} -ea -classpath %classpath ${packageClassName} D:\\PlayGround\\Trapping.res</exec.args>
+                <exec.executable>${profiler.java}</exec.executable>
+            </properties>
+        </action>
+    </actions>

pom.xml

@@ -0,0 +1,57 @@
+<project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+  xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">
+  <modelVersion>4.0.0</modelVersion>
+
+  <groupId>hep.DataForge</groupId>
+  <artifactId>Trapping</artifactId>
+  <version>1.0-SNAPSHOT</version>
+  <packaging>jar</packaging>
+
+  <name>Trapping</name>
+  <url>http://maven.apache.org</url>
+
+  <properties>
+    <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
+  </properties>
+
+
+<build>
+    <plugins>
+        <plugin>
+            <groupId>org.apache.maven.plugins</groupId>
+            <artifactId>maven-resources-plugin</artifactId>
+            <version>2.6</version>
+                
+        </plugin>
+      <plugin>
+        <groupId>org.apache.maven.plugins</groupId>
+        <artifactId>maven-compiler-plugin</artifactId>
+        <version>2.3.2</version>
+        <configuration>
+          <source>1.7</source>
+          <target>1.7</target>
+          <showDeprecation>true</showDeprecation>
+        </configuration>
+      </plugin>
+    </plugins>
+</build>
+
+  <dependencies>
+    <dependency>
+      <groupId>junit</groupId>
+      <artifactId>junit</artifactId>
+      <version>3.8.1</version>
+      <scope>test</scope>
+    </dependency>
+    <dependency>
+      <groupId>org.apache.commons</groupId>
+      <artifactId>commons-math3</artifactId>
+      <version>3.2</version>
+    </dependency>
+    <dependency>
+      <groupId>net.java.dev.jna</groupId>
+      <artifactId>jna</artifactId>
+      <version>3.5.2</version>
+    </dependency>
+  </dependencies>
+</project>

src/main/java/hep/dataforge/trapping/DoubleValue.java

@@ -0,0 +1,33 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package hep.dataforge.trapping;
+
+/**
+ * Класс нужен исключительно чтобы сделать простой доступ к Сишным экспортам
+ * @author Darksnake
+ */
+public class DoubleValue {
+    private double value;
+
+    public DoubleValue(double value) {
+        this.value = value;
+    }
+
+    
+    
+    /**
+     * @return the value
+     */
+    public double getValue() {
+        return value;
+    }
+
+    /**
+     * @param value the value to set
+     */
+    public void setValue(double value) {
+        this.value = value;
+    }
+}

src/main/java/hep/dataforge/trapping/ElectronTrappingSimulator.java

@@ -0,0 +1,244 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package hep.dataforge.trapping;
+
+import java.util.ArrayList;
+import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
+import org.apache.commons.math3.geometry.euclidean.threed.SphericalCoordinates;
+import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
+import org.apache.commons.math3.util.Precision;
+
+/**
+ *
+ * @author Darksnake
+ */
+public class ElectronTrappingSimulator {
+
+    private TrappingRandomGenerator generator = new TrappingRandomGenerator();
+    double Elow = 14000d;
+    double thetaTransport = 24.107064 / 180 * Math.PI;
+    double thetaPinch = 19.481097 / 180 * Math.PI;
+
+    public ElectronTrappingSimulator() {
+    }
+
+    public static enum EndState {
+
+        ACCEPTED,//трэппинговый электрон попал в аксептанс
+        REJECTED,//трэппинговый электрон вылетел через заднюю пробку
+        LOWENERGY,//потерял слишком много энергии
+        PASS,//электрон никогда не запирался и прошел напрямую, нужно для нормировки
+        NONE
+    }
+
+    public void setFields(double Bsource, double Btransport, double Bpinch) {
+        this.thetaTransport = Math.asin(Math.sqrt(Bsource / Btransport));
+        this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch));
+    }
+
+    /**
+     * Симулируем один пробег электрона от начального значения и до вылетания из
+     * иточника или до того момента, как энергия становится меньше Elow.
+     * Считаем, что за один проход источника происходит в среднем существенно
+     * меньше одного столкновения, поэтому выход вперед или назад совершенно
+     * симментричны.
+     *
+     */
+    public SimulaionResult simulateOne(double initEnergy, double initTheta) {
+        assert initEnergy > 0;
+        assert initTheta > 0 && initTheta < Math.PI / 2;
+
+        if (initTheta < this.thetaPinch) {
+            if (generator.heads()) {
+                return new SimulaionResult(EndState.PASS, initEnergy, initTheta, 0);
+            } else {
+                return new SimulaionResult(EndState.REJECTED, initEnergy, initTheta, 0);
+            }
+        } else if(initTheta<this.thetaTransport){
+                return new SimulaionResult(EndState.REJECTED, initEnergy, initTheta, 0);            
+        }
+
+
+        double E = initEnergy;
+        double theta = initTheta;
+        int colNum = 0;
+        EndState state = EndState.NONE;
+
+        boolean stopflag = false;
+
+        while (!stopflag) {
+            colNum++;
+            DoubleValue dE = new DoubleValue(0);
+            DoubleValue dTheta = new DoubleValue(0);
+
+
+            //Вычисляем сечения и нормируем их на полное сечение
+//            double sigmaTotal = Scatter.sigmaTotal(E);
+//            double sigmaIon = Scatter.sigmaion(E) / sigmaTotal;
+//            double sigmaEl = Scatter.sigmael(E) / sigmaTotal;
+//            double sigmaexc = Scatter.sigmaexc(E) / sigmaTotal;
+            double sigmaIon = Scatter.sigmaion(E);
+            double sigmaEl = Scatter.sigmael(E);
+            double sigmaexc = Scatter.sigmaexc(E);
+            double sigmaTotal = sigmaEl + sigmaIon + sigmaexc;
+            sigmaIon /= sigmaTotal;
+            sigmaEl /= sigmaTotal;
+            sigmaexc /= sigmaTotal;
+            //проверяем нормировку
+            assert Precision.equals(sigmaEl + sigmaexc + sigmaIon, 1, 1e-2);
+
+            double alpha = generator.next();
+
+            if (alpha > sigmaEl) {
+                if (alpha > sigmaEl + sigmaexc) {
+                    //ionization case
+                    Scatter.randomion(E, dE, dTheta);
+                } else {
+                    //excitation case
+                    Scatter.randomexc(E, dE, dTheta);
+                }
+            } else {
+                // elastic
+                Scatter.randomel(E, dE, dTheta);
+            }
+
+//            if (alpha < sigmaEl) {
+//                Scatter.randomel(E, dE, dTheta);
+//            } else if (alpha < sigmaexc) {
+//                Scatter.randomexc(E, dE, dTheta);
+//            } else {
+//                Scatter.randomion(E, dE, dTheta);
+//            }
+
+
+
+
+            //Обновляем значени угла и энергии независимо ни от чего
+            E -= dE.getValue();
+            //Изменение угла
+            theta = addTheta(theta, dTheta.getValue() / 180 * Math.PI);
+            //следим чтобы угол был от 0 до 90, если он перекинется через границу, считаем что электрон остается в потоке
+            theta = normalizeTheta(theta);
+
+            if (theta < thetaPinch) {
+                stopflag = true;
+                if (generator.heads()) {
+                    //Учитываем тот факт, что электрон мог вылететь в правильный угол, но назад
+                    state = EndState.ACCEPTED;
+                } else {
+                    state = EndState.REJECTED;
+                }
+            }
+
+            if (theta >= thetaPinch && theta <= thetaTransport) {
+                //Вылет через заднюю пробку эквивалентен отражению от передней.
+                //это верно при низких концентрациях
+                stopflag = true;
+                state = EndState.REJECTED;
+            }
+
+            if (E < Elow) {
+                //Если энергия стала слишком маленькой
+                stopflag = true;
+                state = EndState.LOWENERGY;
+            }
+        }
+        return new SimulaionResult(state, E, theta, colNum);
+
+    }
+
+    /**
+     * Сложение вектора с учетом случайно распределения по фи
+     *
+     * @param theta
+     * @param dTheta
+     * @return
+     */
+    double addTheta(double theta, double dTheta) {
+        //Генерируем случайный фи
+        double phi = generator.next() * 2 * Math.PI;
+        //Создаем начальный вектор в сферических координатах
+        SphericalCoordinates init = new SphericalCoordinates(1, 0, theta+dTheta);
+        // Задаем вращение относительно оси, перпендикулярной исходному вектору 
+        SphericalCoordinates rotate = new SphericalCoordinates(1, 0, theta);
+        // поворачиваем исходный вектора на dTheta
+        Rotation rot = new Rotation(rotate.getCartesian(), phi);
+
+        Vector3D result = rot.applyTo(init.getCartesian());
+
+  //      assert Vector3D.angle(result, rotate.getCartesian()) == dTheta;
+        return Math.acos(result.getZ());
+
+    }
+
+    /**
+     * Симулируем пролет num электронов.
+     *
+     * @param E
+     * @param num
+     * @return
+     */
+    public ArrayList<SimulaionResult> simulateAll(double E, int num) {
+        ArrayList<SimulaionResult> res = new ArrayList();
+        double theta;
+
+        System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", E, num);
+        for (int i = 0; i < num; i++) {
+//            System.out.printf("Running electron number %d", i);
+            theta = this.getRandomTheta();
+            res.add(this.simulateOne(E, theta));
+        }
+        return res;
+    }
+
+    private double normalizeTheta(double theta) {
+        double res = theta;
+        if (theta < 0) {
+            res = -theta;
+        }
+        if (res >= 0 && res <= Math.PI / 2) {
+            return res;
+        } else if (res > Math.PI / 2 && res <= Math.PI) {
+            return Math.PI - res;
+        } else if (res > Math.PI && res <= Math.PI * 3 / 2) {
+            return res - Math.PI;
+        } else {
+            throw new IllegalStateException();
+        }
+    }
+
+    public double getRandomTheta() {
+        double x = generator.next();
+        return Math.acos(x);
+    }
+
+//    /**
+//     * Генерируем случайный угол таким образом, чтобы электрон заведомо был
+//     * траппинговый
+//     *
+//     * @return
+//     */
+//    public double getRandomTrappedTheta() {
+//        double res = 0;
+//        while (res < this.thetaTransport) {
+//            res = this.getRandomTheta();
+//
+//        }
+//        return res;
+//    }
+    public class SimulaionResult {
+
+        public SimulaionResult(EndState state, double E, double theta, int collisionNumber) {
+            this.state = state;
+            this.E = E;
+            this.theta = theta;
+            this.collisionNumber = collisionNumber;
+        }
+        public EndState state;
+        public double E;
+        public double theta;
+        public int collisionNumber;
+    }
+}

src/main/java/hep/dataforge/trapping/Scatter.java

@@ -0,0 +1,100 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package hep.dataforge.trapping;
+
+import com.sun.jna.Library;
+import com.sun.jna.Native;
+import com.sun.jna.ptr.DoubleByReference;
+
+/**
+ *
+ * @author Darksnake
+ */
+public class Scatter {
+
+    public interface LibScatter extends Library {
+
+        public void randomel(double E, DoubleByReference Eloss, DoubleByReference theta);
+
+        public void randomexc(double E, DoubleByReference Eloss, DoubleByReference theta);
+
+        public void randomion(double E, DoubleByReference Eloss, DoubleByReference theta);
+
+        public double sigmael(double E);
+
+        public double sigmaexc(double E);
+
+        public double sigmaion(double E);
+    }
+
+    /**
+     * PENDING переделать, чтобы возвращались нормальные значения
+     * @param E
+     * @param Eloss
+     * @param theta 
+     */
+    static void randomel(double E, DoubleValue Eloss, DoubleValue theta) {
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+
+        DoubleByReference ElossPointer = new DoubleByReference(Eloss.getValue());
+        DoubleByReference thetaPointer = new DoubleByReference(theta.getValue());
+        lib.randomel(E, ElossPointer, thetaPointer);
+        Eloss.setValue(ElossPointer.getValue());
+        theta.setValue(thetaPointer.getValue());
+
+    }
+
+    static void randomexc(double E, DoubleValue Eloss, DoubleValue theta) {
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+
+        DoubleByReference ElossPointer = new DoubleByReference(Eloss.getValue());
+        DoubleByReference thetaPointer = new DoubleByReference(theta.getValue());
+        lib.randomexc(E, ElossPointer, thetaPointer);
+        Eloss.setValue(ElossPointer.getValue());
+        theta.setValue(thetaPointer.getValue());
+
+    }
+
+    static void randomion(double E, DoubleValue Eloss, DoubleValue theta) {
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+
+        DoubleByReference ElossPointer = new DoubleByReference(Eloss.getValue());
+        DoubleByReference thetaPointer = new DoubleByReference(theta.getValue());
+        lib.randomion(E, ElossPointer, thetaPointer);
+        Eloss.setValue(ElossPointer.getValue());
+        theta.setValue(thetaPointer.getValue());
+
+    }
+
+    /**
+     * Все сечения в м^2
+     * @param E
+     * @return 
+     */
+    public static double sigmael(double E) {
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+        return lib.sigmael(E);
+    }
+
+    public static double sigmaexc(double E) {
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+        return lib.sigmaexc(E);
+    }
+
+    public static double sigmaion(double E) {
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+        return lib.sigmaion(E);
+    }
+    
+    /**
+     * Полное сечение с учетом квазиупругих столкновений
+     * @param E
+     * @return 
+     */
+    public static double sigmaTotal(double E){
+        LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class);
+        return lib.sigmael(E)+lib.sigmaexc(E)+lib.sigmaion(E);
+    }
+}

src/main/java/hep/dataforge/trapping/Trapping.java

@@ -0,0 +1,93 @@
+package hep.dataforge.trapping;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.Iterator;
+
+/**
+ * Hello world!
+ *
+ */
+public class Trapping {
+
+    public static void main(String[] args) throws FileNotFoundException {
+        PrintWriter out = null;
+        if ((args != null) && (args[0] != null)) {
+            File file = new File(args[0]);
+            out = new PrintWriter(file);
+        }
+
+        double E = 18000d;
+        System.out.println();
+//        System.setProperty("jna.library.path", "d:\\projects\\Trapping\\target\\classes\\win32-amd64\\libScatter.dll");
+        ElectronTrappingSimulator simulator = new ElectronTrappingSimulator();
+        simulator.setFields(0.6, 3.6, 7.2);
+//        ElectronTrappingSimulator.SimulaionResult result;
+
+        int accepted = 0;
+        int pass = 0;
+        int rejected = 0;
+        int lowE = 0;
+
+        ArrayList<ElectronTrappingSimulator.SimulaionResult> results = simulator.simulateAll(E, 100000);
+
+        for (Iterator<ElectronTrappingSimulator.SimulaionResult> it = results.iterator(); it.hasNext();) {
+            ElectronTrappingSimulator.SimulaionResult res = it.next();
+
+            if (out != null) {
+                out.printf("%g\t%g\t%d\t%s%n", res.E, res.theta * 180 / Math.PI, res.collisionNumber, res.state.toString());
+            }
+            switch (res.state) {
+
+                case ACCEPTED:
+                    System.out.printf("%g\t%g\t%d\t%s%n", res.E, res.theta * 180 / Math.PI, res.collisionNumber, res.state.toString());
+                    accepted++;
+                    break;
+                case LOWENERGY:
+                    lowE++;
+                    break;
+                case PASS:
+                    pass++;
+                    break;
+                case REJECTED:
+                    rejected++;
+                    break;
+                case NONE:
+                    throw new IllegalStateException();
+            }
+        }
+
+        System.out.printf("%nThe total number of events is %d.%n%n", results.size());
+
+        System.out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI);
+        System.out.printf("The transport mirroring angle is %g.%n", simulator.thetaTransport * 180 / Math.PI);
+        System.out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI);
+        System.out.printf("The starting energy is %g.%n", E);
+        System.out.printf("The lower energy boundary is %g.%n%n", simulator.Elow);
+
+        System.out.printf("The total number of ACCEPTED events is %d.%n", accepted);
+        System.out.printf("The total number of PASS events is %d.%n", pass);
+        System.out.printf("The total number of REJECTED events is %d.%n", rejected);
+        System.out.printf("The total number of LOWENERGY events is %d.%n%n", lowE);
+
+        if (out != null) {
+            out.printf("The total number of events is %d.%n%n", results.size());
+
+            out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI);
+            out.printf("The transport mirroring angle is %g.%n", simulator.thetaTransport * 180 / Math.PI);
+            out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI);
+            out.printf("The starting energy is %g.%n", E);
+            out.printf("The lower energy boundary is %g.%n%n", simulator.Elow);
+
+            out.printf("The total number of ACCEPTED events is %d.%n", accepted);
+            out.printf("The total number of PASS events is %d.%n", pass);
+            out.printf("The total number of REJECTED events is %d.%n", rejected);
+            out.printf("The total number of LOWENERGY events is %d.%n%n", lowE);
+            out.close();
+        }
+        
+
+    }
+}

src/main/java/hep/dataforge/trapping/TrappingRandomGenerator.java

@@ -0,0 +1,43 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package hep.dataforge.trapping;
+
+import org.apache.commons.math3.random.JDKRandomGenerator;
+import org.apache.commons.math3.random.RandomGenerator;
+
+/**
+ *
+ * @author Darksnake
+ */
+public class TrappingRandomGenerator {
+    RandomGenerator generator;
+
+    public TrappingRandomGenerator() {
+        this.generator = new JDKRandomGenerator();
+    }
+
+    public TrappingRandomGenerator(RandomGenerator generator) {
+        this.generator = generator;
+    }
+    
+    /**
+     * heads-tails random.
+     * @return 
+     */
+    public boolean heads(){
+        return generator.nextBoolean();
+    }
+    
+    /**
+     * next uniform in [0;1] 
+     * @return 
+     */
+    public double next(){
+        return generator.nextDouble();
+    }
+    
+    
+    
+}

src/test/java/hep/dataforge/trapping/AppTest.java

@@ -0,0 +1,38 @@
+package hep.dataforge.trapping;
+
+import junit.framework.Test;
+import junit.framework.TestCase;
+import junit.framework.TestSuite;
+
+/**
+ * Unit test for simple App.
+ */
+public class AppTest 
+    extends TestCase
+{
+    /**
+     * Create the test case
+     *
+     * @param testName name of the test case
+     */
+    public AppTest( String testName )
+    {
+        super( testName );
+    }
+
+    /**
+     * @return the suite of tests being tested
+     */
+    public static Test suite()
+    {
+        return new TestSuite( AppTest.class );
+    }
+
+    /**
+     * Rigourous Test :-)
+     */
+    public void testApp()
+    {
+        assertTrue( true );
+    }
+}