JUNO Opticks : Summer Progress

• Issue G4Cerenkov_modified::GetAverageNumberOfPhotons : slightly -ve
  • Alternative "S2" implementation developed
• PMTSimParamSvc::get_pmt_ce : angular efficiency > 1.
  • after MiaoYu bugfix, now NNVT+NNVT_HighQE > 1. at low angle
• Geant4 Meeting + Recent Opticks "Snapshot" tags
  • Started tagging to assist Geant4 Opticks example development
• Yuxiang progress on simulation matching
  • Problems seem focussed at front of CD sticks
• QUDARap/QCerenkov
  • Wavelength agreement achieved using double precision
  • Single precision technique using ICDF under testing
• CSGOptiX
  • started integrating OptiX 7 geometry with QUDARap simulation

For sanity need to make the leap to OptiX 7 ASAP : avoid dividing time

Frank-Tamm Formula : Cerenkov Photon Yield /mm at BetaInverse

                                                BetaInverse^2
          N_photon/369.81  =    Integral ( 1 - -----------------  )  dE         where   BetaInverse < ri[E]
                                                   ri(E)^2

                            =   Integral [ 1 ] dE -  BetaInverse^2 * Integral[ 1./(ri[E]*ri[E]) ] dE

G4Cerenkov::BuildThePhysicsTable -> CerenkovAngleIntegrals (misnomer)

• Integral[ 1/ri^2 ] dE cumulative trapezoidal approx. integral over RINDEX[E] of material

G4Cerenkov::GetAverageNumberOfPhotons assumes monotonic RINDEX + only one permitted energy region

636   G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
637
638   G4double dp, ge;
642   if (nMax < BetaInverse)        // ... no photons
649   else if (nMin > BetaInverse) {
650      dp = Pmax - Pmin;
651      ge = CAImax;
660   } else {
661      Pmin = Rindex->GetEnergy(BetaInverse);
662      dp = Pmax - Pmin;
667      G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
668      ge = CAImax - CAImin;
674   }
677   G4double NumPhotons = Rfact * charge/eplus * charge/eplus * (dp - ge * BetaInverse*BetaInverse);

• putting together the split integral leads to -ve NumPhotons when close to RINDEX peak
• G4PhysicsVector::Value applies linear interpolation to cumulative integral of 1/ri^2 <-- POOR APPROX

Alternative "s2" integral approach : more precise, simpler, faster

                 BetaInverse*BetaInverse
Integral [ 1. -  ----------------------- ] (for BetaInverse < RINDEX)
                 RINDEX * RINDEX

Integral [ 1. - cos^2 theta ]

Integral [ sin^2 theta ]

Integral [ s2 ]             ( s2 > 0 )

Do not split the integral, do "s2" integral in one pass. Advantages:

• avoids one level of linear approximation, ca
• cannot give -ve values
• simple one pass code, no separate find_crossings
• also s2 is faster other than when numPhotons is zero
  • AVOID by not calling when BetaInverse > RINDEX_max

G4Cerenkov_modified::GetAverageNumberOfPhotons_s2

