contrib/haskell/src/Hkl/Xrd/OneD.hs

   1 {-# LANGUAGE CPP #-}
   2 {-# LANGUAGE FlexibleInstances #-}
   3 {-# LANGUAGE GADTs #-}
   4 {-# LANGUAGE MultiParamTypeClasses #-}
   5 {-# LANGUAGE OverloadedStrings #-}
   6 {-# LANGUAGE Rank2Types #-}
   7 {-# LANGUAGE StandaloneDeriving #-}
   8 {-# LANGUAGE UnicodeSyntax #-}
   9
  10 module Hkl.Xrd.OneD
  11        ( XrdOneD
  12        , XRDRef(..)
  13        , XrdRefSource(..)
  14        , XRDSample(..)
  15        , Threshold(..)
  16        , XrdNxs(..)
  17        , XrdOneDParams(..)
  18        , XrdSource(..)
  19        , PoniExt(..)
  20          -- reference
  21        , getPoseEdf
  22        , getPoniExtRef
  23          -- integration
  24        , integrate
  25        , substract
  26          -- integrateMulti
  27        , integrateMulti
  28        , substractMulti
  29          -- tools
  30        , dummiesForPy
  31        ) where
  32
  33 import Control.Applicative ((<$>), (<*>), pure)
  34 import Control.Concurrent.Async (mapConcurrently)
  35 import Control.Monad (forM_, forever, void, when, zipWithM_)
  36 import Control.Monad.Morph (hoist)
  37 import Control.Monad.Trans.Maybe (MaybeT, runMaybeT)
  38 import Control.Monad.Trans.State.Strict (StateT, get, put)
  39 import Data.Array.Repa (DIM1, ix1, size)
  40 import Data.Attoparsec.Text (parseOnly)
  41 import qualified Data.List as List (lookup)
  42 import Data.Maybe (fromJust, fromMaybe, isJust)
  43 import Data.Text (Text)
  44 import qualified Data.Text as Text (unlines, pack, intercalate)
  45 import Data.Text.IO (readFile)
  46 import Data.Vector.Storable (concat, head)
  47 import Numeric.LinearAlgebra (fromList)
  48 import Numeric.Units.Dimensional.Prelude (meter, nano, (/~), (*~))
  49 import System.Exit ( ExitCode( ExitSuccess ) )
  50 import System.FilePath ((<.>), (</>), dropExtension, replaceExtension, takeFileName, takeDirectory)
  51 import Text.Printf ( printf )
  52
  53 import Pipes
  54     ( Consumer
  55     , Pipe
  56     , lift
  57     , (>->)
  58     , runEffect
  59     , await
  60     , yield
  61     )
  62 import Pipes.Lift ( evalStateP )
  63 import Pipes.Prelude ( drain, filter, toListM )
  64 import Pipes.Safe ( runSafeT )
  65
  66 import Hkl.C ( Factory ( K6c )
  67              , Geometry ( Geometry )
  68              , geometryDetectorRotationGet
  69              )
  70 import Hkl.DataSource ( DataItem ( DataItemH5 )
  71                       , DataSource( DataSourceH5 )
  72                       , atIndex'
  73                       )
  74 import Hkl.Detector ( Detector ( ZeroD ) )
  75 import Hkl.Edf ( Edf ( Edf )
  76                , edfFromFile
  77                )
  78 import Hkl.Flat ( Flat )
  79 import Hkl.H5 ( lenH5Dataspace )
  80 import Hkl.PyFAI ( AIMethod, Poni
  81                  , PoniExt ( PoniExt )
  82                  , PoniPath
  83                  , Pose ( Pose )
  84                  , move
  85                  , poniP
  86                  , poniToText
  87                  )
  88 import Hkl.Python ( PyVal
  89                   , toPyVal
  90                   )
  91 import Hkl.MyMatrix ( Basis ( HklB )
  92                     , MyMatrix ( MyMatrix )
  93                     )
  94 import Hkl.Nxs ( DataFrameH5 ( DataFrameH5 )
  95                , Nxs ( Nxs )
  96                , XrdOneD
  97                , DataFrameH5Path ( XrdOneDH5Path )
  98                , withDataFrameH5
  99                )
 100 import Hkl.Script ( Gnuplot, Py2
 101                   , Script ( ScriptGnuplot, Py2Script )
 102                   , run
 103                   , scriptSave
 104                   )
 105 import Hkl.Types ( AbsDirPath, SampleName
 106                  , Source ( Source )
 107                  )
 108 import Hkl.Utils ( hasContent )
 109
 110 -- | TODO
 111 -- * When we skip the last frame there is problem.
 112
 113 -- * Let's add a method in order to customize the movement of the poni.
 114
 115 -- | Types
 116
 117 data Threshold = Threshold Int deriving (Show)
 118
 119 instance PyVal Threshold where
 120   toPyVal (Threshold i) = toPyVal i
 121
 122 data XrdRefSource = XrdRefNxs (Nxs XrdOneD) Int
 123                   | XrdRefEdf FilePath PoniPath
 124                   deriving (Show)
 125
 126 data XRDRef = XRDRef SampleName AbsDirPath XrdRefSource
 127             deriving (Show)
 128
 129 data XrdSource = XrdSourceNxs (Nxs XrdOneD)
 130                | XrdSourceEdf [FilePath]
 131                  deriving (Show)
 132
 133 data XrdNxs
 134     = XrdNxs
 135       DIM1 -- bins
 136       DIM1 -- bins for the multibins
 137       (Maybe Threshold) -- threshold use to remove image Intensity
 138       [Int] -- Index of the frames to skip
 139       XrdSource -- data source
 140     deriving (Show)
 141
 142 data XRDSample = XRDSample SampleName AbsDirPath [XrdNxs] -- ^ nxss
 143                deriving (Show)
 144
 145 data XrdOneDParams a = XrdOneDParams PoniExt (Maybe (Flat a)) AIMethod
 146
 147 data DifTomoFrame sh =
 148   DifTomoFrame { difTomoFrameNxs :: Nxs XrdOneD-- ^ nexus of the current frame
 149                , difTomoFrameIdx :: Int -- ^ index of the current frame
 150                , difTomoFrameEOF :: Bool -- ^ is it the eof of the stream
 151                , difTomoFrameGeometry :: Geometry -- ^ diffractometer geometry
 152                , difTomoFramePoniExt :: PoniExt -- ^ the ref poniext
 153                } deriving (Show)
 154
 155 class Frame t where
 156   len :: t -> IO (Maybe Int)
 157   row :: t -> Int -> MaybeT IO (DifTomoFrame DIM1)
 158
 159 instance Frame (DataFrameH5 XrdOneD) where
 160   len (DataFrameH5 _ _ _ (DataSourceH5 _ d) _ _) = lenH5Dataspace d
 161
 162   row d@(DataFrameH5 nxs' _ g d' w ponigen) idx = do
 163     n <- lift $ len d
 164     let eof = fromJust n - 1 == idx
 165     let mu = 0.0
 166     let komega = 0.0
 167     let kappa = 0.0
 168     let kphi = 0.0
 169     gamma <- g `atIndex'` (ix1 0)
 170     delta <- d' `atIndex'` (ix1 idx)
 171     wavelength <- w `atIndex'` (ix1 0)
 172     let source = Source (Data.Vector.Storable.head wavelength *~ nano meter)
 173     let positions = Data.Vector.Storable.concat [mu, komega, kappa, kphi, gamma, delta]
 174     -- print positions
 175     let geometry =  Geometry K6c source positions Nothing
 176     let detector = ZeroD
 177     m <- lift $ geometryDetectorRotationGet geometry detector
 178     let pose = Pose (MyMatrix HklB m)
 179     poniext <- lift $ ponigen pose idx
 180     return $ DifTomoFrame { difTomoFrameNxs = nxs'
 181                           , difTomoFrameIdx = idx
 182                           , difTomoFrameEOF = eof
 183                           , difTomoFrameGeometry = geometry
 184                           , difTomoFramePoniExt = poniext
 185                           }
 186
 187 -- type PipeE e a b m r = EitherT e (Pipe a b m) r
 188
 189 frames :: (Frame a) => Pipe a (DifTomoFrame DIM1) IO ()
 190 frames = do
 191   d <- await
 192   (Just n) <- lift $ len d
 193   forM_ [0..n-1] (\i' -> do
 194                      f <- lift $ runMaybeT $ row d i'
 195                      when (isJust f) (yield (fromJust f)))
 196
 197 frames' :: (Frame a) => [Int] -> Pipe a (DifTomoFrame DIM1) IO ()
 198 frames' is = do
 199   d <- await
 200   forM_ is (\i' -> do
 201               f <- lift $ runMaybeT $ row d i'
 202               when (isJust f) (yield (fromJust f)))
 203
 204 skip :: [Int] -> DifTomoFrame sh -> Bool
 205 skip is' (DifTomoFrame _ i _ _ _) = notElem i is'
 206
 207 -- {-# ANN module "HLint: ignore Use camelCase" #-}
 208
 209
 210 -- import Graphics.Rendering.Chart.Easy
 211 -- import Graphics.Rendering.Chart.Backend.Diagrams
 212
 213 -- plotPonies :: FilePath -> [PoniEntry] -> IO ()
 214 -- plotPonies f entries = toFile def f $ do
 215 --     layout_title .= "Ponies"
 216 --     setColors [opaque blue]
 217 --     let values = map extract entries
 218 --     plot (line "am" [values [0,(0.5)..400]])
 219 --     -- plot (points "am points" (signal [0,7..400]))
 220 --     where
 221 --       extract (PoniEntry _ _ (Length poni1) _ _ _ _ _ _) = poni1
 222
 223 -- | Usual methods
 224
 225 dummiesForPy ∷ Maybe Threshold → String
 226 dummiesForPy mt = unlines [ "# Compute the dummy values for the dynamic mask"
 227                           , "DUMMY=" ++ dummy
 228                           , "DELTA_DUMMY=" ++ delta_dummy
 229                           ]
 230   where
 231     dummy = maybe "None" (\_ → "4294967296") mt -- TODO the default value depends on the number od bits per pixels.
 232     delta_dummy = maybe "None" (\(Threshold t) → show (4294967296 - t)) mt
 233
 234 getScanDir ∷ AbsDirPath → FilePath → FilePath
 235 getScanDir o f = o </> (dropExtension . takeFileName) f
 236
 237 pgen :: AbsDirPath -> FilePath -> Int -> FilePath
 238 pgen o f i = o </> scandir </>  scandir ++ printf "_%02d.poni" i
 239   where
 240     scandir = (dropExtension . takeFileName) f
 241
 242 getPoseEdf :: FilePath -> IO Pose
 243 getPoseEdf f = do
 244   edf@(Edf lambda _) <- edfFromFile f
 245   let mnes = map Text.pack ["_mu", "_keta", "_kap", "_kphi", "nu", "del"]
 246   let source = Source lambda
 247   let positions = fromList $ map (extract edf) mnes
 248   let geometry =  Geometry K6c source positions Nothing
 249   let detector = ZeroD
 250   m <- geometryDetectorRotationGet geometry detector
 251   return $ Pose (MyMatrix HklB m)
 252     where
 253       extract :: Edf -> Text -> Double
 254       extract (Edf _ ms) key = fromMaybe 0.0 (List.lookup key ms)
 255
 256 poniFromFile :: FilePath -> IO Poni
 257 poniFromFile filename = do
 258   content <- Data.Text.IO.readFile filename
 259   return $ case parseOnly poniP content of
 260     Left _     -> error $ "Can not parse the " ++ filename ++ " poni file"
 261     Right poni -> poni
 262
 263 getPoniExtRef :: XRDRef -> IO PoniExt
 264 getPoniExtRef (XRDRef _ output (XrdRefNxs nxs'@(Nxs f _) idx)) = do
 265   poniExtRefs <- runSafeT $
 266                  toListM ( withDataFrameH5 nxs' (gen output f) yield
 267                            >-> hoist lift ( frames' [idx]))
 268   return $ difTomoFramePoniExt (Prelude.last poniExtRefs)
 269   where
 270     gen :: FilePath -> FilePath -> Pose -> Int -> IO PoniExt
 271     gen root nxs'' p idx' = PoniExt
 272                             <$> poniFromFile (root </> scandir ++ printf "_%02d.poni" idx')
 273                             <*> pure p
 274       where
 275         scandir = takeFileName nxs''
 276 getPoniExtRef (XRDRef _ _ (XrdRefEdf e p)) = PoniExt
 277                                              <$> poniFromFile p
 278                                              <*> getPoseEdf e
 279
 280 integrate ∷ XrdOneDParams a → [XRDSample] → IO ()
 281 integrate p ss = void $ mapConcurrently (integrate' p) ss
 282
 283 integrate' ∷ XrdOneDParams a → XRDSample → IO ()
 284 integrate' p (XRDSample _ output nxss) = void $ mapConcurrently (integrate'' p output) nxss
 285
 286 integrate'' ∷ XrdOneDParams a → AbsDirPath → XrdNxs → IO ()
 287 integrate'' p output (XrdNxs b _ mt is (XrdSourceNxs nxs'@(Nxs f _))) = do
 288   print f
 289   runSafeT $ runEffect $
 290     withDataFrameH5 nxs' (gen p) yield
 291     >-> hoist lift (frames
 292                     >-> Pipes.Prelude.filter (skip is)
 293                     >-> savePonies (pgen output f)
 294                     >-> savePy p b mt
 295                     >-> saveGnuplot
 296                     >-> drain)
 297   where
 298     gen :: XrdOneDParams a -> Pose -> Int -> IO PoniExt
 299     gen (XrdOneDParams ref' _ _) m _idx = return $ move ref' m
 300
 301 createPy ∷ XrdOneDParams a → DIM1 → Maybe Threshold → FilePath → DifTomoFrame' sh → (Script Py2, FilePath)
 302 createPy (XrdOneDParams _ mflat m) b mt scriptPath (DifTomoFrame' f poniPath) = (Py2Script (script, scriptPath), output)
 303     where
 304       script = Text.unlines $
 305                map Text.pack ["#!/bin/env python"
 306                              , ""
 307                              , "import numpy"
 308                              , "from h5py import File"
 309                              , "from pyFAI import load"
 310                              , ""
 311                              , "PONIFILE = " ++ toPyVal poniPath
 312                              , "NEXUSFILE = " ++ toPyVal nxs'
 313                              , "IMAGEPATH = " ++ toPyVal i'
 314                              , "IDX = " ++ toPyVal idx
 315                              , "N = " ++ toPyVal (size b)
 316                              , "OUTPUT = " ++ toPyVal output
 317                              , "WAVELENGTH = " ++ toPyVal (w /~ meter)
 318                              , ""
 319                              , "# load the flat"
 320                              , "flat = " ++ toPyVal mflat
 321                              , ""
 322                              , dummiesForPy mt
 323                              , ""
 324                              , "ai = load(PONIFILE)"
 325                              , "ai.wavelength = WAVELENGTH"
 326                              , "ai._empty = numpy.nan"
 327                              , ""
 328                              , "with File(NEXUSFILE, mode='r') as f:"
 329                              , "    img = f[IMAGEPATH][IDX]"
 330                              , ""
 331                              , "    # Compute the mask"
 332                              , "    mask = numpy.zeros_like(img, dtype=bool)"
 333                              , "    mask[:,550:] = True"
 334                              , "    #mask_module[0:50, :] = True"
 335                              , "    #mask_module[910:960, :] = True"
 336                              , "    #mask_module[:,0:10] = True"
 337                              , "    if flat is not None:  # this should be removed for pyFAI >= 0.13.1 it is now done by PyFAI"
 338                              , "        mask = numpy.logical_or(mask, flat == 0.0)"
 339                              , ""
 340                              , "    ai.integrate1d(img, N, filename=OUTPUT, unit=\"2th_deg\", error_model=\"poisson\", correctSolidAngle=False, method=\"" ++ show m ++ "\", mask=mask, flat=flat, dummy=DUMMY, delta_dummy=DELTA_DUMMY)"
 341                              ]
 342       (Nxs nxs' (XrdOneDH5Path (DataItemH5 i' _) _ _ _)) = difTomoFrameNxs f
 343       idx = difTomoFrameIdx f
 344       output = poniPath `replaceExtension` "dat"
 345       (Geometry _ (Source w) _ _) = difTomoFrameGeometry f
 346
 347 -- | Pipes
 348
 349 data DifTomoFrame' sh = DifTomoFrame' { difTomoFrame'DifTomoFrame :: DifTomoFrame sh
 350                                       , difTomoFrame'PoniPath :: FilePath
 351                                       }
 352
 353 savePonies :: (Int -> FilePath) -> Pipe (DifTomoFrame sh) (DifTomoFrame' sh) IO ()
 354 savePonies g = forever $ do
 355   f <- await
 356   let filename = g (difTomoFrameIdx f)
 357   let (PoniExt p _) = difTomoFramePoniExt f
 358   lift $ filename `hasContent` (poniToText p)
 359   yield $ DifTomoFrame' { difTomoFrame'DifTomoFrame = f
 360                         , difTomoFrame'PoniPath = filename
 361                         }
 362
 363 data DifTomoFrame'' sh = DifTomoFrame'' { difTomoFrame''DifTomoFrame' :: DifTomoFrame' sh
 364                                         , difTomoFrame''PySCript :: Script Py2
 365                                         , difTomoFrame''DataPath :: FilePath
 366                                         }
 367
 368 savePy ∷ XrdOneDParams a → DIM1 → Maybe Threshold → Pipe (DifTomoFrame' sh) (DifTomoFrame'' sh) IO ()
 369 savePy p b mt = forever $ do
 370   f@(DifTomoFrame' _difTomoFrame poniPath) <- await
 371   let scriptPath = poniPath `replaceExtension`"py"
 372   let (script, dataPath) = createPy p b mt scriptPath f
 373   ExitSuccess <- lift $ run script True
 374   yield $ DifTomoFrame'' { difTomoFrame''DifTomoFrame' = f
 375                          , difTomoFrame''PySCript = script
 376                          , difTomoFrame''DataPath = dataPath
 377                          }
 378
 379 data DifTomoFrame''' sh = DifTomoFrame''' { difTomoFrame'''DifTomoFrame'' ∷ DifTomoFrame'' sh
 380                                           , difTomoFrame'''GnuplotScript ∷ Script Gnuplot
 381                                           , difTomoFrame'''Curves ∷ [FilePath]
 382                                           }
 383
 384 mkGnuplot ∷ [FilePath] → FilePath → Script Gnuplot
 385 mkGnuplot fs o = ScriptGnuplot (content, o)
 386   where
 387     content = Text.unlines $
 388               ["plot \\"]
 389               ++ [Text.intercalate ",\\\n" [ Text.pack (show f ++ " u 1:2 w l") | f <- fs ]]
 390               ++ ["pause -1"]
 391
 392 saveGnuplot' :: Pipe (DifTomoFrame'' sh) (DifTomoFrame''' sh) (StateT [FilePath] IO) r
 393 saveGnuplot' = forever $ do
 394   curves <- lift get
 395   f@(DifTomoFrame'' (DifTomoFrame' _ poniPath) _ dataPath) <- await
 396   let curves' = curves ++ [dataPath]
 397   let script = mkGnuplot curves' (takeDirectory poniPath </> "plot.gnuplot")
 398   lift . lift $ scriptSave script
 399   lift $ put $! curves'
 400   yield $ DifTomoFrame''' { difTomoFrame'''DifTomoFrame'' = f
 401                           , difTomoFrame'''GnuplotScript = script
 402                           , difTomoFrame'''Curves = curves'
 403                           }
 404
 405 saveGnuplot :: Pipe (DifTomoFrame'' sh) (DifTomoFrame''' sh) IO r
 406 saveGnuplot = evalStateP [] saveGnuplot'
 407
 408 -- substract a sample from another one
 409
 410 substractPy ∷ [FilePath] → [FilePath] → [FilePath] → FilePath → Script Py2
 411 substractPy fs1 fs2 os scriptPath = Py2Script (content, scriptPath)
 412   where
 413     content ∷ Text
 414     content = Text.unlines $
 415               map Text.pack ["#!/bin/env python"
 416                             , ""
 417                             , "import numpy"
 418                             , ""
 419                             , "S1 = " ++ toPyVal fs1
 420                             , "S2 = " ++ toPyVal fs2
 421                             , "OUTPUTS = " ++ toPyVal os
 422                             , ""
 423                             , "def substract(f1, f2, o):"
 424                             , "    a1 = numpy.genfromtxt(f1)"
 425                             , "    a2 = numpy.genfromtxt(f2)"
 426                             , "    res = numpy.copy(a2)"
 427                             , "    res[:,1] -= a1[:,1]"
 428                             , "    # TODO deal with the error propagation"
 429                             , "    numpy.savetxt(output, res)"
 430                             , ""
 431                             , "for (s1, s2, output) in zip(S1, S2, OUTPUTS):"
 432                             , "    substract(s1, s2, output)"
 433                             ]
 434
 435 targetP ∷ (Int → FilePath) → Pipe (DifTomoFrame sh) FilePath IO ()
 436 targetP g = forever $ do
 437   f ← await
 438   let poniPath = g (difTomoFrameIdx f)
 439   let dataPath = poniPath `replaceExtension` "dat"
 440   yield dataPath
 441
 442 target' ∷  XrdOneDParams a → AbsDirPath → XrdNxs → IO (FilePath, [FilePath])
 443 target' p output (XrdNxs _ _ _ is (XrdSourceNxs nxs'@(Nxs f _))) = do
 444   fs ← runSafeT $ toListM $
 445        withDataFrameH5 nxs' (gen p) yield
 446        >-> hoist lift (frames
 447                        >-> Pipes.Prelude.filter (skip is)
 448                        >-> targetP (pgen output f)
 449                       )
 450   return (getScanDir output f, fs)
 451   where
 452     gen :: XrdOneDParams a -> Pose -> Int -> IO PoniExt
 453     gen (XrdOneDParams ref' _ _) m _idx = return $ move ref' m
 454
 455 targets ∷ XrdOneDParams a → XRDSample → IO [(FilePath, [FilePath])]
 456 targets p (XRDSample _ output nxss) = mapConcurrently (target' p output) nxss
 457
 458 substract' ∷ XrdOneDParams a → XRDSample → XRDSample → IO ()
 459 substract' p s1@(XRDSample name _ _) s2 = do
 460   -- compute the output of the s1 sample
 461   -- we take only the first list of the sample
 462   f1s:_ ← targets p s1
 463   -- compute the output of the s2 sample
 464   f2s ← targets p s2
 465   -- do the substraction via a python script and add the gnuplot file
 466   _ ← mapConcurrently (go f1s) f2s
 467   return ()
 468   where
 469     go ∷ (FilePath, [FilePath]) → (FilePath, [FilePath]) → IO ()
 470     go (_, f1) (d, f2) = do
 471       -- compute the substracted output file names take into account
 472       -- that f1 and f2 could have different length
 473       let outputs = [dropExtension f ++ "-" ++ name <.> "dat" | (_, f) ← zip f1 f2]
 474       -- compute the script name
 475       let scriptPath = d </> "substract.py"
 476       let script = substractPy f1 f2 outputs scriptPath
 477       ExitSuccess ← run script False
 478       -- gnuplot
 479       let gnuplotPath = d </> "substract.gnuplot"
 480       scriptSave $ mkGnuplot outputs gnuplotPath
 481       return ()
 482
 483 substract ∷  XrdOneDParams a → XRDSample → [XRDSample] → IO ()
 484 substract p s ss = mapM_ (substract' p s) ss
 485
 486 -- | PyFAI MultiGeometry
 487
 488 integrateMulti ∷ XrdOneDParams a → [XRDSample] → IO ()
 489 integrateMulti p samples = mapM_ (integrateMulti' p) samples
 490
 491 integrateMulti' ∷ XrdOneDParams a → XRDSample → IO ()
 492 integrateMulti' p (XRDSample _ output nxss) = mapM_ (integrateMulti'' p output) nxss
 493
 494 integrateMulti'' ∷ XrdOneDParams a → AbsDirPath → XrdNxs → IO ()
 495 integrateMulti'' p output (XrdNxs _ mb mt is (XrdSourceNxs nxs'@(Nxs f _))) = do
 496   print f
 497   runSafeT $ runEffect $
 498     withDataFrameH5 nxs' (gen p) yield
 499     >-> hoist lift (frames
 500                     >-> Pipes.Prelude.filter (skip is)
 501                     >-> savePonies (pgen output f)
 502                     >-> saveMultiGeometry p mb mt)
 503   where
 504     gen :: XrdOneDParams a -> Pose -> Int -> IO PoniExt
 505     gen (XrdOneDParams ref' _ _)  m _idx = return $ move ref' m
 506
 507 integrateMulti'' p output (XrdNxs b _ mt _ (XrdSourceEdf fs)) = do
 508   -- generate all the ponies
 509   zipWithM_ (go p) fs ponies
 510
 511   -- generate the multi.py python script
 512   let scriptPath = output </> "multi.py"
 513   let (script, _) = createMultiPyEdf p b mt fs ponies scriptPath (output </> "multi.dat")
 514   scriptSave script
 515     where
 516       ponies = [output </> (dropExtension . takeFileName) f ++ ".poni" | f <- fs]
 517
 518       go ∷ XrdOneDParams a → FilePath → FilePath → IO ()
 519       go (XrdOneDParams ref _ _) f o = do
 520         m <- getPoseEdf f
 521         let (PoniExt p' _) = move ref m
 522         o `hasContent` (poniToText p')
 523
 524 createMultiPy ∷ XrdOneDParams a → DIM1 → Maybe Threshold → FilePath → DifTomoFrame' sh → [(Int, FilePath)] → (Script Py2, FilePath)
 525 createMultiPy (XrdOneDParams _ mflat _) b mt scriptPath (DifTomoFrame' f _) idxPonies = (Py2Script (content, scriptPath), output)
 526     where
 527       content = Text.unlines $
 528                map Text.pack ["#!/bin/env python"
 529                              , ""
 530                              , "import numpy"
 531                              , "from h5py import File"
 532                              , "from pyFAI.multi_geometry import MultiGeometry"
 533                              , ""
 534                              , "NEXUSFILE = " ++ toPyVal nxs'
 535                              , "IMAGEPATH = " ++ toPyVal i'
 536                              , "BINS = " ++ toPyVal (size b)
 537                              , "OUTPUT = " ++ toPyVal output
 538                              , "WAVELENGTH = " ++ toPyVal (w /~ meter)
 539                              , "THRESHOLD = " ++ toPyVal mt
 540                              , ""
 541                              , "# load the flat"
 542                              , "flat = " ++ toPyVal mflat
 543                              , ""
 544                              , "# Load all images"
 545                              , "PONIES = " ++ toPyVal ponies
 546                              , "IDXS = " ++ toPyVal idxs
 547                              , ""
 548                              , "# Read all the images"
 549                              , "imgs = []"
 550                              , "with File(NEXUSFILE, mode='r') as f:"
 551                              , "    for idx in IDXS:"
 552                              , "        imgs.append(f[IMAGEPATH][idx])"
 553                              , ""
 554                              , "# Compute the mask"
 555                              , "mask = numpy.zeros_like(imgs[0], dtype=bool)"
 556                              , "mask[:,550:] = True"
 557                              , "if flat is not None:  # this should be removed for pyFAI >= 0.13.1 it is now done by PyFAI"
 558                              , "    mask = numpy.logical_or(mask, flat == 0.0)"
 559                              , "lst_mask = []"
 560                              , "for img in imgs:  # remove all pixels above the threshold"
 561                              , "    if THRESHOLD is not None:"
 562                              , "        mask_t = numpy.where(img > THRESHOLD, True, False)"
 563                              , "        lst_mask.append(numpy.logical_or(mask, mask_t))"
 564                              , "    else:"
 565                              , "        lst_mask.append(mask)"
 566                              , ""
 567                              , "# Integration multi-geometry 1D"
 568                              , "mg = MultiGeometry(PONIES, unit=\"2th_deg\", radial_range=(0,80))"
 569                              , "p = mg.integrate1d(imgs, BINS, lst_mask=lst_mask, lst_flat=flat)"
 570                              , ""
 571                              , "# Save the datas"
 572                              , "numpy.savetxt(OUTPUT, numpy.array(p).T)"
 573                              ]
 574       (Nxs nxs' (XrdOneDH5Path (DataItemH5 i' _) _ _ _)) = difTomoFrameNxs f
 575       output = "multi.dat"
 576       (Geometry _ (Source w) _ _) = difTomoFrameGeometry f
 577       (idxs, ponies) = unzip idxPonies
 578
 579 createMultiPyEdf ∷ XrdOneDParams a → DIM1 → Maybe Threshold → [FilePath] → [FilePath] → FilePath → FilePath → (Script Py2, FilePath)
 580 createMultiPyEdf (XrdOneDParams _ mflat _) b mt edfs ponies scriptPath output = (Py2Script (content, scriptPath), output)
 581     where
 582       content = Text.unlines $
 583                 map Text.pack ["#!/bin/env python"
 584                               , ""
 585                               , "import numpy"
 586                               , "from fabio import open"
 587                               , "from pyFAI.multi_geometry import MultiGeometry"
 588                               , ""
 589                               , "EDFS = " ++ toPyVal edfs
 590                               , "PONIES = " ++ toPyVal ponies
 591                               , "BINS = " ++ toPyVal (size b)
 592                               , "OUTPUT = " ++ toPyVal output
 593                               , "THRESHOLD = " ++ toPyVal mt
 594                               , ""
 595                               , "# load the flat"
 596                               , "flat = " ++ toPyVal mflat
 597                               , ""
 598                               , "# Read all the images"
 599                               , "imgs = [open(edf).data for edf in EDFS]"
 600                               , ""
 601                               , "# Compute the mask"
 602                               , "mask = numpy.zeros_like(imgs[0], dtype=bool)"
 603                               , "if THRESHOLD is not None:"
 604                               , "    for img in imgs:"
 605                               , "        mask_t = numpy.where(img > THRESHOLD, True, False)"
 606                               , "        mask = numpy.logical_or(mask, mask_t)"
 607                               , ""
 608                               , "# Integration multi-geometry 1D"
 609                               , "mg = MultiGeometry(PONIES, unit=\"2th_deg\", radial_range=(0,80))"
 610                               , "p = mg.integrate1d(imgs, BINS, lst_mask=mask)"
 611                               , ""
 612                               , "# Save the datas"
 613                               , "numpy.savetxt(OUTPUT, numpy.array(p).T)"
 614                               ]
 615
 616 saveMulti' ∷ XrdOneDParams a → DIM1 → Maybe Threshold → Consumer (DifTomoFrame' sh) (StateT [(Int, FilePath)] IO) r
 617 saveMulti' p b mt = forever $ do
 618   idxPonies <- lift get
 619   f'@(DifTomoFrame' f@(DifTomoFrame _ idx _ _ _) poniPath) <- await
 620   let directory = takeDirectory poniPath
 621   let filename = directory </> "multi.py"
 622   let (script, _) = createMultiPy p b mt filename f' idxPonies
 623   ExitSuccess ← lift . lift $ if (difTomoFrameEOF f) then (run script True) else return ExitSuccess
 624   lift $ put $! (idxPonies ++ [(idx, poniPath)])
 625
 626 saveMultiGeometry ∷ XrdOneDParams a → DIM1 → Maybe Threshold → Consumer (DifTomoFrame' sh) IO r
 627 saveMultiGeometry p b mt = evalStateP [] (saveMulti' p b mt)
 628
 629
 630 -- substract a sample from another one
 631
 632 targetMulti' ∷  XrdOneDParams a → AbsDirPath → XrdNxs → (FilePath, FilePath)
 633 targetMulti' _ output (XrdNxs _ _ _ _ (XrdSourceNxs (Nxs f _))) = (d, o)
 634   where
 635     d = getScanDir output f
 636     o = d </> "multi.dat"
 637
 638 targetMulti ∷ XrdOneDParams a → XRDSample → [(FilePath, FilePath)]
 639 targetMulti p (XRDSample _ output nxss) = map (targetMulti' p output) nxss
 640
 641 substractMulti' ∷ XrdOneDParams a → XRDSample → XRDSample → IO ()
 642 substractMulti' p s1@(XRDSample name _ _) s2 = do
 643   -- compute the output of the s1 sample
 644   -- we take only the first list of the sample
 645   let f1s:_ = targetMulti p s1
 646   -- compute the output of the s2 sample
 647   let f2s = targetMulti p s2
 648   -- do the substraction via a python script and add the gnuplot file
 649   _ ← mapConcurrently (go f1s) f2s
 650
 651   return ()
 652   where
 653     go ∷ (FilePath, FilePath) → (FilePath, FilePath) → IO ()
 654     go (_, f1) (d, f2) = do
 655       -- compute the substracted output file names
 656       let outputs = dropExtension f2 ++ "-" ++ name <.> "dat"
 657       -- compute the script name
 658       let scriptPath = d </> "multi-substract.py"
 659       let script = substractPy [f1] [f2] [outputs] scriptPath
 660       ExitSuccess ← run script False
 661       -- gnuplot
 662       let gnuplotPath = d </> "multi-substract.gnuplot"
 663       scriptSave $ mkGnuplot [outputs] gnuplotPath
 664       return ()
 665
 666 substractMulti ∷  XrdOneDParams a → XRDSample → [XRDSample] → IO ()
 667 substractMulti p s ss = mapM_ (substractMulti' p s) ss