[bug] fix detp output bug after adding srcid, #163

README.md

@@ -1337,7 +1337,7 @@ to adpot the newly implemented JSON/.jdat format for easy data sharing.
 The .mch file contains a 256 byte binary header, followed by a 2-D numerical array
 of dimensions `#savedphoton * #colcount` as recorded in the header.
 ```
- typedef struct MCXHistoryHeader{
+typedef struct MCXHistoryHeader{
 	char magic[4];                 // magic bits= 'M','C','X','H'
 	unsigned int  version;         // version of the mch file format 
 	unsigned int  maxmedia;        // number of media in the simulation 
@@ -1346,14 +1346,15 @@ of dimensions `#savedphoton * #colcount` as recorded in the header.
 	unsigned int  totalphoton;     // how many total photon simulated 
 	unsigned int  detected;        // how many photons are detected (not necessarily all saved) 
 	unsigned int  savedphoton;     // how many detected photons are saved in this file 
-	float unitinmm;                // what is the voxel size of the simulation 
+	float unitinmm;                // what is the voxel size of the simulation
 	unsigned int  seedbyte;        // how many bytes per RNG seed
-        float normalizer;              // what is the normalization factor
+	float normalizer;              // what is the normalization factor
 	int respin;                    // if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset 
 	unsigned int  srcnum;          // number of sources for simultaneous pattern sources 
 	unsigned int  savedetflag;     // number of sources for simultaneous pattern sources 
-	int reserved[2];               // reserved fields for future extension 
- } History;
+    unsigned int  totalsource;     // total source number when multiple sources are defined
+	int reserved[1];               // reserved fields for future extension 
+} History;
 ```
 When the `-q` flag is set to 1, the detected photon initial seeds are also stored
 following the detected photon data, consisting of a 2-D byte array of `#savedphoton * #seedbyte`.
@@ -1369,36 +1370,36 @@ When `-F jnii` is specified, instead of saving the detected photon into the lega
 a .jdat file is written, which is a pure JSON file. This file contains a hierachical data
 record of the following JSON structure
 ````