G4double G4Cerenkov_modified::GetAverageNumberOfPhotons_s2(
    const G4double charge, const G4double beta, const G4Material*, G4MaterialPropertyVector* Rindex) const
{
    G4double BetaInverse = 1./beta;
    G4double s2integral(0.) ;

    for(unsigned i=0 ; i < Rindex->GetVectorLength()-1 ; i++)
    {
        G4double en_0 = Rindex->Energy(i)   ; G4double en_1 = Rindex->Energy(i+1) ;
        G4double ri_0 = (*Rindex)[i]        ; G4double ri_1 = (*Rindex)[i+1] ;
        G4double ct_0 = BetaInverse/ri_0    ; G4double ct_1 = BetaInverse/ri_1 ;
        G4double s2_0 = (1.-ct_0)*(1.+ct_0) ; G4double s2_1 = (1.-ct_1)*(1.+ct_1) ;

        G4bool cross = s2_0*s2_1 < 0. ;
        G4double en_cross =  cross ? en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0) : -1. ;
        // linear crossing more precision than s2 zeros

        if( s2_0 <= 0. and s2_1 <= 0. )  // no CK
        {
            // noop
        }
        else if( s2_0 < 0. and s2_1 > 0. )  // s2 becomes +ve within the bin  left edge triangle
        {
            s2integral +=  (en_1 - en_cross)*s2_1*0.5 ;
        }
        else if( s2_0 >= 0. and s2_1 >= 0. )   // s2 +ve across full bin    trapezoid
        {
            s2integral += (en_1 - en_0)*(s2_0 + s2_1)*0.5 ;
        }
        else if( s2_0 > 0. and s2_1 < 0. )  // s2 becomes -ve within the bin   right edge triangle
        {
            s2integral +=  (en_cross - en_0)*s2_0*0.5 ;
        }
    }
    const G4double Rfact = 369.81/(eV * cm);
    return Rfact * charge/eplus * charge/eplus * s2integral ;
}

scan_GetAverageNumberOfPhotons_cputime

test_GetAverageNumberOfPhotons_plot

scan_GetAverageNumberOfPhotons_plot_1.0000_2.0000

scan_GetAverageNumberOfPhotons_plot_1.4536_1.7930

scan_GetAverageNumberOfPhotons_plot_1.7000_1.8000

scan_GetAverageNumberOfPhotons_difference_plot

PMTSimParamSvc::get_pmt_ce

Angular scan/plot code to check efficiency as function of theta

• https://github.com/simoncblyth/j/blob/main/PMTAngular/PMTAngular.py
• https://github.com/simoncblyth/j/blob/main/PMTAngular/PMTAngular.sh

Following Miaoyu bug fix, efficiencies > 1. for NNVT and NNVT_HighQE

• a good problem to have !
• trivial kludge fix ? : std::min( 1.0 , ce )

PMTAngular_fig1

PMTAngular_fig2

Geant4 Meeting + Recent Opticks "Snapshot" tags

Aug 20 : meeting with Geant4 : Fermilab + Warwick(Ben Morgan, Software Management)

• regarding what is needed for Opticks advanced example in next Geant4
• requested Opticks to start providing regular snapshot tags

Recent Tags

• https://github.com/simoncblyth/opticks/tags
• PLANNED : https://github.com/simoncblyth/opticks/releases/tag/v0.1.3
  • Sept 2nd 2021
  • OPTICKS_VERSION_NUMBER 13
  • CManager FIX
• https://github.com/simoncblyth/opticks/releases/tag/v0.1.2
  • Aug 30 2021
  • OPTICKS_VERSION_NUMBER 12
  • fixed Cerenkov wavelenth bug : reported by Yuxiang
• https://github.com/simoncblyth/opticks/releases/tag/v0.1.1
  • Aug 28 2021
  • OPTICKS_VERSION_NUMBER 11
  • requested by Fermilab Geant4 team

Yuxiang Progress on Simulation Matching

• noticed Cerenkov wavelength bug, now fixed
• confirmed suspicion of CD Sticks geometry problem
  • removing "setupCD_sticks" gives agreement in history categories
• found history category "SI BT BT BT BT AB" 0x4cccc2 with strong Opticks excess when "setupCD_sticks" is present

NEXT:

• shine light on the discrepant geometry using input photons to see nature of problem

Opticks Packages (OptiX <= 6.5)

epsilon:~ blyth$ opticks-deps
[2021-09-01 18:17:04,496] p10852 {/Users/blyth/opticks/bin/CMakeLists.py:165} INFO - home /Users/blyth/opticks
      API_TAG :        reldir :         bash- :     Proj.name : dep Proj.names
