simulPET/README_resModelling.md

Resolution modeling

The helper function is loaded by:

root[].L getResolution.C

Individual resolution calculation

According to literature, the resolution contribution from spatial resolution of the crystals is calculated from sum of sines, which is coded in getResolution(ctheta0,N) where ctheta0 is the smalles cosine theta still captured by the system geometry (assumming coverage from 1), and N is number of splits in cos theta used in numeric integration.

root[] double detRes=2.;//resolution in mm
root[] double sigmaRSquared=detRes*detRes/getResolution(ctheta0,40);

In the same manner, resolution contribution for timing is calculated:

root[] double tRes=50.;//resolution in ps
root[] double xtRes=3e11*tRes*1e-12;//resolution in mm
root[] double sigmaTSquared=xtRes*xtRes/getResolutionT(ctheta0,40);

To generate individual graphs, do

TGraph gr(20);
for (int i=0;i<gr.GetN();i++) { 
    double ctheta0=1-(i+0.5)/(gr.GetN()); 
    double ttheta0=sqrt(1-pow(ctheta0,2))/ctheta0;
    gr.SetPoint(i,ttheta0,1./sqrt(getResolutionT(ctheta0,40)));
}
//equivalent for getResolution

Combined resolution calculation

root[] double tRes=50.;//resolution in ps
root[] double xtRes=3e11*tRes*1e-12;//resolution in mm
root[] double xRes=2;//resolution in mm
root [] for (int i=0;i<gr.GetN();i++) { 
    double ctheta0=1-(i+0.5)/(gr.GetN()); 
    double ttheta0=sqrt(1-pow(ctheta0,2))/ctheta0;
    double sigmaT=getResolutionT(ctheta0,40)/xtRes/xtRes;
    double sigmaX=getResolution(ctheta0,40)/xRes/xRes;
    gr.SetPoint(i,ttheta0,sqrt(1./sigmaT+1./sigmaX));
}