def test2(self):
     "neutron"
     # creating instances
     from mcni import neutron
     n = neutron(r=(0,0,0), v=(0,0,3000))
     n = neutron(r=(0,0,0), v=(0,0,3000), time=1000.)
     n = neutron(r=(0,0,0), v=(0,0,3000), time=1000., prob=10.)
     n = neutron(r=(0,0,0), v=(0,0,3000), s=(0,1), time=1000., prob=10.)
     # printing
     print n
     print n.state
     print n.state.velocity
     # setting values
     n.state.position = mcni.position(1,2,3)
     self.assertEqual(tuple(n.state.position), (1,2,3))
     n.state.velocity = mcni.velocity(1,2,3)
     self.assertEqual(tuple(n.state.velocity), (1,2,3))
     n.state.spin = mcni.spin(1,0)
     self.assertEqual(n.state.spin.s1, 1)
     self.assertEqual(n.state.spin.s2, 0)
     n.time = 999
     self.assertEqual(n.time, 999)
     n.probability = 888
     self.assertEqual(n.probability, 888)
     return

    def test_scatterercopy(self):
        '''scatterercopy'''
        # create a weird shape
        from mccomposite.geometry import primitives
        shape = primitives.block( (1,1,1) )

        #create pure python representation of scatterer composite
        composite1 = mccomposite.composite( shape )
        nprinter = NeutronPrinter( shape )
        composite1.addElement( nprinter )
        #create a copy
        copy = mccomposite.scatterercopy( composite1 )
        
        #create a larget composite
        shape = primitives.block( (1,1,2) )
        composite = mccomposite.composite( shape )
        composite.addElement( composite1, (0,0,-0.5) )
        composite.addElement( copy, (0,0,+0.5) )

        #render the c++ representation
        ccomposite = mccomposite.scattererEngine( composite )

        ev = mcni.neutron( r = (0,0,-5), v = (0,0,1) )
        ccomposite.scatter(ev)
        
        return

    def test1(self):
        'mccomponents.sample.samplecomponent isotropic kernel, multiple-scattering'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,3000), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'Al', 'sampleassemblies/Al-isotropickernel/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,1), (0,0,0) )

        N0 = 1
        neutrons = mcni.neutron_buffer(N0)

        mcni.simulate( instrument, geometer, neutrons, multiple_scattering=True)

        N = len(neutrons)

        for i in range(N):
            neutron = neutrons[i]
            print neutron
            continue

        return

    def test3(self):
        '''create pure python representation of a homogeneous scatterer with
        composite kernel. render the c++ computation engine of that kernel.
        '''
        #shape
        from mccomposite.geometry import primitives
        shape = primitives.block( (1,1,1) )

        #kernel
        nprinter = NeutronPrinter( )
        
        #composite kernel
        composite_kernel = hs.compositeKernel()
        composite_kernel.addElement( nprinter )

        #scatterer
        scatterer = hs.homogeneousScatterer(
            shape, composite_kernel)

        #render the c++ representation
        cscatterer = hs.scattererEngine( scatterer )

        for i in range(10):
            ev = mcni.neutron( r = (0,0,-5), v = (0,0,1) )
            cscatterer.scatter(ev)
            continue
        return

    def testCompositeScatteringKernel(self):
        'CompositeScatteringKernel'
        shape = mccompositebp.Block(1,1,1)

        from neutron_printer3 import cKernel as Printer
        printer = Printer( )

        kernels = mccomponentsbp.pointer_vector_Kernel(0)
        kernels.append( printer )
        
        weights = mccomponentsbp.vector_double(0)
        weights.append(1.)

        rotmats = mccomponentsbp.vector_rotmat(0)
        rotmat = mcnibp.RotationMatrix_double(1,0,0, 0,1,0, 0,0,1)
        rotmats.append(rotmat)
        
        average=False
        kernelcomposite = mccomponentsbp.CompositeScatteringKernel( 
            kernels, weights, rotmats, average)

        mcweights = mccomponentsbp.MCWeights_AbsorptionScatteringTransmission()
        scatterer = mccomponentsbp.HomogeneousNeutronScatterer(
            shape, kernelcomposite, mcweights )

        for i in range(10):
            ev = mcni.neutron( r = (0,0,-5), v = (0,0,1) )
            scatterer.scatter(ev)
            continue
        return

 def test1(self):
     'shape positioning: plate perp to beam'
     # source
     from mcni.components.MonochromaticSource import MonochromaticSource
     import mcni
     neutron = mcni.neutron(r=(0,0,-1), v=(0,0,1000), prob=1)
     source = MonochromaticSource('s', neutron, dx=0.07, dy=0.09, dE=0)
     # sample
     from mccomponents.sample import samplecomponent
     scatterer = samplecomponent('sa', 'plate/sampleassembly.xml' )
     # neutrons
     N = 1000
     neutrons = mcni.neutron_buffer(N)
     neutrons = source.process(neutrons)
     # find neutrons out of target
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     missing = (x>0.03) | (x<-0.03) | (y>0.04) | (y<-0.04)
     # print neutrons
     scatterer.process(neutrons)
     # print neutrons
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     assert (z[missing] < -.9).all()
     hit = arr[np.logical_not(missing), :3]
     x,y,z = hit.T
     assert (z > -.1).all()
     assert (np.isclose(np.abs(x), 0.03) | np.isclose(np.abs(y), 0.04) | np.isclose(np.abs(z), 0.005)).all()
     return

 def test2(self):
     'shape positioning: cylinder with axis along beam'
     # source
     from mcni.components.MonochromaticSource import MonochromaticSource
     import mcni
     neutron = mcni.neutron(r=(0,0,-1), v=(0,0,1000), prob=1)
     source = MonochromaticSource('s', neutron, dx=0.09, dy=0.09, dE=0)
     # sample
     from mccomponents.sample import samplecomponent
     scatterer = samplecomponent('sa', 'cyl-along-beam/sampleassembly.xml' )
     # neutrons
     N = 1000
     neutrons = mcni.neutron_buffer(N)
     neutrons = source.process(neutrons)
     # find neutrons out of target
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     missing = x*x+y*y > 0.04**2
     # print neutrons
     scatterer.process(neutrons)
     # print neutrons
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     assert (z[missing] < -.9).all()
     hit = arr[np.logical_not(missing), :3]
     x,y,z = hit.T
     assert (z > -.1).all()
     assert (np.isclose((x*x + y*y)**.5, 0.04) | np.isclose(np.abs(z), 0.005)).all()
     return

 def process(self, neutrons):
     import mcni
     for i in range(len(neutrons)):
         info.log("loop #%d" % i)
         neutrons[i] = mcni.neutron( r = ( 0,0,0 ), v = (0,0,random.random()) )
         continue
     return neutrons

    def process(self, neutrons):
        import mcni

        for i in range(len(neutrons)):
            neutrons[i] = mcni.neutron(r=(1, 2, 3), v=(1, 2, 3))
            continue
        return neutrons

    def test(self):
        'neutron_storage.Storage: write and then read'

        path = 'test-storage'
        if os.path.exists(path):
            os.remove( path )
        
        from mcni.neutron_storage.Storage import Storage

        #open storage for writing
        s = Storage( path, 'w' )

        #create neutrons
        import mcni
        neutrons = mcni.neutron_buffer( 7 )
        neutrons[5] = mcni.neutron( v = (8,9,10) )

        #write neutrons
        s.write( neutrons )

        #delete the storage to make sure it flushes all neutrons
        del s

        #open the storage for reading
        sr = Storage( path, 'r')
        neutrons = sr.read()
        self.assertEqual( len(neutrons), 7 )

        self.assertAlmostEqual( neutrons[5].state.velocity[0] , 8 )
        return

    def test3(self):
        'neutron_storage.Storage: wrap-reading (nread>ntotal)'

        path = 'test-storage-3'
        if os.path.exists(path):
            os.remove( path )
        
        from mcni.neutron_storage.Storage import Storage

        #open storage for writing
        s = Storage( path, 'w' )

        #create neutrons
        import mcni
        neutrons = mcni.neutron_buffer( 7 )
        for i in range(7):
            neutrons[i] = mcni.neutron( v = (i,0,0) )

        #write 
        s.write( neutrons )

        # flush
        del s
        
        #open the storage for reading
        sr = Storage( path, 'r', packetsize=10)

        neutrons = sr.read()
        self.assertEqual( len(neutrons), 10 )
        self.assertAlmostEqual( neutrons[5].state.velocity[0] , 5 )
        self.assertAlmostEqual( neutrons[9].state.velocity[0] , 2 )

        neutrons = sr.read()
        self.assertAlmostEqual( neutrons[0].state.velocity[0] , 3 )
        return

    def test5(self):
        'neutron_storage.Storage: wrap-reading (nread>>ntotal)'

        path = 'test-storage-4'
        if os.path.exists(path):
            os.remove( path )
        
        from mcni.neutron_storage.Storage import Storage

        #open storage for writing
        s = Storage( path, 'w' )

        #create neutrons
        import mcni
        neutrons = mcni.neutron_buffer( 5 )
        for i in range(5):
            neutrons[i] = mcni.neutron( v = (i,0,0) )

        #write 
        s.write( neutrons )

        # flush
        del s
        
        #open the storage for reading
        sr = Storage( path, 'r')

        neutrons = sr.read(100)
        self.assertEqual( len(neutrons), 100 )
        self.assertAlmostEqual( neutrons[3].state.velocity[0] , 3 )
        self.assertAlmostEqual( neutrons[4].state.velocity[0] , 4 )
        self.assertAlmostEqual( neutrons[6].state.velocity[0] , 1 )
        self.assertAlmostEqual( neutrons[7].state.velocity[0] , 2 )

        return

    def test1(self):
        'complex. pack, detector, pixel hierarchy'
        
        mca = md.eventModeMCA(
            outfilename,
            (npacks, ndetsperpack, npixelsperdet,) )
        cylinder = operations.subtract( primitives.cylinder( sample2det * 1.1, detlength ),
                                        primitives.cylinder( sample2det * 0.9, detlength ) )
        ds = md.detectorSystem( cylinder, tofparams, mca )

        pack0 = makepack()
        packs = [pack0]
        for i in range(1, npacks):
            packs.append( mccomposite.scatterercopy( pack0, id = i ) )
            continue

        for i in range( npacks ):
            z = 0 * meter
            angle = (i-packindexat0)* 5./180 * N.pi
            x = sample2det * math.cos(angle)
            y = sample2det * math.sin(angle)
            ds.addElement( packs[i], (x,y,z) )
            continue

        cds = mh.scattererEngine( ds, coordinate_system = "InstrumentScientist" )

        for i in range(nevents):
            if i%1000 == 0: print i
            ev = mcni.neutron( r = (-L1,0,0), v = (vi,0,0) )
            cds.scatter(ev)
            continue

        return

    def test(self):
        "wrap SNS_source_r1"
        from mcstas2 import componentfactory

        factory = componentfactory(category, componentname)
        component = factory(
            "component",
            S_filename="source_sct521_bu_17_1.dat",
            width=0.1,
            height=0.12,
            dist=2.5,
            xw=0.1,
            yh=0.12,
            Emin=50,
            Emax=70,
        )

        import mcni

        neutrons = mcni.neutron_buffer(5)
        for i in range(5):
            neutrons[i] = mcni.neutron(r=(0, 0, -1), v=(0, 0, 3000), time=0, prob=1)
            continue
        component.process(neutrons)
        print neutrons
        return

    def test(self):
        E_Q = "Q*Q/3."
        S_Q = "1"
        Qmin = 0; Qmax = 10
        absorption_coefficient = scattering_coefficient = 1.
        kernel = mccomponentsbp.create_E_Q_Kernel(
            E_Q, S_Q,
            Qmin, Qmax,
            absorption_coefficient,
            scattering_coefficient,
            )

        ei = 500 # meV
        from mcni.utils import conversion
        vil = conversion.e2v(ei)
        vi = (0,0,vil)

        import numpy.linalg as nl
        import numpy as np
        for i in range(10):
            event = mcni.neutron( 
                r = (0,0,0), v = vi, 
                prob = 1, time = 0 )
            kernel.scatter( event );
            vf = np.array(event.state.velocity)
            diffv = vi - vf
            Q = conversion.v2k(nl.norm(diffv))
            ef = conversion.v2e(nl.norm(vf))
            E = ei - ef
            # print E, Q, event
            E1 = eval(E_Q)
            self.assertAlmostEqual(E, E1)
            continue

        return

    def test1(self):
        'detector hierarchy from xml'
        from instrument.nixml import parse_file
        instrument = parse_file( 'ARCS.xml' )

        import instrument.geometers as ig
        instrument.geometer.changeRequestCoordinateSystem(
            ig.coordinateSystem( coordinate_system ) )
        
        assignLocalGeometers( instrument, coordinate_system = coordinate_system )
        
        detectorSystem = instrument.getDetectorSystem()

        tofparams = 0, 10e-3, 1e-4
        detectorSystem.tofparams = tofparams
        dims = getDetectorHierarchyDimensions( instrument )
        dims = [ dim for name, dim in dims ]
        mca = md.eventModeMCA(  outfilename, dims )
        detectorSystem.mca = mca
        
        cds = mh.scattererEngine( detectorSystem, coordinate_system = coordinate_system )

        for i in range(nevents):
            if i%1000 == 0: print i
            ev = mcni.neutron( r = (0,0,0), v = (1500,0,2000) )
            cds.scatter(ev)
            continue

        instrument.geometer = instrument.global_geometer
        return

    def test1(self):
        'mccomponents.sample.samplecomponent'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,4149.48), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'V-constantE', 'sampleassemblies/V-constantE/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,1), (0,0,0) )
        
        neutrons = mcni.neutron_buffer( 1 )
        
        from mcni.pyre_support.ConsoleNeutronTracer import ConsoleNeutronTracer
        tracer = ConsoleNeutronTracer()
        mcni.simulate( 
            instrument, geometer, neutrons, 
            multiple_scattering=True,
            tracer = tracer
            )

        return

    def test1(self):
        'mccomponents.sample.samplecomponent: DGSSXResKernel'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,3000), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'Al', 'sampleassemblies/Al-DGSSXResKernel/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,6), (0,0,0) )

        N0 = 10
        neutrons = mcni.neutron_buffer(N0)

        mcni.simulate( instrument, geometer, neutrons )

        N = len(neutrons)

        for i in range(10):
            neutron = neutrons[i]
            # print neutron
            self.assert_(np.allclose(neutron.state.velocity, [3000, 0, 0], atol=20))
            continue

        return

    def test_copy(self):
        '''copy'''
        print "This test creates two identical blocks, each of which "\
              "does not interact with neutrons. They print the info "\
              "about the neutrons passing thru them, however. "\
              "This test then send one neutron through these two "\
              "blocks, so we should see two printings of neutron info, "\
              "differing only on time-of-flight."
        # create a shape
        from mccomposite.geometry import primitives
        smallblock = primitives.block( (1,1,1) )

        #create pure python representation of scatterer composite
        composite1 = mccomposite.composite( smallblock )
        import UseNeutronPrinter2
        nprinter = UseNeutronPrinter2.NeutronPrinter( smallblock )
        composite1.addElement( nprinter )
        #create a copy
        copy = Copy.Copy( composite1 )
        
        #create a larget composite
        largeblock = primitives.block( (1,1,2) )
        composite = mccomposite.composite( largeblock )
        composite.addElement( composite1, (0,0,-0.5) )
        #composite.addElement( nprinter, (0,0,-0.5) )
        composite.addElement( copy, (0,0,+0.5) )

        #render the c++ representation
        ccomposite = mccomposite.scattererEngine( composite )

        ev = mcni.neutron( r = (0,0,-5), v = (0,0,1) )

        ccomposite.scatter(ev)
        return

    def test1(self):
        'mccomponents.sample.samplecomponent: IsotropicKernel'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,3000), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'Al', 'sampleassemblies/Al-isotropickernel/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,1), (0,0,0) )

        N0 = 1000
        neutrons = mcni.neutron_buffer(N0)

        mcni.simulate( instrument, geometer, neutrons )

        N = len(neutrons)

        for i in range(10):
            neutron = neutrons[i]
            print neutron
            continue

        # N should be about 2/3 of N0. this is determined by
        # the mc weights in Al-isotropic-kernel-plate-scatterer.xml
        self.assert_( N < 0.72*N0 and N > 0.6*N0)

        return