10        OKCONF :        okconf :        okconf :        OKConf : OpticksCUDA OptiX G4
20        SYSRAP :        sysrap :        sysrap :        SysRap : OKConf NLJSON PLog
30          BRAP :      boostrap :          brap :      BoostRap : Boost BoostAsio NLJSON PLog SysRap Threads
40           NPY :           npy :           npy :           NPY : PLog GLM BoostRap
50        OKCORE :   optickscore :           okc :   OpticksCore : NPY
60          GGEO :          ggeo :          ggeo :          GGeo : OpticksCore
90         OKGEO :    opticksgeo :           okg :    OpticksGeo : OpticksCore GGeo


100       CUDARAP :       cudarap :       cudarap :       CUDARap : SysRap OpticksCUDA
110         THRAP :     thrustrap :         thrap :     ThrustRap : OpticksCore CUDARap
120         OXRAP :      optixrap :         oxrap :      OptiXRap : OKConf OptiX OpticksGeo ThrustRap
130          OKOP :          okop :          okop :          OKOP : OptiXRap


140        OGLRAP :        oglrap :        oglrap :        OGLRap : ImGui OpticksGLEW BoostAsio OpticksGLFW OpticksGeo
150          OKGL :     opticksgl :          okgl :     OpticksGL : OGLRap OKOP
160            OK :            ok :            ok :            OK : OpticksGL
165            X4 :         extg4 :            x4 :         ExtG4 : G4 GGeo OpticksXercesC CLHEP
170          CFG4 :          cfg4 :          cfg4 :          CFG4 : G4 ExtG4 OpticksXercesC OpticksGeo ThrustRap
180          OKG4 :          okg4 :          okg4 :          OKG4 : OK CFG4
190          G4OK :          g4ok :          g4ok :          G4OK : CFG4 ExtG4 OKOP
200          None :   integration :   integration :   Integration :
epsilon:~ blyth$

• new GPU/CPU geometry model required : CSGFoundry
• All GPU code needed to be replaced for OptiX 7
  • -> Four packages need to be almost entirely relaced (large fraction of Opticks)

New Opticks Packages (OptiX >= 7.0)

CSG

Basis geometry model

CSG_GGeo

Translation of GGeo into CSG model

QUDARap

Simulation Implementation, excluding geometry

• Scintillation generation
• Cerenkov generation

CSGOptiX

OptiX 7 geometry, depending on:

• CSG : geometry model
• QUDARap :simulation

QUDARap : new Heart of Opticks Simulation

• More modular approach than old OptiXRap/cu/generate.cu
• simple CPU/GPU "counterpart pattern" organization

	CPU	GPU header
context steering	QSim.hh	qsim.h
curandState setup	QRng.hh	qrng.h
property interpolation	QProp.hh	qprop.h
event handling	QEvent.hh	qevent.h
Cerenkov generation	QCerenkov.hh
Scintillation generation	QScint.hh
texture handling	QTex.hh

1. CPU : data access, uploads/downloads, prepares GPU instance with device pointers
2. GPU : simple header

Aims of counterpart code organization:

• facilitate fine-grained modular simulation testing
• keep bulk of GPU code in simple to test headers (some can be tested on CPU)
• QUDARap does not depend on OptiX -> more flexible -> simpler testing

QUDARap/QCerenkov Progress : single precision solution ?

• Cerenkov wavelength discrepancy
  • generally correct shape, but poor chi2
  • achieved good chi2, but only by using double precision
    • having to use double precision is a performance problem
• Single precision approach, using ~1000 x inverse-CDF for each BetaInverse
  • python prototype works (ana/rindex.py)
  • C++ CPU implementation done, currently testing
    • QUDARap/QCerenkov + NP::pdomain
    • https://github.com/simoncblyth/np/

QUDARap-CSGOptiX integration

CSGOptiX focus on geomerty

• uses CSG/CSGFoundry geometry model
• high level simulation, details from QUDARap via CUDA headers qsim.h etc..