- {
-   "MCXData": {
+{
+   "MCXData":{
        "Info":{
            "Version":
-	   "MediaNum":
-	   "DetNum":
-	   ...
-	   "Media":{
-	      ...
-	   }
+           "MediaNum":
+           "DetNum":
+           ...
+           "Media":{
+               ...
+           }
        },
        "PhotonData":{
            "detid":
-	   "nscat":
-	   "ppath":
-	   "mom":
-	   "p":
-	   "v":
-	   "w0":
+           "nscat":
+           "ppath":
+           "mom":
+           "p":
+           "v":
+           "w0":
        },
        "Trajectory":{
            "photonid":
-	   "p":
-	   "w0":
+           "p":
+           "w0":
        },
        "Seed":[
            ...
        ]
    }
- }
+}
 ````
 where "Info" is required, and other subfields are optional depends on users' input.
 Each subfield in this file may contain JData 1-D or 2-D array constructs to allow 

README.txt

@@ -1328,21 +1328,22 @@ The .mch file contains a 256 byte binary header, followed by a 2-D numerical arr
 of dimensions #savedphoton * #colcount as recorded in the header.
 
  typedef struct MCXHistoryHeader{
-	char magic[4];                 /**< magic bits= 'M','C','X','H' */
-	unsigned int  version;         /**< version of the mch file format */
-	unsigned int  maxmedia;        /**< number of media in the simulation */
-	unsigned int  detnum;          /**< number of detectors in the simulation */
-	unsigned int  colcount;        /**< how many output files per detected photon */
-	unsigned int  totalphoton;     /**< how many total photon simulated */
-	unsigned int  detected;        /**< how many photons are detected (not necessarily all saved) */
-	unsigned int  savedphoton;     /**< how many detected photons are saved in this file */
-	float unitinmm;                /**< what is the voxel size of the simulation */
-	unsigned int  seedbyte;        /**< how many bytes per RNG seed */
-        float normalizer;              /**< what is the normalization factor */
-	int respin;                    /**< if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset */
-	unsigned int  srcnum;          /**< number of sources for simultaneous pattern sources */
-	unsigned int  savedetflag;     /**< number of sources for simultaneous pattern sources */
-	int reserved[2];               /**< reserved fields for future extension */
+	char magic[4];                 // magic bits= 'M','C','X','H'
+	unsigned int  version;         // version of the mch file format 
+	unsigned int  maxmedia;        // number of media in the simulation 
+	unsigned int  detnum;          // number of detectors in the simulation 
+	unsigned int  colcount;        // how many output files per detected photon 
+	unsigned int  totalphoton;     // how many total photon simulated 
+	unsigned int  detected;        // how many photons are detected (not necessarily all saved) 
+	unsigned int  savedphoton;     // how many detected photons are saved in this file 
+	float unitinmm;                // what is the voxel size of the simulation
+	unsigned int  seedbyte;        // how many bytes per RNG seed
+	float normalizer;              // what is the normalization factor
+	int respin;                    // if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset 
+	unsigned int  srcnum;          // number of sources for simultaneous pattern sources 
+	unsigned int  savedetflag;     // number of sources for simultaneous pattern sources 
+    unsigned int  totalsource;     // total source number when multiple sources are defined
+	int reserved[1];               // reserved fields for future extension 
  } History;
 
 When the `-q` flag is set to 1, the detected photon initial seeds are also stored
@@ -1360,29 +1361,29 @@ a .jdat file is written, which is a pure JSON file. This file contains a hierach
 record of the following JSON structure
 
  {
-   "MCXData": {
+   "MCXData":{
        "Info":{
            "Version":
-	   "MediaNum":
-	   "DetNum":
-	   ...
-	   "Media":{
-	      ...
-	   }
+           "MediaNum":
+           "DetNum":
+           ...
+           "Media":{
+               ...
+           }
        },
        "PhotonData":{
            "detid":
-	   "nscat":
-	   "ppath":
-	   "mom":
-	   "p":
-	   "v":
-	   "w0":
+           "nscat":
+           "ppath":
+           "mom":
+           "p":
+           "v":
+           "w0":
        },
        "Trajectory":{
            "photonid":
-	   "p":
-	   "w0":
+           "p":
+           "w0":
        },
        "Seed":[
            ...

pmcx/Makefile

@@ -0,0 +1,6 @@
+all:
+	python3 -m pip wheel .
+clean:
+	rm -rf *.whl
+install:
+	python3 -m pip install --force-reinstall pmcx-*.whl

pmcx/README.md

@@ -4,7 +4,7 @@
 
 - Copyright: (C) Matin Raayai Ardakani (2022-2023) <raayaiardakani.m at northeastern.edu>, Qianqian Fang (2019-2023) <q.fang at neu.edu>, Fan-Yu Yen (2023) <yen.f at northeastern.edu>
 - License: GNU Public License V3 or later
-- Version: 0.2.10
+- Version: 0.2.11
 - URL: https://pypi.org/project/pmcx/
 - Github: https://github.com/fangq/mcx
 

pmcx/pmcx/__init__.py

@@ -49,7 +49,7 @@
 # from .files import loadmc2, loadmch, load, save
 from .bench import bench
 
-__version__ = "0.2.10"
+__version__ = "0.2.11"
 
 __all__ = (
     "gpuinfo",

pmcx/pmcx/utils.py

@@ -357,6 +357,7 @@ def detphoton(detp, medianum, savedetflag, issaveref=None, srcnum=None):
                  newdetp['p'] or ['v']: exit position and direction, when cfg.issaveexit=1 (3)
                  newdetp['w0']: photon initial weight at launch time (3)
                  newdetp['s']: exit Stokes parameters for polarized photon (4)
+                 newdetp['srcid']: the ID of the source when multiple sources are defined (1)
     """
     newdetp = {}
     c0 = 0
@@ -367,6 +368,9 @@ def detphoton(detp, medianum, savedetflag, issaveref=None, srcnum=None):
             newdetp["w0"] = detp[0, :].transpose()
         else:
             newdetp["detid"] = detp[0, :].astype(int).transpose()
+            if np.any(newdetp["detid"] > 65535):
+                newdetp["srcid"] = np.right_shift(newdetp["detid"], 16, dtype=np.int32)
+                newdetp["detid"] = np.bitwise_and(newdetp["detid"], 0xff, dtype=np.int32)
         c0 = 1
 
     length = medianum

pmcx/setup.py

@@ -123,7 +123,7 @@ def build_extension(self, ext):
 setup(
     name="pmcx",
     packages=["pmcx"],
-    version="0.2.10",
+    version="0.2.11",
     requires=["numpy"],
     license="GPLv3+",
     author="Matin Raayai Ardakani, Qianqian Fang, Fan-Yu Yen",

src/mcx_core.cu

@@ -374,11 +374,11 @@ __device__ inline void savedetphoton(float n_det[], uint* detectedphoton, float*
             baseaddr *= gcfg->reclen;
 
             if (SAVE_DETID(gcfg->savedetflag)) {
-                n_det[baseaddr++] = detid;
-
                 if (gcfg->extrasrclen && gcfg->srcid <= 0) {
-                    n_det[baseaddr++] = (((int)ppath[gcfg->w0offset - 1]) << 16);
+                    detid |= (((int)ppath[gcfg->w0offset - 1]) << 16);
                 }
+
+                n_det[baseaddr++] = detid;
             }
 
             for (i = 0; i < gcfg->partialdata; i++) {

src/mcx_utils.c

@@ -1024,12 +1024,13 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c
     }
 
     if (cfg->his.detected == 0  && cfg->his.savedphoton) {
-        char colnum[] = {1, 3, 1};
-        char* dtype[] = {"uint32", "single", "single"};
-        char* dname[] = {"photonid", "p", "w0"};
+        char colnum[] = {1, 3, 1, 1};
+        char* dtype[] = {"uint32", "single", "single", "uint32"};
+        char* dname[] = {"photonid", "p", "w0", "srcid"};
+        int activecol = sizeof(colnum) - (cfg->his.totalsource == 1);
         cJSON_AddItemToObject(obj, "Trajectory", dat = cJSON_CreateObject());
 
-        for (int id = 0; id < sizeof(colnum); id++) {
+        for (int id = 0; id < activecol; id++) {
             uint dims[2] = {count, colnum[id]};
             float* buf = (float*)calloc(dims[0] * dims[1], sizeof(float));
 
@@ -1048,9 +1049,9 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c
             col += dims[1];
         }
     } else {
-        char colnum[] = {1, cfg->his.maxmedia, cfg->his.maxmedia, cfg->his.maxmedia, 3, 3, 1};
-        char* dtype[] = {"uint32", "uint32", "single", "single", "single", "single", "single"};
-        char* dname[] = {"detid", "nscat", "ppath", "mom", "p", "v", "w0"};
+        char colnum[] = {1, cfg->his.maxmedia, cfg->his.maxmedia, cfg->his.maxmedia, 3, 3, 1, 4};
+        char* dtype[] = {"uint32", "uint32", "single", "single", "single", "single", "single", "single"};
+        char* dname[] = {"detid", "nscat", "ppath", "mom", "p", "v", "w0", "s"};
         cJSON_AddItemToObject(obj, "PhotonData", dat = cJSON_CreateObject());
 
         for (int id = 0; id < sizeof(colnum); id++) {
@@ -1059,13 +1060,18 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c
                 void* val = NULL;
                 float* fbuf = NULL;
                 uint*  ibuf = NULL;
+                int hassrcid = ((cfg->his.totalsource > 1) && id == 1);
 
                 if (!strcmp(dtype[id], "uint32")) {
                     ibuf = (uint*)calloc(dims[0] * dims[1], sizeof(uint));
 
                     for (int i = 0; i < dims[0]; i++)
                         for (int j = 0; j < dims[1]; j++) {
                             ibuf[i * dims[1] + j] = ppath[i * cfg->his.colcount + col + j];
+
+                            if (hassrcid) {
+                                ibuf[i * dims[1] + j] &= 0xFFFF;
+                            }
                         }
 
                     val = (void*)ibuf;
@@ -1080,6 +1086,23 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c
                     val = (void*)fbuf;
                 }
 
+                if (hassrcid) {
+                    uint* srcid = (uint*)calloc(dims[0] * dims[1], sizeof(uint));
+
+                    for (int i = 0; i < dims[0] * dims[1]; i++) {
+                        srcid[i] = ibuf[i] >> 16;
+                        ibuf[i] &= 0xFFFF;
+                    }
+
+                    cJSON_AddItemToObject(dat, "srcid", sub = cJSON_CreateObject());
+
+                    if (mcx_jdataencode((void*)srcid, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, 0, cfg)) {
+                        MCX_ERROR(-1, "error when converting to JSON");
+                    }
+
+                    free(srcid);
+                }
+
                 cJSON_AddItemToObject(dat, dname[id], sub = cJSON_CreateObject());
 
                 if (mcx_jdataencode(val, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, 0, cfg)) {

utils/loadmch.m

@@ -28,7 +28,7 @@
 %                    [detid(1) nscat(M) ppath(M) mom(M) p(3) v(3) w0(1) s(4)]
 %        header: file header info, a structure has the following fields
 %                [version,medianum,detnum,recordnum,totalphoton,detectedphoton,
-%                 savedphoton,lengthunit,seedbyte,normalizer,respin,srcnum,savedetflag]
+%                 savedphoton,lengthunit,seedbyte,normalizer,respin,srcnum,savedetflag,totalsource]
 %        photonseed: (optional) if the mch file contains a seed section, this
 %                returns the seed data for each detected photon. Each row of 
 %                photonseed is a byte array, which can be used to initialize a  
@@ -70,7 +70,8 @@
 	respin=fread(fid,1,'int');
 	srcnum=fread(fid,1,'uint');
 	savedetflag=fread(fid,1,'uint');
-	junk=fread(fid,2,'uint');
+	totalsource=fread(fid,1,'uint');
+	junk=fread(fid,1,'uint');
 
     detflag=dec2bin(bitand(savedetflag,(2^8-1)))-'0';
     if(strcmp(endian,'ieee-le'))
@@ -119,5 +120,5 @@
                        'recordnum',header(4),'totalphoton',header(5),...
                        'detectedphoton',header(6),'savedphoton',header(7),...
                        'lengthunit',header(8),'seedbyte',seedbyte,'normalizer',normalizer,...
-		       'respin',respin,'srcnum',srcnum,'savedetflag',savedetflag);
+		       'respin',respin,'srcnum',srcnum,'savedetflag',savedetflag,'totalsource',totalsource);
 end

utils/mcxdetphoton.m

@@ -38,6 +38,10 @@
         newdetp.w0=detp(1,:)';
     else
         newdetp.detid=int32(detp(1,:))';
+        if(any(newdetp.detid>65535))
+            newdetp.srcid=bitshift(newdetp.detid,-16);
+            newdetp.detid=bitand(newdetp.detid,int32(hex2dec('ffff')));
+        end
     end
     c0=2;
 end