Minimimize non-geometry code here, as OptiX dependency demands a very specific organization.

CSGOptiXSimulation : cross-section from planar gensteps

• test OptiX 7 geometry (CSGFoundry) + QUDARap genstep handling using 2d planar gensteps
• geometry x-z plane cross-section from thousands of intersects, acrylic + PMTs + Tyvek visible
• colors from boundary index of intersected geometry
• LATER (NOV 2021) : FOUND FUNNY PMT SHAPE ABOVE IS FROM A NOW FIXED BUG WITH ELLIPSOID IN CSG COMBINATION : SCALE TRANSFORM WAS BEING STOMPED

CSGOptiXSimulation : cross-section from planar gensteps : closeup

• zooming in you can see gaps between intersects
• grid of points are starting points of the 2d planar gensteps : inner surfaces may be missed if no rays can get to them
• geometry x-z plane cross-section from thousands of intersects

CSGOptiXSimulation : PMT geometry overlap+coincidence problems ?

• zooming closer : definitely looks like a geometry overlap problem
  • OR perhaps Opticks geometry translation bug

G4Cerenkov : Photon Energy/Wavelength + Cone angle generation

168 G4VParticleChange*
169 G4Cerenkov_modified::PostStepDoIt(
      const G4Track& aTrack, const G4Step& aStep)
...
252   G4double Pmin = Rindex->GetMinLowEdgeEnergy();
253   G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
254   G4double dp = Pmax - Pmin;
...
268   G4double maxCos = BetaInverse / nMax;
270   G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
...
315   for (G4int i = 0; i < fNumPhotons; i++) {
317
318       G4double rand;
319       G4double sampledEnergy, sampledRI;
320       G4double cosTheta, sin2Theta;
...
324       do {
325          rand = G4UniformRand();
326          sampledEnergy = Pmin + rand * dp;
327          sampledRI = Rindex->Value(sampledEnergy);
334          cosTheta = BetaInverse / sampledRI;
342          sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
343          rand = G4UniformRand();
344
346       } while (rand*maxSin2 > sin2Theta);

https://bitbucket.org/simoncblyth/opticks/src/master/examples/Geant4/CerenkovStandalone/G4Cerenkov_modifiedTest.cc

G4Cerenkov_modifiedTest_SKIP_CONTINUE

1.5/1.8 = .83, 1.5/1.62 = 0.92

ana/ck.py rejection sampling

Cerenkov energy/cone angle sample : resorting to double precision

574 template <typename T>
575 inline QCTX_METHOD void qctx<T>::cerenkov_photon_expt(
        quad4& p, unsigned id, curandStateXORWOW& rng, int print_id )
576 {
577     double BetaInverse = 1.5 ;
578     double Pmin = 1.55 ;
579     double Pmax = 15.5 ;
580     double nMax = 1.793 ;
581     double maxCos = BetaInverse / nMax;
582     double maxSin2 = ( 1. - maxCos )*( 1. + maxCos );
583
584     double u0 ;
585     double u1 ;
586     double energy ;
587     double sampledRI ;
588     double cosTheta ;
589     double sin2Theta ;
590     double u_mxs2_s2 ;
592     unsigned loop = 0u ;
593
594     do {
596         u0 = curand_uniform_double(&rng) ;
598         energy = Pmin + u0*(Pmax - Pmin) ;
600         sampledRI = prop->interpolate( 0u, energy );
602         cosTheta = BetaInverse / sampledRI ;
604         sin2Theta = (1. - cosTheta)*(1. + cosTheta);
606         u1 = curand_uniform_double(&rng) ;
608         u_mxs2_s2 = u1*maxSin2 - sin2Theta ;
610         loop += 1 ;
612     } while ( u_mxs2_s2 > 0. );

ck_photon_enprop 100 very poor chi2

ck_photon_expt 100 matched

ck_photon_enprop 1001 better but still poor chi2

ck_photon_expt 1001 matched