Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 112 → Rev 113

/petrecofmf/PETProjDataMgr.h
File deleted
/petrecofmf/m16_reverse.map
File deleted
/petrecofmf/Makefile
File deleted
/petrecofmf/plot.cxx
File deleted
/petrecofmf/sumpedestals_zero.dat
File deleted
/petrecofmf/m16.map
File deleted
/petrecofmf/sumpedestals.dat
File deleted
/petrecofmf/modules.map
File deleted
/petrecofmf/pedestals.cxx
File deleted
/petrecofmf/readdata.C
File deleted
/petrecofmf/xpath2.c
File deleted
/petrecofmf/fotovrh.cxx
File deleted
/petrecofmf/calibration.cxx
File deleted
/petrecofmf/pedestals_zero.dat
File deleted
/petrecofmf/pedestals.dat
File deleted
/petrecofmf/README
File deleted
/petrecofmf/fotovrh.dat
File deleted
/petrecofmf/calibration_all.root
Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream
Property changes:
Deleted: svn:mime-type
## -1 +0,0 ##
-application/octet-stream
\ No newline at end of property
Index: petrecofmf/config_zero.xml
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
/petrecofmf/config_zero.xml
Property changes:
Deleted: svn:mime-type
## -1 +0,0 ##
-application/xml
\ No newline at end of property
Index: petrecofmf/PETProjDataMgr.C
===================================================================
--- petrecofmf/PETProjDataMgr.C (revision 112)
+++ petrecofmf/PETProjDataMgr.C (nonexistent)
@@ -1,612 +0,0 @@
-//
-// ../bin/FBP2D FBP2D.par
-//
-// ********************************************************************
-// * License and Disclaimer *
-// * *
-// * The GAMOS software is copyright of the Copyright Holders of *
-// * the GAMOS Collaboration. It is provided under the terms and *
-// * conditions of the GAMOS Software License, included in the file *
-// * LICENSE and available at http://fismed.ciemat.es/GAMOS/license .*
-// * These include a list of copyright holders. *
-// * *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work make any representation or warranty, express or implied, *
-// * regarding this software system or assume any liability for its *
-// * use. Please see the license in the file LICENSE and URL above *
-// * for the full disclaimer and the limitation of liability. *
-// * *
-// * This code implementation is the result of the scientific and *
-// * technical work of the GAMOS collaboration. *
-// * By using, copying, modifying or distributing the software (or *
-// * any work based on the software) you agree to acknowledge its *
-// * use in resulting scientific publications, and indicate your *
-// * acceptance of all terms of the GAMOS Software license. *
-// ********************************************************************
-//
-#include "PETProjDataMgr.h"
-#include "TH2F.h"
-#include "TH3F.h"
-#include "TMath.h"
-#include "TRandom3.h"
-
-#include <iomanip>
-#include <iostream>
-#include <cstdlib>
-
-using namespace std;
-
-//----------------------------------------------------------------------
-PETProjDataMgr* PETProjDataMgr::m_Instance = 0;
-
-//----------------------------------------------------------------------
-PETProjDataMgr* PETProjDataMgr::GetInstance()
-{
- if( !m_Instance ) {
- m_Instance = new PETProjDataMgr;
- }
-
- return m_Instance;
-
-}
-
-//-----------------------------------------------------------------------
-PETProjDataMgr::PETProjDataMgr()
-{
-
- /*
-// ARGUMENTS:
-" cout << " -------------------------- \n"
-" Arguments convention: \n"
-" -a Axial FOV (mm), <theDist_axial=100.0> \n"
-" -d Diameter Transaxial FOV (mm), <m_RingDiameter=300.0> \n"
-" -i Type of the input file (by default: 0 = Arce_binary), <typeINfile=0> \n"
-" \n"
-" -m Maximum ring difference (by default: -1 = m_NOfPlanes), <m_MaxRingDifferenceiff> \n"
-" -n Name of output file, <m_Filename> \n"
-" -p Axial number of planes, <m_NOfPlanes> \n"
-" -r Number of bins, \"distancias\", <m_NOfBins> \n"
- // -s Span TO DO:span !!!!!
-" -t Number of angular views, \"direcciones\", <m_NOfAngles> \n"
-" -v Verbosity (by default: 0=silent, 3=debug), <verbos> \n"
-" -x Maximum number of coincidences to be stored (by default: -1 = no limit), <Max_Coinci> \n"
-" -o Output type (by default: 0 = mcc Interfile, 1 = STIR Interfile), <OutType> \n"
-
-" \n"
-" PET Reconstruction. CIEMAT 2009-11 \n"
-" mario.canadas@ciemat.es \n"
- " -------------------------- \n";
- */
- m_rnd= new TRandom3();
- m_AxialDistance = (9-1)*2.25; // Axial pixel dimension*NOfPlanes
- m_RingDiameter = 120.0; // notranji premer peta
- m_NOfPlanes = 9; // stevilo ravnin
- m_NOfBins = 128; // stevilo binov v razdalji
- m_nch = 128; // stevilo padov okoli in okoli
- m_NOfAngles = TMath::Abs(m_nch)/2; // stevilo kotov = stevilo padov okoli in okoli /2
-
- m_MaxRingDifference = -1;
- //m_MaxRingDifference = 3; // najvecja razdalja med padi
- // toDo: theSpan = int(GmParameterMgr::GetInstance()->GetNumericValue("PET:ProjData:Span",1));
-
-
- m_OutFormat = 1; // 1.. projections 0 .. image
- m_Filename = "sino3D";
- m_Debug=1;
-
- if (m_MaxRingDifference==-1) m_MaxRingDifference=m_NOfPlanes-1;
-
-
- m_TotalAxialPlanes=m_NOfPlanes*m_NOfPlanes;
- if (m_OutFormat==1) m_TotalAxialPlanes= (2*m_NOfPlanes-1 - m_MaxRingDifference)*m_MaxRingDifference + m_NOfPlanes; // total number of Axial planes (segments*planes) in STIR format
-
- /*--- Initialize sino3D ---*/
- m_projections = new SINO_TYPE**[m_NOfBins];
- for(int i=0;i<m_NOfBins;i++){
- m_projections[i] = new SINO_TYPE*[m_NOfAngles];
- for(int j=0;j<m_NOfAngles;j++){
- m_projections[i][j] = new SINO_TYPE[m_TotalAxialPlanes]; /// ! If m_OutFormat==1 (STIR output):Matrix size depends on the MAX_Ring_Difference
- for(int k=0;k<m_TotalAxialPlanes;k++){
- m_projections[i][j][k]=0;
- }
- }
- }
-
- m_TotalProjectionCoincidences=0;
- m_TotalCoincidences=0;
- //OutputType = "pet";
-}
-
-//-----------------------------------------------------------------------
-void PETProjDataMgr::SetProjection( int axialplane, TH2F * h)
-{
-for(int i=0;i<m_NOfBins;i++){
- for(int j=0;j<m_NOfAngles;j++){
- m_projections[i][j][axialplane]=h->GetBinContent(i+1,j+1);
- }
-}
-
-}
-//-----------------------------------------------------------------------
-
-//-----------------------------------------------------------
-// from Gate
-//------------------------------------------------------------
-double ComputeSinogramS(double X1, double Y1, double X2, double Y2)
-{
- double s;
-
- double denom = (Y1-Y2) * (Y1-Y2) + (X2-X1) * (X2-X1);
-
- if (denom!=0.) {
- denom = sqrt(denom);
- s = ( X1 * (Y2-Y1) + Y1 * (X1-X2) ) / denom;
- } else {
- s = 0.;
- }
-
- double theta;
- if ((X1-X2)!=0.) {
- theta=atan((X1-X2) /(Y1-Y2));
- } else {
- theta=3.1416/2.;
- }
- if ((theta > 0.) && ((X1-X2) > 0.)) s = -s;
- if ((theta < 0.) && ((X1-X2) < 0.)) s = -s;
- if ( theta < 0.) {
- theta = theta+3.1416;
- s = -s;
- }
- return s;
-}
-
-
-void PETProjDataMgr::AddEvent( const TVector3& pos1 , const TVector3& pos2)
-{
- int z1_i, z2_i;
- //for discretization on the crystal: int x1_i, x2_i, y1_i, y2_i;
-
- double z1_abs=pos1.z()+m_AxialDistance/2;
- double z2_abs=pos2.z()+m_AxialDistance/2;
- double a, b, phi, dis;
- int phi_i, dis_i;
- int ring_diff;
-
- double _PI=2*asin(1);
-
- m_TotalCoincidences++;
-
- z1_i=(int)(m_NOfPlanes* z1_abs/m_AxialDistance); //round --> m_NOfPlanes+1 ...
- z2_i=(int)(m_NOfPlanes* z2_abs/m_AxialDistance);
-
- // control; if z_i out of range: return
-
- if ( (pos1.x()==pos2.x()) && (pos1.y()==pos2.y()) ) {
-#ifndef GAMOS_NO_VERBOSE
- if( m_Debug ) {
- cout << "PETProjDataMgr::AddEvent:WARNING! Event_1 == Event_2 ; x= " << pos2.x() << ", y= " << pos2.y() << ", z= " << pos2.z() << endl;
- }
-#endif
- return;
- }
-
- if ( (z1_i<0) || (z2_i<0) || (z1_i>= m_NOfPlanes) || (z2_i>= m_NOfPlanes) ) {
-#ifndef GAMOS_NO_VERBOSE
- if( m_Debug ) {
- cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Axial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
- }
-#endif
- return;
- }
-
- ring_diff = (int)fabs(z1_i-z2_i);
-
- // max ring difference; control:
- if (ring_diff > m_MaxRingDifference) {
-#ifndef GAMOS_NO_VERBOSE
- if( m_Debug ) {
- cout <<"PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Max. Ring Diff.): " << ring_diff << ">" << m_MaxRingDifference << " x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
- }
-#endif
- return;
- }
-
- a=(double)(pos2.y()- pos1.y());
- b=(double)(pos2.x()- pos1.x());
-
- if (a==0.0){
- phi=_PI*0.5;
- }
- else{
- phi=atan(b/a);
- }
-
- if (phi<0) phi = phi +_PI;
-
- dis=pos1.x()*cos(phi) - pos1.y()*sin(phi);
- //dis=ComputeSinogramS(pos1.x(), pos1.y(), pos2.x(), pos2.x());
- // control; transaxial FOV
- if ( fabs(dis) > m_RingDiameter*0.5 ) {
-#ifndef GAMOS_NO_VERBOSE
- if( m_Debug ) {
- cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Transaxial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
- }
-#endif
- return;
- }
-
- dis = dis + m_RingDiameter*0.5;
-
- // discret values:
- phi_i=(int)round( (double)(m_NOfAngles-1)*phi/_PI );
- dis_i=(int)round( (double)(m_NOfBins-1)*dis/(double)m_RingDiameter );
-
- if ((phi_i>=m_NOfAngles) || (dis_i>=m_NOfBins)) return; // only possible "=" because 'round' check it..
-
- // OLD: (SRRB included) sino3D[dis_i][phi_i][ (z1_i+z2_i)+ring_diff*(m_NOfPlanes-1) ]++;
-
- int Zpos;
-
- if (m_OutFormat==0) {
- Zpos = (z1_i*m_NOfPlanes + z2_i);
- }
- else{
-
- if (z1_i>=z2_i) { // SIN Max Ring_Diff: Zpos= ( ((m_NOfPlanes-ring_diff)*(m_NOfPlanes-1-ring_diff))/2 + z2_i );
-
- Zpos= ( ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff)*(m_MaxRingDifference - ring_diff))/2 + z2_i);
-
- }else{
- Zpos= ( (m_TotalAxialPlanes) - ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff +1)*(m_MaxRingDifference - ring_diff +1))/2 + z1_i );
-
- }
- }
-
- m_projections[dis_i][phi_i][ Zpos ]++;
- m_TotalProjectionCoincidences++;
-
-#ifndef GAMOS_NO_VERBOSE
- if( m_Debug >1) {
- cout << "PETProjDataMgr::AddEvent: x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
- cout << "PETProjDataMgr::AddEvent: Sinogram pos.: distance(s)= " << dis_i << ", angular view(phi)= " << phi_i << " ; Zpos=" << Zpos <<"; Segment (Ring diff.) = " << ring_diff << endl;
- }
-#endif
-
-
-}
-
-//-----------------------------------------------------------------------
-PETProjDataMgr::~PETProjDataMgr()
-{
- int i,j;
-
- for(i=0;i<m_NOfBins;i++){
- for(j=0;j<m_NOfAngles;j++){
- free(m_projections[i][j]);
-
- }
- free(m_projections[i]);
-
- }
- free(m_projections);
-
-}
-
-/* TO DO: call lm_to_sino3D program
-//-----------------------------------------------------------------------
-void PETIOMgr::ReadFile()
-{
- if( !theFileIn ) OpenFileIn();
-
- PETOutput po;
- G4bool bEof;
- for(;;) {
- po = ReadEvent( bEof );
- // theFileIn->read(reinterpret_cast<char *>(&po),sizeof(PetOutput));
- if(bDumpCout) PrintEvent(" PETOutput: ", po, bCartesian);
- if( bEof ) break;
- }
-}
-//-----------------------------------------------------------------------
-PETOutput PETProjDataMgr::ReadEvent( G4bool& bEof )
-{
- if( theFileIn == 0 ){
- G4Exception("PETIOMgr::ReadEvent, file not opened, call OpenFileIn() first ");
- }
-
- PETOutput po;
- fread (&po, sizeof(struct PETOutput),1,theFileIn);
- if ( feof( theFileIn ) ) {
- bEof = TRUE;
- } else {
- bEof = FALSE;
- }
-
- return po;
-
-}
-*/
-
-
-//-----------------------------------------------------------------------
-void PETProjDataMgr::WriteInterfile()
-{
-
- char name_hv[512];
- char name_v[512];
-
- if (m_OutFormat==0){
-
- strcpy(name_hv, m_Filename);
- strcpy(name_v,m_Filename);
- strcat(name_hv, ".hv");
- strcat(name_v, ".v");
-
- fp=fopen(name_hv, "w");
-
- fprintf (fp, "!INTERFILE := \n");
- fprintf (fp, "name of data file := %s\n", name_v);
- fprintf (fp, "!GENERAL DATA := \n");
- fprintf (fp, "!GENERAL IMAGE DATA :=\n");
- fprintf (fp, "!type of data := tomographic\n");
- fprintf (fp, "!version of keys := 3.3\n");
- fprintf (fp, "!data offset in bytes := 0\n");
- fprintf (fp, "imagedata byte order := littleendian\n");
- fprintf (fp, "!PET STUDY (General) :=\n");
- fprintf (fp, "!PET data type := 3D-Sinogram\n");
- fprintf (fp, "process status := Reconstructed\n");
- fprintf (fp, "!number format := unsigned short\n");
- fprintf (fp, "!number of bytes per pixel := 2\n");
- fprintf (fp, "number of dimensions := 3\n");
- fprintf (fp, "matrix axis label [1] := x\n");
- fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);
- fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(m_RingDiameter/(m_NOfBins-1.0)));
-
- fprintf (fp, "matrix axis label [2] := y\n");
- fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);
-
- fprintf (fp, "scaling factor (degree/pixel) [2] := %f\n",(float)(360./m_NOfAngles));
-
- fprintf (fp, "matrix axis label [3] := z\n");
- fprintf (fp, "!matrix size [3] := %i\n",m_NOfPlanes*m_NOfPlanes);
- fprintf (fp, "scaling factor (mm/pixel) [3] := %f\n",(float)(m_AxialDistance/(m_NOfPlanes-1.0)));
-
- fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);
- fprintf (fp, "number of time frames := 1\n");
- fprintf (fp, "image scaling factor[1] := 1\n");
- fprintf (fp, "data offset in bytes[1] := 0\n");
- fprintf (fp, "quantification units := 1\n");
- fprintf (fp, "!END OF INTERFILE := \n");
-
- fclose(fp);
- //(size_t)(m_NOfBins*m_NOfAngles*m_NOfPlanes*m_NOfPlanes);
-
- }else{
-
- strcpy(name_hv, m_Filename);
- strcpy(name_v,m_Filename);
-
- strcat(name_hv, ".hs"); // STIR extension: .hs .s
- strcat(name_v, ".s");
-
- fp=fopen(name_hv, "w");
-
- fprintf (fp, "!INTERFILE := \n");
- fprintf (fp, "name of data file := %s\n",name_v);
- fprintf (fp, "!GENERAL DATA := \n");
- fprintf (fp, "!GENERAL IMAGE DATA :=\n");
- fprintf (fp, "!type of data := PET\n");
- // fprintf (fp, "!version of keys := 3.3\n"); STIR format is not 3.3 (almost but not completely), ERROR in STIR if it is not removed
- // fprintf (fp, "!data offset in bytes := 0\n");
- fprintf (fp, "imagedata byte order := littleendian\n");
- fprintf (fp, "!PET STUDY (General) :=\n");
- fprintf (fp, "!PET data type := Emission\n");
- fprintf (fp, "applied corrections := {arc correction}\n"); // {none}\n");
- // fprintf (fp, "process status := Reconstructed\n");
- fprintf (fp, "!number format := unsigned integer\n");
- fprintf (fp, "!number of bytes per pixel := 2\n");
-
- fprintf (fp, "number of dimensions := 4\n");
- fprintf (fp, "matrix axis label [4] := segment\n");
- fprintf (fp, "!matrix size [4] := %i\n",m_MaxRingDifference*2 + 1);
- // fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(d_FOV/(m_NOfBins-1)));
- fprintf (fp, "matrix axis label [3] := axial coordinate\n");
- fprintf (fp, "!matrix size [3] := { ");
- if (m_MaxRingDifference==0)
- {
- fprintf (fp, "%i}\n", m_NOfPlanes);
- }else{
- for(int m=m_NOfPlanes-m_MaxRingDifference;m<=m_NOfPlanes;m++) fprintf (fp, "%i,", m);
- for(int m=m_NOfPlanes-1;m>m_NOfPlanes-m_MaxRingDifference;m--) fprintf (fp, "%i,", m);
- fprintf (fp, "%i}\n", m_NOfPlanes-m_MaxRingDifference);
- }
- fprintf (fp, "matrix axis label [2] := view\n");
- fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);
- fprintf (fp, "matrix axis label [1] := tangential coordinate\n");
- fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);
-
- fprintf (fp, "minimum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable)
- fprintf (fp, "%i", -m_MaxRingDifference);
- for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m);
- fprintf (fp, "}\n");
- fprintf (fp, "maximum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable)
- fprintf (fp, "%i", -m_MaxRingDifference);
- for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m);
- fprintf (fp, "}\n");
-
- fprintf (fp, "inner ring diameter (cm) := %f\n", m_RingDiameter/10); // STIR Required parameter, now assigned to FOV (not detectors)
- fprintf (fp, "average depth of interaction (cm) := 0.0001\n");
- fprintf (fp, "default bin size (cm) := %f\n",0.1*((float)m_RingDiameter/((float)m_NOfBins-1.0)) );
- fprintf (fp, "number of rings := %i\n",m_NOfPlanes );
- fprintf (fp, "distance between rings (cm) := %f\n", 0.1*((float)m_AxialDistance/(float)(m_NOfPlanes-1)) ); // Axial pixel dimension
-
- fprintf (fp, "number of detectors per ring := %i\n",m_NOfAngles*2 );
- // fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);
- fprintf (fp, "number of time frames := 1\n");
- fprintf (fp, "image scaling factor[1] := 1\n");
- fprintf (fp, "data offset in bytes[1] := 0\n");
- fprintf (fp, "quantification units := 1\n");
- fprintf (fp, "!END OF INTERFILE := \n");
-
- fclose(fp);
-
- }
- m_Buffer = (SINO_TYPE*) malloc( m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE));
-
- long unsigned int cont=0;
- int i,j,k;
-
- for(k=0;k<m_TotalAxialPlanes;k++){
- for(j=0;j<m_NOfAngles;j++){
- for(i=0;i<m_NOfBins;i++){
- m_Buffer[cont]=m_projections[i][j][k];
- cont++;
- }
- }
- }
-
- fp=fopen(name_v, "wb");
-
- //cout << 4096*sizeof(SINO_TYPE) << endl;
- int nb=fwrite(m_Buffer,1,m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE), fp);
- fclose(fp);
-
-#ifndef GAMOS_NO_VERBOSE
- cout << "PETProjDataMgr::WriteInterfile: File name: " << m_Filename << endl;
- cout << "PETProjDataMgr::WriteInterfile: Numer of bytes written: " << nb << endl;
- cout << "PETProjDataMgr::WriteInterfile: Planes = " << m_NOfPlanes << "; bins = " << m_NOfBins << "; ang_views = " << m_NOfAngles << endl;
- cout << "PETProjDataMgr::WriteInterfile: Dimensions (mm): Transaxial FOV = " << m_RingDiameter << "; Axial FOV = " << m_AxialDistance << " ; Transaxial_pix = " << m_RingDiameter/(m_NOfBins-1) <<"; Plane width = " << m_AxialDistance/(m_NOfPlanes-1) << endl; // Image Axial Pixel(ssrb) == 0.5*(Plane_Width);
- cout << "... " << endl;
-
- cout << "PETProjDataMgr::WriteInterfile: Total Coinci: " << m_TotalCoincidences << endl;
- cout << "PETProjDataMgr::WriteInterfile: Sino3D Coinci: " << m_TotalProjectionCoincidences << endl;
-#endif
-
-}
-
-
-TVector3 PETProjDataMgr::Hits2Digits(const TVector3 &r){
- if (!m_nch) return r;
- float smear=0.5;
-
- if (m_nch<0) smear=m_rnd->Rndm();
-
- double angle = TMath::ATan2(r.X(),r.Y()); // vrne kot med -pi in pi
- if (angle<0) angle+=TMath::TwoPi();
-
- angle= (int(angle/TMath::TwoPi()*TMath::Abs(m_nch))+smear)*TMath::TwoPi()/TMath::Abs(m_nch);
- //(m_rnd->Rndm()-0.5)*m_AxialDistance;
- return TVector3(sin(angle), cos(angle),0); // z coordinata ni cisto v redu
-
-}
-
-int PETProjDataMgr::FwdProject(double x,double y, double z, int nmax, TH1*h){
-TVector3 r(x,y,z);
-int h2d=h->InheritsFrom("TH2F");
-double tfac=m_RingDiameter*m_RingDiameter/4-r.Mag2();
-double rfac= m_AxialDistance/m_RingDiameter;
-for (int i=0;i<nmax;i++){
-
- double phi= m_rnd->Rndm()*TMath::Pi();
- TVector3 s(1,0,0);
- s.SetPhi(phi);
- double sign = (m_rnd->Rndm()>0.5)? 1 : 0;
- double theta = TMath::ACos(m_rnd->Rndm()*rfac);
- theta+=sign*TMath::Pi();
-
- s.SetTheta(theta);
- double t=r*s;
- TVector3 rx=r-t*s;
-
- double d=TMath::Sqrt(t*t+tfac);
-
- TVector3 r1=rx+d*s;
- TVector3 r2=rx-d*s;
-
- //r1=Hits2Digits(r1);
- //r2=Hits2Digits(r2);
-
- AddEvent( r1 , r2);
- if (h!=NULL){
- TVector3 s1=r2-r1;
- double s1len= s1.Mag();
- int niter=int (100*s1len/m_RingDiameter);
- for (int j=0;j<niter;j++){
- r2=r1+m_rnd->Rndm()*s1;
- if (h2d) ((TH2F *) h)->Fill(r2.X(),r2.Y());
- else ((TH3F *) h)->Fill(r2.X(),r2.Y(),r2.Z());
- }
- }
-}
-return 0;
-}
-
-int PETProjDataMgr::FwdProject(TH2F *img, TH2F *h){
-
-for (int i=0;i<img->GetNbinsX();i++) {
- double x_=img->GetXaxis()->GetBinCenter( i+1 );
- for (int j=0;j<img->GetNbinsY();j++) {
- double y_=img->GetYaxis()->GetBinCenter( j+1 );
- double density= img->GetBinContent(i+1,j+1);
- if (density>0) FwdProject(x_,y_,m_AxialDistance*(m_rnd->Rndm()-0.5), density,h);
- }
-}
-return 0;
-}
-
-int PETProjDataMgr::FwdProject(TH3F *img, TH3F *h){
-
-for (int i=0;i<img->GetNbinsX();i++) {
- double x_=img->GetXaxis()->GetBinCenter( i+1 );
- for (int j=0;j<img->GetNbinsY();j++) {
- double y_=img->GetYaxis()->GetBinCenter( j+1 );
- for (int k=0;k<img->GetNbinsZ();k++) {
- double z_=img->GetZaxis()->GetBinCenter( k+1 );
- double density= img->GetBinContent(i+1,j+1,k+1);
- if (density>0) FwdProject(x_,y_,z_, density,h);
- }
- }
-}
-return 0;
-}
-
-TH2F *PETProjDataMgr::Phantom(int kaj){
-TH2F *img= new TH2F("img","Original Image",100,-50,50,100,-50,50);
-
-// izberi sliko 0: kroglice, 1: point source 2: central ball
-switch (kaj){
-
-case 0:
-for (int i=0;i<img->GetNbinsX();i++) {
- for (int j=0;j<img->GetNbinsY();j++) {
- double x_=img->GetXaxis()->GetBinCenter( i+1 );
- double y_=img->GetYaxis()->GetBinCenter( j+1 );
- double density=1000;
- if ((x_*x_+y_*y_)<6) img->SetBinContent(i+1,j+1,density);
-
- density=500; if ((x_-25)*(x_-25)+y_*y_<12) img->SetBinContent(i+1,j+1,density);
- density=2000; if ((y_-25)*(y_-25)+x_*x_<2) img->SetBinContent(i+1,j+1,density);
- }
-}
-break;
-
-case 2:
-for (int i=0;i<img->GetNbinsX();i++) {
- for (int j=0;j<img->GetNbinsY();j++) {
- double x_=img->GetXaxis()->GetBinCenter( i+1 );
- double y_=img->GetYaxis()->GetBinCenter( j+1 );
- double density=1000;
- if ((x_*x_+y_*y_)<12.5) img->SetBinContent(i+1,j+1,density);
- }
-}
-break;
-
-case 1:
-img->Fill(25,25,10000);
-break;
-
-}
-
-return img;
-}
Index: petrecofmf/pedestals_new.dat
===================================================================
--- petrecofmf/pedestals_new.dat (revision 112)
+++ petrecofmf/pedestals_new.dat (nonexistent)
@@ -1,64 +0,0 @@
-0 339.13
-1 354.26
-2 327.546
-3 282.69
-4 312.863
-5 365.795
-6 317.055
-7 316.198
-8 348.804
-9 337.947
-10 385.535
-11 341.414
-12 351.877
-13 346.92
-14 339.204
-15 334.106
-16 331.704
-17 333.163
-18 302.789
-19 319.067
-20 325.209
-21 348.417
-22 337.329
-23 326.671
-24 351.699
-25 348.799
-26 308.549
-27 286.519
-28 332.446
-29 324.031
-30 351.92
-31 342.048
-32 359.25
-33 325.419
-34 334.892
-35 322.25
-36 293.95
-37 354.971
-38 304.267
-39 370.495
-40 374.192
-41 347.656
-42 340.723
-43 311.266
-44 301.842
-45 320.301
-46 387.6
-47 318.137
-48 321.878
-49 348.859
-50 342.913
-51 332.14
-52 346.995
-53 350.787
-54 337.321
-55 329.6
-56 346.56
-57 303.365
-58 298.587
-59 302.891
-60 306.159
-61 293.716
-62 324.878
-63 341.973
Index: petrecofmf/config.xml
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
/petrecofmf/config.xml
Property changes:
Deleted: svn:mime-type
## -1 +0,0 ##
-application/xml
\ No newline at end of property
Index: petrecofmf/ph511.dat
===================================================================
--- petrecofmf/ph511.dat (revision 112)
+++ petrecofmf/ph511.dat (nonexistent)
@@ -1,84 +0,0 @@
-0 1700
-1 1900
-2 1680
-3 1700
-4 1870
-5 1680
-6 1430
-7 1460
-8 1680
-9 1700
-10 1713
-11 1650
-12 1630
-13 1830
-14 1750
-15 1740
-16 2440
-17 2320
-18 2640
-19 2500
-20 2790
-21 2600
-22 2900
-23 2820
-24 2660
-25 2470
-26 2310
-27 2330
-28 2817
-29 2860
-30 2140
-31 2140
-32 2120
-33 2410
-34 1990
-35 2480
-36 2200
-37 2100
-38 2040
-39 2060
-40 2090
-41 1970
-42 2032
-43 1970
-44 2054
-45 2120
-46 2076
-47 1920
-48 2206
-49 2420
-50 1930
-51 2290
-52 2400
-53 2425
-54 1950
-55 2400
-56 2076
-57 2010
-58 1792
-59 1830
-60 2000
-61 2032
-62 1640
-63 1614
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: petrecofmf/sinoread.cxx
===================================================================
--- petrecofmf/sinoread.cxx (revision 112)
+++ petrecofmf/sinoread.cxx (nonexistent)
@@ -1,152 +0,0 @@
-#include <stdio.h>
-#include "TH2F.h"
-#include "TCanvas.h"
-#include "TMath.h"
-#include "TStyle.h"
-#include "TColor.h"
-
-void mypalette(){
-
- Int_t MyPalette[1000];
- double stops[] = {0.00, 0.34, 0.61, 0.84, 1.00};
- double red[] = {1.00, 0.84, 0.61, 0.34, 0.00};
- double green[] = {1.00, 0.84, 0.61, 0.34, 0.00};
- double blue[] = {1.00, 0.84, 0.61, 0.34, 0.00};
-
-
- Int_t gs = TColor::CreateGradientColorTable(5, stops, red, green, blue, 999);
- //for (int i=0;i<999;i++) MyPalette[i] = gs+998-i;
-
- for (int i=0;i<999;i++) MyPalette[i] = gs+i;
- gStyle->SetPalette(999, MyPalette);
-}
-
-int imageread(char *fname,int nchx=28, int nchy=56, int nchz=0, int format=0, int size=0, char *opt=""){
-FILE *fp= fopen(fname,"read");
-const int maxbuf=360000;
-float *fbuf= new float[maxbuf];
-unsigned short*ibuf=((unsigned short*) &fbuf[0]);
-if (size==0) {
- if (format) size=sizeof(float);
- else size=sizeof(unsigned short);
-}
-int ncount=0;
-
-if(fp){
-
-gStyle->SetPalette(100);
-mypalette();
-TCanvas *c=new TCanvas("c","c",1500,1500);
-int sizx=sqrt(nchz);
-c->Divide(sizx,TMath::Ceil(nchz/double(sizx)) );
-char hname[128];
-
-while (!feof(fp)){
-
-
-int nb=fread(fbuf,size,nchx*nchy,fp);
-printf("nb=%d\n",nb);
-if (nb>0){
- sprintf(hname,"slice %d",ncount);
- TH2F *h=new TH2F("h",hname,nchx,-0.5,nchx-0.5,nchy,-0.5,nchy-0.5);
- float val=0;
- for (int i=0 ; i<nb;i++){
- //h->SetBinContent(i%nchx+1, i/nchx+1, fbuf[i]);
- if (format) val=fbuf[i]; else val=float(ibuf[i]);
- h->SetBinContent(i/nchy+1,i%nchy+1, val);
- }
- c->cd(ncount+1);
- h->SetMinimum(0);
- h->Draw(opt);
-}
-
-ncount++;
-}
-fclose(fp);
-c->Modified();
-c->Update();
-}
-
-delete fbuf;
-return 0;
-}
-
-int interfileread(char *fname, char *opt=""){
- FILE *fp=fopen(fname,"r");
- if (!fp) {
- printf("Cannot open %s\n",fname);
- return 0;
- }
- int j=0;
- int ndim=400;
- char *line= new char[ndim];
- int indx[4]={-1,-1,-1,-1};
- int dim[5];
- char *cmd,*val,*token;
- char fimage[100];
- while (fgets(line,ndim,fp)!=NULL) {
- char *v= strstr(line,":=");
- if (v!=NULL) {
- cmd=line;
- v[0]=0;
- val = &v[2];
- val[strlen(val)-1]=0;
- printf("%d #%s# %s\n",j++, val, cmd);
- int axis=-1;
- if (strstr(cmd,"matrix axis label")!= NULL){
-
- sscanf(cmd,"matrix axis label [%d]", &axis);
- if(axis>0 && axis<5){
- if (strstr(val," x")!= NULL) indx[axis-1]=1;
- if (strstr(val," y")!= NULL) indx[axis-1]=0;
- if (strstr(val," z")!= NULL) indx[axis-1]=3;
-
- if (strstr(val,"view") != NULL) indx[axis-1]=0;
- if (strstr(val,"tangential")!= NULL) indx[axis-1]=1;
- if (strstr(val,"segment") != NULL) indx[axis-1]=2;
- if (strstr(val,"axial") != NULL) indx[axis-1]=3;
- printf("-->> %s axis=%d indx=%d\n",val, axis,indx[axis-1]);
-
- }
- }
- if (strstr(cmd,"!matrix size")!= NULL){
- sscanf(cmd,"!matrix size [%d]", &axis);
- if(axis>0 && axis<5){
- if (indx[axis-1]<3){ dim[indx[axis-1]]=atoi(val);}
- else {
- int sum=0;
- token = strtok(val,"{,");
- sum+=atoi(token);
- while ( (token=strtok(NULL,"{,")) ) sum+=atoi(token);
- printf("total dim=%d\n",sum);
- dim[3]=sum;
- };
- }
- printf("axis=%d indx=%d\n",axis,indx[axis-1]);
- }
- if (strstr(cmd,"name of data file")!= NULL){
- sscanf(val,"%s",fimage);
- printf("image file=%s\n",fimage);
- }
- if (strstr(cmd,"!number of bytes per pixel")!= NULL){
- dim[4]=atoi(val);
- printf("bytes per pixel=%d\n",dim[4]);
- }
- if (strstr(cmd,"!number format")!= NULL){
- if (strstr(val,"float")!=NULL){ dim[5]=1; };
- if (strstr(val,"integer")!=NULL){ dim[5]=0; };
- printf("number format=%d\n",dim[5]);
- }
-
- }
- }
- for (int i=0;i<6;i++) printf ("dim[%d]=%d\n", i,dim[i]);
- fclose(fp);
- imageread(fimage,dim[0], dim[1], dim[3], dim[5], dim[4], opt);
- return 0;
-}
-
-
-int sinoread(char *fname, char *opt=""){
- return interfileread(fname,opt);
-}
Index: petrecofmf/FBP2D.sh
===================================================================
--- petrecofmf/FBP2D.sh (revision 112)
+++ petrecofmf/FBP2D.sh (nonexistent)
@@ -1,58 +0,0 @@
-#!/bin/bash
-#
-export STIRPATH=~rok/pet/reco/STIR2.1/STIR
-#
-export PATH=$PATH:$STIRPATH/bin
-EXE=FBP2D
-IMAGEVIEW=manip_image
-HVNAME=$2.hv
-
-TEMPFILE=`tempfile`
-rm -f $TEMPFILE
-echo -n "fbp2dparameters :=
-
-input file := $1
-output filename prefix := $2
-
-;input file := sino3D.hs
-;output filename prefix := output
-
-
-; output image parameters
-; zoom defaults to 1
-zoom := 1
-; image size defaults to whole FOV
-; ;;;;;;;; zakomentirano 8.1.2013 xy output image size (in pixels) := 100
-xy output image size (in pixels) := 200
-
-; can be used to call SSRB first
-; default means: call SSRB only if no axial compression is already present
-;num segments to combine with ssrb := -1
-
-; filter parameters, default to pure ramp
-alpha parameter for ramp filter := 0.5
-cut-off for ramp filter (in cycles) := 0.3
-
-; allow less padding. DO NOT USE
-; (unless you're sure that the object occupies only half the FOV)
-;Transaxial extension for FFT:=1
-
-; back projector that could be used (defaults to interpolating backprojector)
-; Back projector type:= some type
-
-; display data during processing for debugging purposes
-Display level := 1
-end :=
-" >> $TEMPFILE
-
-echo $EXE $TEMPFILE
-$EXE $TEMPFILE
-
-
-$IMAGEVIEW $HVNAME
-root src/sinoread.cxx+\(\"$HVNAME\",\"colz\"\)
-
-echo "za ponovni ogled slike"
-echo $IMAGEVIEW $HVNAME
-echo ali
-echo root src/sinoread.cxx+\(\"$HVNAME\",\"colz\"\)
/petrecofmf/FBP2D.sh
Property changes:
Deleted: svn:executable
Index: praktikum/petrecofmf/Makefile
===================================================================
--- praktikum/petrecofmf/Makefile (nonexistent)
+++ praktikum/petrecofmf/Makefile (revision 113)
@@ -0,0 +1,60 @@
+
+ROOTINC=$(shell root-config --incdir )
+ROOTLIB=$(shell root-config --libs )
+ROOTGLIBS=$(shell root-config --libs --glibs )
+
+LIBS=$(ROOTLIB) -L./ -lm
+
+
+XMLCFLAGS =$(shell xml2-config --cflags )
+XMLLIBS =$(shell xml2-config --libs )
+
+
+INC=-I. -I$(ROOTINC) $(XMLCFLAGS)
+LIBS=$(ROOTGLIBS) $(XMLLIBS)
+
+SRC = .
+INC1 = -I. -I../lib -I/usr/include
+DBG =
+CFLAGS = $(DBG) $(INC1) -Wall -g
+
+
+
+
+BIN= bin
+LIB=lib
+
+TARGET = $(BIN)/petreco
+FILES = $(SRC)/readdata.C $(SRC)/PETProjDataMgr.C
+
+LIBFILE = $(LIB)/libfmfpetreco.a
+
+
+
+OBJ_FILES = $(SRC)/readdata.o
+
+$(TARGET): $(FILES)
+ $(CXX) $(INC) -DMAIN $(FILES) $(CFLAGS) -o $(TARGET) $(LIBS) -lstdc++
+
+
+
+.C.o:
+ $(CXX) $(INC) -c $<
+ ar r $(LIBFILE) $@
+
+.c.o:
+ $(CXX) -c $<
+ ar r $(LIBFILE) $@
+
+
+xpath2: $(SRC)/xpath2.c
+ $(CXX) -o $(BIN)/xpath2 `xml2-config --cflags` $(SRC)/xpath2.c `xml2-config --libs`
+
+
+
+clean:
+ rm Dict.C $(OBJ_FILES) $(TARGET) $(LIBFILE)
+
+
+tgz:
+ tar czvf fmfpet.tgz Makefile README $(SRC)/*.C $(SRC)/*.h $(SRC)/*.c $(SRC)/*.cxx $(SRC)/*.sh *.xml *.dat *.map *.root
Index: praktikum/petrecofmf/README
===================================================================
--- praktikum/petrecofmf/README (nonexistent)
+++ praktikum/petrecofmf/README (revision 113)
@@ -0,0 +1,2 @@
+http://www-f9.ijs.si/wiki/Main/PetRekonstrukcija
+
Index: praktikum/petrecofmf/sumpedestals_zero.dat
===================================================================
--- praktikum/petrecofmf/sumpedestals_zero.dat (nonexistent)
+++ praktikum/petrecofmf/sumpedestals_zero.dat (revision 113)
@@ -0,0 +1,4 @@
+0 6000
+1 5264.29
+2 5394.42
+3 5261.11
Index: praktikum/petrecofmf/m16.map
===================================================================
--- praktikum/petrecofmf/m16.map (nonexistent)
+++ praktikum/petrecofmf/m16.map (revision 113)
@@ -0,0 +1,17 @@
+#CH IX IY
+15 0 0
+14 0 1
+13 1 0
+12 1 1
+11 0 3
+10 0 2
+9 1 3
+8 1 2
+7 2 3
+6 2 2
+5 3 3
+4 3 2
+3 2 0
+2 2 1
+1 3 0
+0 3 1
Index: praktikum/petrecofmf/sumpedestals.dat
===================================================================
--- praktikum/petrecofmf/sumpedestals.dat (nonexistent)
+++ praktikum/petrecofmf/sumpedestals.dat (revision 113)
@@ -0,0 +1,4 @@
+0 5413.22
+1 5264.29
+2 5394.42
+3 5261.11
Index: praktikum/petrecofmf/modules.map
===================================================================
--- praktikum/petrecofmf/modules.map (nonexistent)
+++ praktikum/petrecofmf/modules.map (revision 113)
@@ -0,0 +1,5 @@
+#module ID (connector) r(mm) angle(deg)
+0 61 0
+1 61 22.5
+2 61 180
+3 61 202.5
Index: praktikum/petrecofmf/pedestals_zero.dat
===================================================================
--- praktikum/petrecofmf/pedestals_zero.dat (nonexistent)
+++ praktikum/petrecofmf/pedestals_zero.dat (revision 113)
@@ -0,0 +1,64 @@
+0 570
+1 570
+2 570
+3 570
+4 570
+5 570
+6 570
+7 570
+8 570
+9 570
+10 570
+11 570
+12 570
+13 570
+14 570
+15 570
+16 331.704
+17 333.163
+18 302.789
+19 319.067
+20 325.209
+21 348.417
+22 337.329
+23 326.671
+24 351.699
+25 348.799
+26 308.549
+27 286.519
+28 332.446
+29 324.031
+30 351.92
+31 342.048
+32 359.25
+33 325.419
+34 334.892
+35 322.25
+36 293.95
+37 354.971
+38 304.267
+39 370.495
+40 374.192
+41 347.656
+42 340.723
+43 311.266
+44 301.842
+45 320.301
+46 387.6
+47 318.137
+48 321.878
+49 348.859
+50 342.913
+51 332.14
+52 346.995
+53 350.787
+54 337.321
+55 329.6
+56 346.56
+57 303.365
+58 298.587
+59 302.891
+60 306.159
+61 293.716
+62 324.878
+63 341.973
Index: praktikum/petrecofmf/pedestals.dat
===================================================================
--- praktikum/petrecofmf/pedestals.dat (nonexistent)
+++ praktikum/petrecofmf/pedestals.dat (revision 113)
@@ -0,0 +1,64 @@
+0 339.764
+1 356.526
+2 328.481
+3 283.352
+4 307.915
+5 365.996
+6 316.9
+7 316.175
+8 348.054
+9 338.349
+10 385.901
+11 342.318
+12 353.552
+13 348.903
+14 340.133
+15 335.819
+16 337.356
+17 334.646
+18 322.328
+19 320.371
+20 304.205
+21 347.213
+22 336.6
+23 325.069
+24 351.203
+25 347.073
+26 307.865
+27 285.498
+28 332.387
+29 325.668
+30 352.567
+31 343.887
+32 360.176
+33 328.103
+34 337.231
+35 324.995
+36 294.699
+37 354.216
+38 304.069
+39 369.572
+40 375.065
+41 347.03
+42 341.669
+43 311.298
+44 338.887
+45 322.734
+46 388.779
+47 320.116
+48 323.118
+49 350.691
+50 344.091
+51 334.575
+52 349.252
+53 351.313
+54 338.096
+55 330.029
+56 348.157
+57 304.036
+58 300.07
+59 304.476
+60 307.644
+61 296.693
+62 327.384
+63 345.163
Index: praktikum/petrecofmf/fotovrh.dat
===================================================================
--- praktikum/petrecofmf/fotovrh.dat (nonexistent)
+++ praktikum/petrecofmf/fotovrh.dat (revision 113)
@@ -0,0 +1,64 @@
+0 3196.160841
+1 1711.854801
+2 1635.564402
+3 1597.419203
+4 1711.854801
+5 1584.704136
+6 1355.832941
+7 1432.123339
+8 1533.843871
+9 1648.279469
+10 1597.419203
+11 1521.128804
+12 1470.268539
+13 1686.424668
+14 1622.849336
+15 1660.994535
+16 2522.262320
+17 2220.457458
+18 2423.898521
+19 2334.893056
+20 2551.049185
+21 2449.328654
+22 2703.629982
+23 2614.624517
+24 2487.473853
+25 2347.608122
+26 2156.882126
+27 2233.172524
+28 2589.194384
+29 2233.172524
+30 1864.435598
+31 1915.295864
+32 2852.854047
+33 1915.295864
+34 1915.295864
+35 2156.882126
+36 1991.586262
+37 1966.156129
+38 1915.295864
+39 1902.580797
+40 1940.725996
+41 1877.150664
+42 1940.725996
+43 1877.150664
+44 1928.010930
+45 1940.725996
+46 1991.586262
+47 1928.010930
+48 2852.854047
+49 2080.591727
+50 1813.575332
+51 1978.871196
+52 2245.887591
+53 2080.591727
+54 1813.575332
+55 1940.725996
+56 1902.580797
+57 1889.865731
+58 1673.709601
+59 1699.139734
+60 1813.575332
+61 1877.150664
+62 1508.413738
+63 1686.424668
Index: praktikum/petrecofmf/calibration_all.root
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream
/praktikum/petrecofmf/calibration_all.root
Property changes:
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property
Index: praktikum/petrecofmf/config_zero.xml
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
/praktikum/petrecofmf/config_zero.xml
Property changes:
Added: svn:mime-type
## -0,0 +1 ##
+application/xml
\ No newline at end of property
Index: praktikum/petrecofmf/config.xml
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
/praktikum/petrecofmf/config.xml
Property changes:
Added: svn:mime-type
## -0,0 +1 ##
+application/xml
\ No newline at end of property
Index: praktikum/petrecofmf/pedestals_new.dat
===================================================================
--- praktikum/petrecofmf/pedestals_new.dat (nonexistent)
+++ praktikum/petrecofmf/pedestals_new.dat (revision 113)
@@ -0,0 +1,64 @@
+0 339.13
+1 354.26
+2 327.546
+3 282.69
+4 312.863
+5 365.795
+6 317.055
+7 316.198
+8 348.804
+9 337.947
+10 385.535
+11 341.414
+12 351.877
+13 346.92
+14 339.204
+15 334.106
+16 331.704
+17 333.163
+18 302.789
+19 319.067
+20 325.209
+21 348.417
+22 337.329
+23 326.671
+24 351.699
+25 348.799
+26 308.549
+27 286.519
+28 332.446
+29 324.031
+30 351.92
+31 342.048
+32 359.25
+33 325.419
+34 334.892
+35 322.25
+36 293.95
+37 354.971
+38 304.267
+39 370.495
+40 374.192
+41 347.656
+42 340.723
+43 311.266
+44 301.842
+45 320.301
+46 387.6
+47 318.137
+48 321.878
+49 348.859
+50 342.913
+51 332.14
+52 346.995
+53 350.787
+54 337.321
+55 329.6
+56 346.56
+57 303.365
+58 298.587
+59 302.891
+60 306.159
+61 293.716
+62 324.878
+63 341.973
Index: praktikum/petrecofmf/ph511.dat
===================================================================
--- praktikum/petrecofmf/ph511.dat (nonexistent)
+++ praktikum/petrecofmf/ph511.dat (revision 113)
@@ -0,0 +1,84 @@
+0 1700
+1 1900
+2 1680
+3 1700
+4 1870
+5 1680
+6 1430
+7 1460
+8 1680
+9 1700
+10 1713
+11 1650
+12 1630
+13 1830
+14 1750
+15 1740
+16 2440
+17 2320
+18 2640
+19 2500
+20 2790
+21 2600
+22 2900
+23 2820
+24 2660
+25 2470
+26 2310
+27 2330
+28 2817
+29 2860
+30 2140
+31 2140
+32 2120
+33 2410
+34 1990
+35 2480
+36 2200
+37 2100
+38 2040
+39 2060
+40 2090
+41 1970
+42 2032
+43 1970
+44 2054
+45 2120
+46 2076
+47 1920
+48 2206
+49 2420
+50 1930
+51 2290
+52 2400
+53 2425
+54 1950
+55 2400
+56 2076
+57 2010
+58 1792
+59 1830
+60 2000
+61 2032
+62 1640
+63 1614
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: praktikum/petrecofmf/m16_reverse.map
===================================================================
--- praktikum/petrecofmf/m16_reverse.map (nonexistent)
+++ praktikum/petrecofmf/m16_reverse.map (revision 113)
@@ -0,0 +1,17 @@
+#CH IX IY
+0 0 0
+1 0 1
+2 1 0
+3 1 1
+4 0 3
+5 0 2
+6 1 3
+7 1 2
+8 2 3
+9 2 2
+10 3 3
+11 3 2
+12 2 0
+13 2 1
+14 3 0
+15 3 1
Index: praktikum/petrecofmf/plot.cxx
===================================================================
--- praktikum/petrecofmf/plot.cxx (nonexistent)
+++ praktikum/petrecofmf/plot.cxx (revision 113)
@@ -0,0 +1,186 @@
+//---------------------------------------------------
+int plot2d(char *opt="colz", float max=0){
+char hname[50];
+ char hn[50];
+ TCanvas *c;
+ TH1I *gICh = (TH1I * ) gDirectory->Get("mapping");
+ for (int i=0;i<4;i++) {
+ sprintf(hname,"cz%d",i);
+ sprintf(hn,"Corr Pmt %d",i);
+ c =new TCanvas(hname,hn, 1000,1000);
+ c->Divide(4,4);
+ for (int j=0;j<16;j++){
+ c->cd(gICh->GetBinContent(j+1)+1);
+ sprintf(hn,"corrch%d",i*16+j);
+ TH1F*h= (TH1F * ) gDirectory->Get(hn);
+ if (max>0) h->SetMaximum(max);
+ h->Draw(opt);
+ }
+ c->Modified();
+ c->Update();
+ }
+return 0;
+}
+
+//---------------------------------------------------
+int plot( char *fname=NULL, float max=0){
+ if (fname!=NULL){
+ cout << "Rootname:" << fname << endl;
+ delete gROOT->GetListOfFiles()->FindObject(fname); // clear memory of file name
+ if( gSystem->AccessPathName(fname) ) {
+ cout << endl << "File: " << fname << " not there!!!" << endl << endl;
+ return 0;
+ }
+ TFile *rootfile = new TFile(fname);
+ }
+ //if (buf) delete buf;
+ char hn[255];
+
+ TCanvas *c1;
+ for (int i=0;i<4;i++){
+ sprintf(hn,"csumadc%d",i);
+ c1=new TCanvas(hn,hn,1200,1200);
+ c1->Divide(9,9,0,0);
+ sprintf(hn,"sumadc%d",i);
+ TH2F *h = ((TH2F * )gDirectory->Get(hn));
+ for (int j=0;j<81;j++){
+ c1->cd(1+j);
+ TH1D *h0 = h->ProjectionX(" ",1+j,1+j);
+ h0->DrawCopy();
+ }
+ c1->Modified();
+ c1->Update();
+ }
+
+
+
+ TCanvas *c2=new TCanvas("c2","c2",800,800);
+ c2->Divide(2,2);
+ for (int i=0;i<4;i++){
+ c2->cd(1+i);
+ sprintf(hn,"pmt1%d",i);
+ ((TH1F * ) gDirectory->Get(hn))->Draw("colz");
+ }
+ c2->Modified();
+ c2->Update();
+
+ TCanvas *c6=new TCanvas("c6","c6",800,800);
+ c6->Divide(2,2);
+ for (int i=0;i<4;i++){
+ c6->cd(1+i);
+ sprintf(hn,"pmt3%d",i);
+ ((TH1F * ) gDirectory->Get(hn))->Draw("colz");
+ }
+ c6->Modified();
+ c6->Update();
+
+ TCanvas *c5=new TCanvas("c5","c5",800,800);
+
+
+ ((TH1F * ) gDirectory->Get("globalxy"))->Draw("colz");
+
+ c5->Modified();
+ c5->Update();
+
+
+ TCanvas *c3=new TCanvas("c3","c3",800,800);
+ c3->Divide(2,2);
+ for (int i=0;i<4;i++){
+ c3->cd(1+i);
+ sprintf(hn,"pmt2%d",i);
+ ((TH1F * ) gDirectory->Get(hn))->Draw("colz");
+ }
+
+
+ c3->Modified();
+ c3->Update();
+
+ TCanvas *c4=new TCanvas("c4","c4",800,800);
+ c4->Divide(4,4);
+ for (int i=0;i<16;i++){
+ c4->cd(1+i)->SetLogy();
+ sprintf(hn,"ach%d",i);
+ ((TH1F * ) gDirectory->Get(hn))->Draw();
+ }
+ c4->Modified();
+ c4->Update();
+ plot2d("colz",max);
+ return 0;
+}
+
+//---------------------------------------------------
+int plotadc(){
+ char hn[255];
+ TCanvas *c;
+ char hname[50];
+ char hn[50];
+ TH1I *gICh = (TH1I * ) gDirectory->Get("mapping");
+ for (int i=0;i<4;i++) {
+ sprintf(hname,"cx%d",i);
+ sprintf(hn,"Pmt %d",i);
+ c =new TCanvas(hname,hn, 1000,1000);
+ c->Divide(4,4);
+ for (int j=0;j<16;j++){
+ c->cd(gICh->GetBinContent(j+1)+1)->SetLogy();
+ sprintf(hn,"ach%d",i*16+j);
+ ((TH1F * ) gDirectory->Get(hn))->Draw();
+ }
+ c->Modified();
+ c->Update();
+ }
+
+ return 0;
+}
+
+//---------------------------------------------------
+int plotsingles(char *opt="colz", float max=0){
+ char hname[50];
+ char hn[50];
+ TCanvas *c;
+ TH1I *gICh = (TH1I * ) gDirectory->Get("mapping");
+ for (int i=0;i<4;i++) {
+ sprintf(hname,"cy%d",i);
+ sprintf(hn,"Single Pmt %d",i);
+ c =new TCanvas(hname,hn, 1000,1000);
+ c->Divide(4,4);
+ for (int j=0;j<16;j++){
+ c->cd(gICh->GetBinContent(j+1)+1);
+ sprintf(hn,"singlech%d",i*16+j);
+ TH1F*h= (TH1F * ) gDirectory->Get(hn);
+ if (max>0) gHistoSingle[i*16+j]->SetMaximum(max);
+ h->Draw(opt);
+ }
+ c->Modified();
+ c->Update();
+ }
+
+ return 0;
+}
+
+
+
+
+//---------------------------------------------------
+int plotcut(char *opt="", float max=0){
+char hname[50];
+ char hn[50];
+ TCanvas *c;
+ TH1I *gICh = (TH1I * ) gDirectory->Get("mapping");
+ for (int i=0;i<4;i++) {
+ sprintf(hname,"cz%d",i);
+ sprintf(hn,"Cut Pmt %d",i);
+ c =new TCanvas(hname,hn, 1000,1000);
+ c->Divide(4,4);
+ for (int j=0;j<16;j++){
+ c->cd(gICh->GetBinContent(j+1)+1)->SetLogy();
+ sprintf(hn,"corrch%d",i);
+ TH1F*h= (TH1F * ) gDirectory->Get(hn);
+ if (max>0) h->SetMaximum(max);
+ h->Draw(opt);
+ }
+ c->Modified();
+ c->Update();
+ }
+return 0;
+}
+//---------------------------------------------------
Index: praktikum/petrecofmf/PETProjDataMgr.C
===================================================================
--- praktikum/petrecofmf/PETProjDataMgr.C (nonexistent)
+++ praktikum/petrecofmf/PETProjDataMgr.C (revision 113)
@@ -0,0 +1,612 @@
+//
+// ../bin/FBP2D FBP2D.par
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The GAMOS software is copyright of the Copyright Holders of *
+// * the GAMOS Collaboration. It is provided under the terms and *
+// * conditions of the GAMOS Software License, included in the file *
+// * LICENSE and available at http://fismed.ciemat.es/GAMOS/license .*
+// * These include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GAMOS collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the GAMOS Software license. *
+// ********************************************************************
+//
+#include "PETProjDataMgr.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TMath.h"
+#include "TRandom3.h"
+
+#include <iomanip>
+#include <iostream>
+#include <cstdlib>
+
+using namespace std;
+
+//----------------------------------------------------------------------
+PETProjDataMgr* PETProjDataMgr::m_Instance = 0;
+
+//----------------------------------------------------------------------
+PETProjDataMgr* PETProjDataMgr::GetInstance()
+{
+ if( !m_Instance ) {
+ m_Instance = new PETProjDataMgr;
+ }
+
+ return m_Instance;
+
+}
+
+//-----------------------------------------------------------------------
+PETProjDataMgr::PETProjDataMgr()
+{
+
+ /*
+// ARGUMENTS:
+" cout << " -------------------------- \n"
+" Arguments convention: \n"
+" -a Axial FOV (mm), <theDist_axial=100.0> \n"
+" -d Diameter Transaxial FOV (mm), <m_RingDiameter=300.0> \n"
+" -i Type of the input file (by default: 0 = Arce_binary), <typeINfile=0> \n"
+" \n"
+" -m Maximum ring difference (by default: -1 = m_NOfPlanes), <m_MaxRingDifferenceiff> \n"
+" -n Name of output file, <m_Filename> \n"
+" -p Axial number of planes, <m_NOfPlanes> \n"
+" -r Number of bins, \"distancias\", <m_NOfBins> \n"
+ // -s Span TO DO:span !!!!!
+" -t Number of angular views, \"direcciones\", <m_NOfAngles> \n"
+" -v Verbosity (by default: 0=silent, 3=debug), <verbos> \n"
+" -x Maximum number of coincidences to be stored (by default: -1 = no limit), <Max_Coinci> \n"
+" -o Output type (by default: 0 = mcc Interfile, 1 = STIR Interfile), <OutType> \n"
+
+" \n"
+" PET Reconstruction. CIEMAT 2009-11 \n"
+" mario.canadas@ciemat.es \n"
+ " -------------------------- \n";
+ */
+ m_rnd= new TRandom3();
+ m_AxialDistance = (9-1)*2.25; // Axial pixel dimension*NOfPlanes
+ m_RingDiameter = 120.0; // notranji premer peta
+ m_NOfPlanes = 9; // stevilo ravnin
+ m_NOfBins = 128; // stevilo binov v razdalji
+ m_nch = 128; // stevilo padov okoli in okoli
+ m_NOfAngles = TMath::Abs(m_nch)/2; // stevilo kotov = stevilo padov okoli in okoli /2
+
+ m_MaxRingDifference = -1;
+ //m_MaxRingDifference = 3; // najvecja razdalja med padi
+ // toDo: theSpan = int(GmParameterMgr::GetInstance()->GetNumericValue("PET:ProjData:Span",1));
+
+
+ m_OutFormat = 1; // 1.. projections 0 .. image
+ m_Filename = "sino3D";
+ m_Debug=1;
+
+ if (m_MaxRingDifference==-1) m_MaxRingDifference=m_NOfPlanes-1;
+
+
+ m_TotalAxialPlanes=m_NOfPlanes*m_NOfPlanes;
+ if (m_OutFormat==1) m_TotalAxialPlanes= (2*m_NOfPlanes-1 - m_MaxRingDifference)*m_MaxRingDifference + m_NOfPlanes; // total number of Axial planes (segments*planes) in STIR format
+
+ /*--- Initialize sino3D ---*/
+ m_projections = new SINO_TYPE**[m_NOfBins];
+ for(int i=0;i<m_NOfBins;i++){
+ m_projections[i] = new SINO_TYPE*[m_NOfAngles];
+ for(int j=0;j<m_NOfAngles;j++){
+ m_projections[i][j] = new SINO_TYPE[m_TotalAxialPlanes]; /// ! If m_OutFormat==1 (STIR output):Matrix size depends on the MAX_Ring_Difference
+ for(int k=0;k<m_TotalAxialPlanes;k++){
+ m_projections[i][j][k]=0;
+ }
+ }
+ }
+
+ m_TotalProjectionCoincidences=0;
+ m_TotalCoincidences=0;
+ //OutputType = "pet";
+}
+
+//-----------------------------------------------------------------------
+void PETProjDataMgr::SetProjection( int axialplane, TH2F * h)
+{
+for(int i=0;i<m_NOfBins;i++){
+ for(int j=0;j<m_NOfAngles;j++){
+ m_projections[i][j][axialplane]=h->GetBinContent(i+1,j+1);
+ }
+}
+
+}
+//-----------------------------------------------------------------------
+
+//-----------------------------------------------------------
+// from Gate
+//------------------------------------------------------------
+double ComputeSinogramS(double X1, double Y1, double X2, double Y2)
+{
+ double s;
+
+ double denom = (Y1-Y2) * (Y1-Y2) + (X2-X1) * (X2-X1);
+
+ if (denom!=0.) {
+ denom = sqrt(denom);
+ s = ( X1 * (Y2-Y1) + Y1 * (X1-X2) ) / denom;
+ } else {
+ s = 0.;
+ }
+
+ double theta;
+ if ((X1-X2)!=0.) {
+ theta=atan((X1-X2) /(Y1-Y2));
+ } else {
+ theta=3.1416/2.;
+ }
+ if ((theta > 0.) && ((X1-X2) > 0.)) s = -s;
+ if ((theta < 0.) && ((X1-X2) < 0.)) s = -s;
+ if ( theta < 0.) {
+ theta = theta+3.1416;
+ s = -s;
+ }
+ return s;
+}
+
+
+void PETProjDataMgr::AddEvent( const TVector3& pos1 , const TVector3& pos2)
+{
+ int z1_i, z2_i;
+ //for discretization on the crystal: int x1_i, x2_i, y1_i, y2_i;
+
+ double z1_abs=pos1.z()+m_AxialDistance/2;
+ double z2_abs=pos2.z()+m_AxialDistance/2;
+ double a, b, phi, dis;
+ int phi_i, dis_i;
+ int ring_diff;
+
+ double _PI=2*asin(1);
+
+ m_TotalCoincidences++;
+
+ z1_i=(int)(m_NOfPlanes* z1_abs/m_AxialDistance); //round --> m_NOfPlanes+1 ...
+ z2_i=(int)(m_NOfPlanes* z2_abs/m_AxialDistance);
+
+ // control; if z_i out of range: return
+
+ if ( (pos1.x()==pos2.x()) && (pos1.y()==pos2.y()) ) {
+#ifndef GAMOS_NO_VERBOSE
+ if( m_Debug ) {
+ cout << "PETProjDataMgr::AddEvent:WARNING! Event_1 == Event_2 ; x= " << pos2.x() << ", y= " << pos2.y() << ", z= " << pos2.z() << endl;
+ }
+#endif
+ return;
+ }
+
+ if ( (z1_i<0) || (z2_i<0) || (z1_i>= m_NOfPlanes) || (z2_i>= m_NOfPlanes) ) {
+#ifndef GAMOS_NO_VERBOSE
+ if( m_Debug ) {
+ cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Axial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
+ }
+#endif
+ return;
+ }
+
+ ring_diff = (int)fabs(z1_i-z2_i);
+
+ // max ring difference; control:
+ if (ring_diff > m_MaxRingDifference) {
+#ifndef GAMOS_NO_VERBOSE
+ if( m_Debug ) {
+ cout <<"PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Max. Ring Diff.): " << ring_diff << ">" << m_MaxRingDifference << " x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
+ }
+#endif
+ return;
+ }
+
+ a=(double)(pos2.y()- pos1.y());
+ b=(double)(pos2.x()- pos1.x());
+
+ if (a==0.0){
+ phi=_PI*0.5;
+ }
+ else{
+ phi=atan(b/a);
+ }
+
+ if (phi<0) phi = phi +_PI;
+
+ dis=pos1.x()*cos(phi) - pos1.y()*sin(phi);
+ //dis=ComputeSinogramS(pos1.x(), pos1.y(), pos2.x(), pos2.x());
+ // control; transaxial FOV
+ if ( fabs(dis) > m_RingDiameter*0.5 ) {
+#ifndef GAMOS_NO_VERBOSE
+ if( m_Debug ) {
+ cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Transaxial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
+ }
+#endif
+ return;
+ }
+
+ dis = dis + m_RingDiameter*0.5;
+
+ // discret values:
+ phi_i=(int)round( (double)(m_NOfAngles-1)*phi/_PI );
+ dis_i=(int)round( (double)(m_NOfBins-1)*dis/(double)m_RingDiameter );
+
+ if ((phi_i>=m_NOfAngles) || (dis_i>=m_NOfBins)) return; // only possible "=" because 'round' check it..
+
+ // OLD: (SRRB included) sino3D[dis_i][phi_i][ (z1_i+z2_i)+ring_diff*(m_NOfPlanes-1) ]++;
+
+ int Zpos;
+
+ if (m_OutFormat==0) {
+ Zpos = (z1_i*m_NOfPlanes + z2_i);
+ }
+ else{
+
+ if (z1_i>=z2_i) { // SIN Max Ring_Diff: Zpos= ( ((m_NOfPlanes-ring_diff)*(m_NOfPlanes-1-ring_diff))/2 + z2_i );
+
+ Zpos= ( ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff)*(m_MaxRingDifference - ring_diff))/2 + z2_i);
+
+ }else{
+ Zpos= ( (m_TotalAxialPlanes) - ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff +1)*(m_MaxRingDifference - ring_diff +1))/2 + z1_i );
+
+ }
+ }
+
+ m_projections[dis_i][phi_i][ Zpos ]++;
+ m_TotalProjectionCoincidences++;
+
+#ifndef GAMOS_NO_VERBOSE
+ if( m_Debug >1) {
+ cout << "PETProjDataMgr::AddEvent: x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
+ cout << "PETProjDataMgr::AddEvent: Sinogram pos.: distance(s)= " << dis_i << ", angular view(phi)= " << phi_i << " ; Zpos=" << Zpos <<"; Segment (Ring diff.) = " << ring_diff << endl;
+ }
+#endif
+
+
+}
+
+//-----------------------------------------------------------------------
+PETProjDataMgr::~PETProjDataMgr()
+{
+ int i,j;
+
+ for(i=0;i<m_NOfBins;i++){
+ for(j=0;j<m_NOfAngles;j++){
+ free(m_projections[i][j]);
+
+ }
+ free(m_projections[i]);
+
+ }
+ free(m_projections);
+
+}
+
+/* TO DO: call lm_to_sino3D program
+//-----------------------------------------------------------------------
+void PETIOMgr::ReadFile()
+{
+ if( !theFileIn ) OpenFileIn();
+
+ PETOutput po;
+ G4bool bEof;
+ for(;;) {
+ po = ReadEvent( bEof );
+ // theFileIn->read(reinterpret_cast<char *>(&po),sizeof(PetOutput));
+ if(bDumpCout) PrintEvent(" PETOutput: ", po, bCartesian);
+ if( bEof ) break;
+ }
+}
+//-----------------------------------------------------------------------
+PETOutput PETProjDataMgr::ReadEvent( G4bool& bEof )
+{
+ if( theFileIn == 0 ){
+ G4Exception("PETIOMgr::ReadEvent, file not opened, call OpenFileIn() first ");
+ }
+
+ PETOutput po;
+ fread (&po, sizeof(struct PETOutput),1,theFileIn);
+ if ( feof( theFileIn ) ) {
+ bEof = TRUE;
+ } else {
+ bEof = FALSE;
+ }
+
+ return po;
+
+}
+*/
+
+
+//-----------------------------------------------------------------------
+void PETProjDataMgr::WriteInterfile()
+{
+
+ char name_hv[512];
+ char name_v[512];
+
+ if (m_OutFormat==0){
+
+ strcpy(name_hv, m_Filename);
+ strcpy(name_v,m_Filename);
+ strcat(name_hv, ".hv");
+ strcat(name_v, ".v");
+
+ fp=fopen(name_hv, "w");
+
+ fprintf (fp, "!INTERFILE := \n");
+ fprintf (fp, "name of data file := %s\n", name_v);
+ fprintf (fp, "!GENERAL DATA := \n");
+ fprintf (fp, "!GENERAL IMAGE DATA :=\n");
+ fprintf (fp, "!type of data := tomographic\n");
+ fprintf (fp, "!version of keys := 3.3\n");
+ fprintf (fp, "!data offset in bytes := 0\n");
+ fprintf (fp, "imagedata byte order := littleendian\n");
+ fprintf (fp, "!PET STUDY (General) :=\n");
+ fprintf (fp, "!PET data type := 3D-Sinogram\n");
+ fprintf (fp, "process status := Reconstructed\n");
+ fprintf (fp, "!number format := unsigned short\n");
+ fprintf (fp, "!number of bytes per pixel := 2\n");
+ fprintf (fp, "number of dimensions := 3\n");
+ fprintf (fp, "matrix axis label [1] := x\n");
+ fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);
+ fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(m_RingDiameter/(m_NOfBins-1.0)));
+
+ fprintf (fp, "matrix axis label [2] := y\n");
+ fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);
+
+ fprintf (fp, "scaling factor (degree/pixel) [2] := %f\n",(float)(360./m_NOfAngles));
+
+ fprintf (fp, "matrix axis label [3] := z\n");
+ fprintf (fp, "!matrix size [3] := %i\n",m_NOfPlanes*m_NOfPlanes);
+ fprintf (fp, "scaling factor (mm/pixel) [3] := %f\n",(float)(m_AxialDistance/(m_NOfPlanes-1.0)));
+
+ fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);
+ fprintf (fp, "number of time frames := 1\n");
+ fprintf (fp, "image scaling factor[1] := 1\n");
+ fprintf (fp, "data offset in bytes[1] := 0\n");
+ fprintf (fp, "quantification units := 1\n");
+ fprintf (fp, "!END OF INTERFILE := \n");
+
+ fclose(fp);
+ //(size_t)(m_NOfBins*m_NOfAngles*m_NOfPlanes*m_NOfPlanes);
+
+ }else{
+
+ strcpy(name_hv, m_Filename);
+ strcpy(name_v,m_Filename);
+
+ strcat(name_hv, ".hs"); // STIR extension: .hs .s
+ strcat(name_v, ".s");
+
+ fp=fopen(name_hv, "w");
+
+ fprintf (fp, "!INTERFILE := \n");
+ fprintf (fp, "name of data file := %s\n",name_v);
+ fprintf (fp, "!GENERAL DATA := \n");
+ fprintf (fp, "!GENERAL IMAGE DATA :=\n");
+ fprintf (fp, "!type of data := PET\n");
+ // fprintf (fp, "!version of keys := 3.3\n"); STIR format is not 3.3 (almost but not completely), ERROR in STIR if it is not removed
+ // fprintf (fp, "!data offset in bytes := 0\n");
+ fprintf (fp, "imagedata byte order := littleendian\n");
+ fprintf (fp, "!PET STUDY (General) :=\n");
+ fprintf (fp, "!PET data type := Emission\n");
+ fprintf (fp, "applied corrections := {arc correction}\n"); // {none}\n");
+ // fprintf (fp, "process status := Reconstructed\n");
+ fprintf (fp, "!number format := unsigned integer\n");
+ fprintf (fp, "!number of bytes per pixel := 2\n");
+
+ fprintf (fp, "number of dimensions := 4\n");
+ fprintf (fp, "matrix axis label [4] := segment\n");
+ fprintf (fp, "!matrix size [4] := %i\n",m_MaxRingDifference*2 + 1);
+ // fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(d_FOV/(m_NOfBins-1)));
+ fprintf (fp, "matrix axis label [3] := axial coordinate\n");
+ fprintf (fp, "!matrix size [3] := { ");
+ if (m_MaxRingDifference==0)
+ {
+ fprintf (fp, "%i}\n", m_NOfPlanes);
+ }else{
+ for(int m=m_NOfPlanes-m_MaxRingDifference;m<=m_NOfPlanes;m++) fprintf (fp, "%i,", m);
+ for(int m=m_NOfPlanes-1;m>m_NOfPlanes-m_MaxRingDifference;m--) fprintf (fp, "%i,", m);
+ fprintf (fp, "%i}\n", m_NOfPlanes-m_MaxRingDifference);
+ }
+ fprintf (fp, "matrix axis label [2] := view\n");
+ fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);
+ fprintf (fp, "matrix axis label [1] := tangential coordinate\n");
+ fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);
+
+ fprintf (fp, "minimum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable)
+ fprintf (fp, "%i", -m_MaxRingDifference);
+ for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m);
+ fprintf (fp, "}\n");
+ fprintf (fp, "maximum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable)
+ fprintf (fp, "%i", -m_MaxRingDifference);
+ for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m);
+ fprintf (fp, "}\n");
+
+ fprintf (fp, "inner ring diameter (cm) := %f\n", m_RingDiameter/10); // STIR Required parameter, now assigned to FOV (not detectors)
+ fprintf (fp, "average depth of interaction (cm) := 0.0001\n");
+ fprintf (fp, "default bin size (cm) := %f\n",0.1*((float)m_RingDiameter/((float)m_NOfBins-1.0)) );
+ fprintf (fp, "number of rings := %i\n",m_NOfPlanes );
+ fprintf (fp, "distance between rings (cm) := %f\n", 0.1*((float)m_AxialDistance/(float)(m_NOfPlanes-1)) ); // Axial pixel dimension
+
+ fprintf (fp, "number of detectors per ring := %i\n",m_NOfAngles*2 );
+ // fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);
+ fprintf (fp, "number of time frames := 1\n");
+ fprintf (fp, "image scaling factor[1] := 1\n");
+ fprintf (fp, "data offset in bytes[1] := 0\n");
+ fprintf (fp, "quantification units := 1\n");
+ fprintf (fp, "!END OF INTERFILE := \n");
+
+ fclose(fp);
+
+ }
+ m_Buffer = (SINO_TYPE*) malloc( m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE));
+
+ long unsigned int cont=0;
+ int i,j,k;
+
+ for(k=0;k<m_TotalAxialPlanes;k++){
+ for(j=0;j<m_NOfAngles;j++){
+ for(i=0;i<m_NOfBins;i++){
+ m_Buffer[cont]=m_projections[i][j][k];
+ cont++;
+ }
+ }
+ }
+
+ fp=fopen(name_v, "wb");
+
+ //cout << 4096*sizeof(SINO_TYPE) << endl;
+ int nb=fwrite(m_Buffer,1,m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE), fp);
+ fclose(fp);
+
+#ifndef GAMOS_NO_VERBOSE
+ cout << "PETProjDataMgr::WriteInterfile: File name: " << m_Filename << endl;
+ cout << "PETProjDataMgr::WriteInterfile: Numer of bytes written: " << nb << endl;
+ cout << "PETProjDataMgr::WriteInterfile: Planes = " << m_NOfPlanes << "; bins = " << m_NOfBins << "; ang_views = " << m_NOfAngles << endl;
+ cout << "PETProjDataMgr::WriteInterfile: Dimensions (mm): Transaxial FOV = " << m_RingDiameter << "; Axial FOV = " << m_AxialDistance << " ; Transaxial_pix = " << m_RingDiameter/(m_NOfBins-1) <<"; Plane width = " << m_AxialDistance/(m_NOfPlanes-1) << endl; // Image Axial Pixel(ssrb) == 0.5*(Plane_Width);
+ cout << "... " << endl;
+
+ cout << "PETProjDataMgr::WriteInterfile: Total Coinci: " << m_TotalCoincidences << endl;
+ cout << "PETProjDataMgr::WriteInterfile: Sino3D Coinci: " << m_TotalProjectionCoincidences << endl;
+#endif
+
+}
+
+
+TVector3 PETProjDataMgr::Hits2Digits(const TVector3 &r){
+ if (!m_nch) return r;
+ float smear=0.5;
+
+ if (m_nch<0) smear=m_rnd->Rndm();
+
+ double angle = TMath::ATan2(r.X(),r.Y()); // vrne kot med -pi in pi
+ if (angle<0) angle+=TMath::TwoPi();
+
+ angle= (int(angle/TMath::TwoPi()*TMath::Abs(m_nch))+smear)*TMath::TwoPi()/TMath::Abs(m_nch);
+ //(m_rnd->Rndm()-0.5)*m_AxialDistance;
+ return TVector3(sin(angle), cos(angle),0); // z coordinata ni cisto v redu
+
+}
+
+int PETProjDataMgr::FwdProject(double x,double y, double z, int nmax, TH1*h){
+TVector3 r(x,y,z);
+int h2d=h->InheritsFrom("TH2F");
+double tfac=m_RingDiameter*m_RingDiameter/4-r.Mag2();
+double rfac= m_AxialDistance/m_RingDiameter;
+for (int i=0;i<nmax;i++){
+
+ double phi= m_rnd->Rndm()*TMath::Pi();
+ TVector3 s(1,0,0);
+ s.SetPhi(phi);
+ double sign = (m_rnd->Rndm()>0.5)? 1 : 0;
+ double theta = TMath::ACos(m_rnd->Rndm()*rfac);
+ theta+=sign*TMath::Pi();
+
+ s.SetTheta(theta);
+ double t=r*s;
+ TVector3 rx=r-t*s;
+
+ double d=TMath::Sqrt(t*t+tfac);
+
+ TVector3 r1=rx+d*s;
+ TVector3 r2=rx-d*s;
+
+ //r1=Hits2Digits(r1);
+ //r2=Hits2Digits(r2);
+
+ AddEvent( r1 , r2);
+ if (h!=NULL){
+ TVector3 s1=r2-r1;
+ double s1len= s1.Mag();
+ int niter=int (100*s1len/m_RingDiameter);
+ for (int j=0;j<niter;j++){
+ r2=r1+m_rnd->Rndm()*s1;
+ if (h2d) ((TH2F *) h)->Fill(r2.X(),r2.Y());
+ else ((TH3F *) h)->Fill(r2.X(),r2.Y(),r2.Z());
+ }
+ }
+}
+return 0;
+}
+
+int PETProjDataMgr::FwdProject(TH2F *img, TH2F *h){
+
+for (int i=0;i<img->GetNbinsX();i++) {
+ double x_=img->GetXaxis()->GetBinCenter( i+1 );
+ for (int j=0;j<img->GetNbinsY();j++) {
+ double y_=img->GetYaxis()->GetBinCenter( j+1 );
+ double density= img->GetBinContent(i+1,j+1);
+ if (density>0) FwdProject(x_,y_,m_AxialDistance*(m_rnd->Rndm()-0.5), density,h);
+ }
+}
+return 0;
+}
+
+int PETProjDataMgr::FwdProject(TH3F *img, TH3F *h){
+
+for (int i=0;i<img->GetNbinsX();i++) {
+ double x_=img->GetXaxis()->GetBinCenter( i+1 );
+ for (int j=0;j<img->GetNbinsY();j++) {
+ double y_=img->GetYaxis()->GetBinCenter( j+1 );
+ for (int k=0;k<img->GetNbinsZ();k++) {
+ double z_=img->GetZaxis()->GetBinCenter( k+1 );
+ double density= img->GetBinContent(i+1,j+1,k+1);
+ if (density>0) FwdProject(x_,y_,z_, density,h);
+ }
+ }
+}
+return 0;
+}
+
+TH2F *PETProjDataMgr::Phantom(int kaj){
+TH2F *img= new TH2F("img","Original Image",100,-50,50,100,-50,50);
+
+// izberi sliko 0: kroglice, 1: point source 2: central ball
+switch (kaj){
+
+case 0:
+for (int i=0;i<img->GetNbinsX();i++) {
+ for (int j=0;j<img->GetNbinsY();j++) {
+ double x_=img->GetXaxis()->GetBinCenter( i+1 );
+ double y_=img->GetYaxis()->GetBinCenter( j+1 );
+ double density=1000;
+ if ((x_*x_+y_*y_)<6) img->SetBinContent(i+1,j+1,density);
+
+ density=500; if ((x_-25)*(x_-25)+y_*y_<12) img->SetBinContent(i+1,j+1,density);
+ density=2000; if ((y_-25)*(y_-25)+x_*x_<2) img->SetBinContent(i+1,j+1,density);
+ }
+}
+break;
+
+case 2:
+for (int i=0;i<img->GetNbinsX();i++) {
+ for (int j=0;j<img->GetNbinsY();j++) {
+ double x_=img->GetXaxis()->GetBinCenter( i+1 );
+ double y_=img->GetYaxis()->GetBinCenter( j+1 );
+ double density=1000;
+ if ((x_*x_+y_*y_)<12.5) img->SetBinContent(i+1,j+1,density);
+ }
+}
+break;
+
+case 1:
+img->Fill(25,25,10000);
+break;
+
+}
+
+return img;
+}
Index: praktikum/petrecofmf/pedestals.cxx
===================================================================
--- praktikum/petrecofmf/pedestals.cxx (nonexistent)
+++ praktikum/petrecofmf/pedestals.cxx (revision 113)
@@ -0,0 +1,71 @@
+int pedestals(char *fname=NULL){
+ char fpede[255]="pedestals.dat";
+
+
+ if (fname!=NULL){
+ cout << "Rootname:" << fname << endl;
+ delete gROOT->GetListOfFiles()->FindObject(fname); // clear memory of file name
+ if( gSystem->AccessPathName(fname) ) {
+ cout << endl << "File: " << fname << " not there!!!" << endl << endl;
+ return 0;
+ }
+ TFile *rootfile = new TFile(fname);
+ }
+
+ TH1I *gICh = (TH1I * ) gDirectory->Get("mapping");
+ char hname[50];
+ char hn[50];
+ TCanvas *c;
+
+
+
+ fp=fopen(fpede,"w");
+ for (int i=0;i<4;i++) {
+ sprintf(hname,"cy%d",i);
+ sprintf(hn,"Single Pmt %d",i);
+ c =new TCanvas(hname,hn, 1000,1000);
+ c->Divide(4,4);
+ for (int j=0;j<16;j++){
+ c->cd(gICh->GetBinContent(j+1)+1)->SetLogy(1);
+
+ sprintf(hn,"ach%d",i*16+j);
+ TH1F*h= (TH1F * ) gDirectory->Get(hn);
+ float mean=h->GetMean();
+ h->Fit("gaus","","",mean-200,mean+100);
+
+ TF1 *res = h->GetFunction("gaus");
+
+ fprintf(fp, "%d %g\n",i*16+j,res->GetParameter(1));
+
+ }
+ c->Modified();
+ c->Update();
+ }
+ fclose(fp);
+
+ printf("pedestals are stored in file %s\n", fpede);
+
+
+ /*
+ char fsumpede[255]="sumpedestals.dat";
+ FILE *fp=fopen(fsumpede,"w");
+ c =new TCanvas("c","c", 600,600);
+ c->Divide(2,2);
+ for (int i=0;i<4;i++){
+ sprintf(hn,"pmt%d",i);
+ c->cd(i+1);
+ TH1F *h= ((TH1F * )gDirectory->Get(hn));
+ if (h) {
+ h->Fit("gaus", "","",3000, 6000);
+ TF1 *res = h->GetFunction("gaus");
+ fprintf(fp, "%d %g\n",i,res->GetParameter(1));
+ }
+ }
+ c->Modified();
+ c->Update();
+ fclose(fp);
+ printf("pedestals are stored in file %s\n", fsumpede);
+ */
+
+return 0;
+}
Index: praktikum/petrecofmf/sinoread.cxx
===================================================================
--- praktikum/petrecofmf/sinoread.cxx (nonexistent)
+++ praktikum/petrecofmf/sinoread.cxx (revision 113)
@@ -0,0 +1,152 @@
+#include <stdio.h>
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TStyle.h"
+#include "TColor.h"
+
+void mypalette(){
+
+ Int_t MyPalette[1000];
+ double stops[] = {0.00, 0.34, 0.61, 0.84, 1.00};
+ double red[] = {1.00, 0.84, 0.61, 0.34, 0.00};
+ double green[] = {1.00, 0.84, 0.61, 0.34, 0.00};
+ double blue[] = {1.00, 0.84, 0.61, 0.34, 0.00};
+
+
+ Int_t gs = TColor::CreateGradientColorTable(5, stops, red, green, blue, 999);
+ //for (int i=0;i<999;i++) MyPalette[i] = gs+998-i;
+
+ for (int i=0;i<999;i++) MyPalette[i] = gs+i;
+ gStyle->SetPalette(999, MyPalette);
+}
+
+int imageread(char *fname,int nchx=28, int nchy=56, int nchz=0, int format=0, int size=0, char *opt=""){
+FILE *fp= fopen(fname,"read");
+const int maxbuf=360000;
+float *fbuf= new float[maxbuf];
+unsigned short*ibuf=((unsigned short*) &fbuf[0]);
+if (size==0) {
+ if (format) size=sizeof(float);
+ else size=sizeof(unsigned short);
+}
+int ncount=0;
+
+if(fp){
+
+gStyle->SetPalette(100);
+mypalette();
+TCanvas *c=new TCanvas("c","c",1500,1500);
+int sizx=sqrt(nchz);
+c->Divide(sizx,TMath::Ceil(nchz/double(sizx)) );
+char hname[128];
+
+while (!feof(fp)){
+
+
+int nb=fread(fbuf,size,nchx*nchy,fp);
+printf("nb=%d\n",nb);
+if (nb>0){
+ sprintf(hname,"slice %d",ncount);
+ TH2F *h=new TH2F("h",hname,nchx,-0.5,nchx-0.5,nchy,-0.5,nchy-0.5);
+ float val=0;
+ for (int i=0 ; i<nb;i++){
+ //h->SetBinContent(i%nchx+1, i/nchx+1, fbuf[i]);
+ if (format) val=fbuf[i]; else val=float(ibuf[i]);
+ h->SetBinContent(i/nchy+1,i%nchy+1, val);
+ }
+ c->cd(ncount+1);
+ h->SetMinimum(0);
+ h->Draw(opt);
+}
+
+ncount++;
+}
+fclose(fp);
+c->Modified();
+c->Update();
+}
+
+delete fbuf;
+return 0;
+}
+
+int interfileread(char *fname, char *opt=""){
+ FILE *fp=fopen(fname,"r");
+ if (!fp) {
+ printf("Cannot open %s\n",fname);
+ return 0;
+ }
+ int j=0;
+ int ndim=400;
+ char *line= new char[ndim];
+ int indx[4]={-1,-1,-1,-1};
+ int dim[5];
+ char *cmd,*val,*token;
+ char fimage[100];
+ while (fgets(line,ndim,fp)!=NULL) {
+ char *v= strstr(line,":=");
+ if (v!=NULL) {
+ cmd=line;
+ v[0]=0;
+ val = &v[2];
+ val[strlen(val)-1]=0;
+ printf("%d #%s# %s\n",j++, val, cmd);
+ int axis=-1;
+ if (strstr(cmd,"matrix axis label")!= NULL){
+
+ sscanf(cmd,"matrix axis label [%d]", &axis);
+ if(axis>0 && axis<5){
+ if (strstr(val," x")!= NULL) indx[axis-1]=1;
+ if (strstr(val," y")!= NULL) indx[axis-1]=0;
+ if (strstr(val," z")!= NULL) indx[axis-1]=3;
+
+ if (strstr(val,"view") != NULL) indx[axis-1]=0;
+ if (strstr(val,"tangential")!= NULL) indx[axis-1]=1;
+ if (strstr(val,"segment") != NULL) indx[axis-1]=2;
+ if (strstr(val,"axial") != NULL) indx[axis-1]=3;
+ printf("-->> %s axis=%d indx=%d\n",val, axis,indx[axis-1]);
+
+ }
+ }
+ if (strstr(cmd,"!matrix size")!= NULL){
+ sscanf(cmd,"!matrix size [%d]", &axis);
+ if(axis>0 && axis<5){
+ if (indx[axis-1]<3){ dim[indx[axis-1]]=atoi(val);}
+ else {
+ int sum=0;
+ token = strtok(val,"{,");
+ sum+=atoi(token);
+ while ( (token=strtok(NULL,"{,")) ) sum+=atoi(token);
+ printf("total dim=%d\n",sum);
+ dim[3]=sum;
+ };
+ }
+ printf("axis=%d indx=%d\n",axis,indx[axis-1]);
+ }
+ if (strstr(cmd,"name of data file")!= NULL){
+ sscanf(val,"%s",fimage);
+ printf("image file=%s\n",fimage);
+ }
+ if (strstr(cmd,"!number of bytes per pixel")!= NULL){
+ dim[4]=atoi(val);
+ printf("bytes per pixel=%d\n",dim[4]);
+ }
+ if (strstr(cmd,"!number format")!= NULL){
+ if (strstr(val,"float")!=NULL){ dim[5]=1; };
+ if (strstr(val,"integer")!=NULL){ dim[5]=0; };
+ printf("number format=%d\n",dim[5]);
+ }
+
+ }
+ }
+ for (int i=0;i<6;i++) printf ("dim[%d]=%d\n", i,dim[i]);
+ fclose(fp);
+ imageread(fimage,dim[0], dim[1], dim[3], dim[5], dim[4], opt);
+ return 0;
+}
+
+
+int sinoread(char *fname, char *opt=""){
+ return interfileread(fname,opt);
+}
Index: praktikum/petrecofmf/readdata.C
===================================================================
--- praktikum/petrecofmf/readdata.C (nonexistent)
+++ praktikum/petrecofmf/readdata.C (revision 113)
@@ -0,0 +1,793 @@
+/*
+
+vhodne datoteke:
+
+Mappingi:
+m16.map channel mapping
+modules.map mapping modulov
+
+Kalibracija:
+pedestals.dat ... adc suma
+ph511.dat ... adc fotovrhov vrhov
+
+
+// sumpedestals.dat vsota adc po modulih
+
+./readdata -c config.xml -i input.dat -o output -n nevents -w writentuple -d debuglevel
+
+*/
+#include <stdlib.h>
+#include <stdio.h>
+#include "TNtuple.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TFile.h"
+#include "TMath.h"
+#include "TRandom3.h"
+#include <zlib.h>
+#include <iostream>
+#include <fstream>
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <vector>
+
+#include <libxml/tree.h>
+#include <libxml/parser.h>
+#include <libxml/xpath.h>
+#include <libxml/xpathInternals.h>
+
+#include "PETProjDataMgr.h"
+
+
+// en krog je 4790 korakov
+#define NSTEPS 4790
+// nt->Draw("px1/max1:py1/max1>>hmy(200,-1,4,200,-1,4)","max1>500","colz")
+//#define MAXCH 64
+#define MAXCH 72
+
+
+#define RUNREC_ID 1
+#define EVTREC_ID 2
+#define POSREC_ID 3
+
+typedef struct {
+ unsigned int num_events,num_channels,pedestal,xy;
+ int nx,x0,dx,ny,y0,dy;
+} RUNREC;
+RUNREC *runrec; // start header: appears once at the start of the file
+
+typedef struct {
+ int mikro_pos_x,set_pos_x;
+ int num_iter_y,mikro_pos_y,set_pos_y;
+} POSREC;
+POSREC *posrec; // position header: appears at every change of position
+
+class Channel {
+public:
+ Channel(){};
+ ~Channel(){};
+ int ix;
+ int iy;
+ int idx;
+ void SetIxIy(int a,int b){ix=a;iy=b;idx=ix+4*iy;};
+};
+
+class Module {
+public:
+ Module(){};
+ ~Module(){};
+ double r;
+ double phi;
+ void SetRPhi(double rr, double pphi){r=rr;phi=pphi*TMath::Pi()/180.;};
+};
+
+class Geometry {
+
+public:
+
+
+ ~Geometry(){ m_calibration->Close(); };
+ Channel ch[16];
+ Module module[4];
+
+
+ TRandom3 *m_rnd; // random generator
+ TH1I *m_crystalid[4];// kalibracijski histogrami
+ char * m_modulesmapName; // ime datoteke s podatki o pozicijah modulov
+ char * m_channelsmapName; // ime datoteke s podatki o pozicijah elektronskih kanalov
+ //char * m_sumpedestalsName;
+ char * m_pedestalsName; // ime datoteke s podatki o pedestalih
+ char * m_photopeakName; // ime datoteke s podatki o vrhovih
+ char * m_calibrationrootName; // ime datoteke s kalibracijskimi histogrami
+ TFile *m_calibration; // datoteke s kalibracijskimi histogrami
+
+ double m_dx; // pitch x kristalov
+ double m_dy; // pitch y kristalov
+ int m_nx; // stevilo kristalov v x smeri
+ int m_ny; // stevilo kristalov v y smeri
+
+ float gPedestals[MAXCH];
+ float gPeak[MAXCH];
+
+ //float gSumPedestals[MAXCH];
+ double m_peakscaling; // scaling za adcje
+ double m_threshold; // threshold za upostevanje zadetkov na kanalih
+ int m_debug;
+
+
+ //---------------------------------------------------
+ int SetThreshold(int threshold){
+ m_threshold=threshold;
+ return 0;
+ };
+ double GetThreshold(){
+ return m_threshold;
+ }
+
+
+ int readfile(const char *fname, float *x, int nmax, int defaultvalue=0){
+ int id;
+ float ix;
+ std::ifstream f(fname);
+ if (f.is_open()){
+ std::cout << "Reading ... " << fname << std::endl;
+ while (f.good()){
+ f >> id >> ix;
+ if (id<nmax) x[id]=ix;
+ if (m_debug) std::cout << fname << " " << id << " " << ix <<" "<< std::endl;
+ }
+ f.close();
+ } else {
+ std::cout << "Cannot read " << fname << std::endl;
+ std::cout << "Default value will be used :" << defaultvalue << std::endl;
+ for (int i=0;i<nmax;i++){
+ x[i]=defaultvalue;
+ }
+ }
+ return 0;
+ };
+
+ void ReadChannelMap(const char *fname){
+ int id,ix,iy;
+ char line[400];
+ std::ifstream f(fname);
+ if (f.is_open()){
+ std::cout << "Reading ... " << fname << std::endl;
+ f.getline(line,400);
+ while (f.good()){
+ f >> id >> ix >> iy;
+ if (id<16) ch[id].SetIxIy(ix,iy);
+ if (m_debug) std::cout << fname << " " << id << " " << ix <<" "<< iy << std::endl;
+ }
+ f.close();
+ } else {
+ std::cout << "Cannot read " << fname << std::endl;
+ }
+
+ };
+
+ void ReadModuleMap(const char *fname){
+ int id;
+ double r,phi;
+ char line[400];
+ std::ifstream f(fname);
+ if (f.is_open()){
+ std::cout << "Reading ... " << fname << std::endl;
+ f.getline(line,400);
+ while (f.good()){
+ f >> id >> r >> phi;
+ if (id<4) module[id].SetRPhi(r,phi);
+ if (m_debug) std::cout << fname << " " << id << " " << r <<" "<< phi << std::endl;
+ }
+ f.close();
+ } else {
+ std::cout << "Cannot read " << fname << std::endl;
+ }
+
+ };
+
+ int getvalue(xmlXPathContextPtr xpathCtx, const char* xpathExpr, std::vector<xmlChar *> &retval) {
+
+ xmlXPathObjectPtr xpathObj;
+ xmlNodeSetPtr nodes;
+ xmlChar * ret=NULL;
+ int size;
+ int i;
+ assert(xpathExpr);
+ /* Evaluate xpath expression */
+ xpathObj = xmlXPathEvalExpression(BAD_CAST xpathExpr, xpathCtx);
+ if(xpathObj == NULL) {
+ fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);
+ return retval.size();
+ }
+
+ nodes = xpathObj->nodesetval;
+ size = (nodes) ? nodes->nodeNr : 0;
+ for(i = 0; i< size; i++) {
+ ret = xmlNodeGetContent(nodes->nodeTab[i]);
+ if(ret) {
+ retval.push_back(ret);
+ printf("[%d] %s Return value: %s\n", i, xpathExpr, ret);
+ }
+ }
+
+ /* Cleanup of XPath data */
+ xmlXPathFreeObject(xpathObj);
+
+ return retval.size();
+ }
+
+ char * getvalue(xmlXPathContextPtr xpathCtx, const char* xpathExpr, int which) {
+
+ xmlXPathObjectPtr xpathObj;
+ xmlNodeSetPtr nodes;
+ xmlChar * ret=NULL;
+ int size;
+ int i;
+ assert(xpathExpr);
+
+ xpathObj = xmlXPathEvalExpression(BAD_CAST xpathExpr, xpathCtx);
+ if(xpathObj == NULL) {
+ fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);
+ return NULL;
+ }
+
+ nodes = xpathObj->nodesetval;
+ size = (nodes) ? nodes->nodeNr : 0;
+ for(i = 0; i< size; i++) {
+ ret = xmlNodeGetContent(nodes->nodeTab[i]);
+ if(ret) {
+ printf("%s=%s\n", xpathExpr, ret);
+ if (which == i) break;
+ }
+ }
+
+ xmlXPathFreeObject(xpathObj);
+ return (char *) ret;
+
+ }
+
+ int readxmlconfig(const char *fname){
+ xmlInitParser();
+ LIBXML_TEST_VERSION
+
+ /* Load XML document */
+ xmlDocPtr doc = xmlParseFile(fname);
+ if (doc == NULL) {
+ fprintf(stderr, "Error: unable to parse file \"%s\"\n", fname);
+ return(-1);
+ } else {
+ std::cout << "Reading ... " << fname << std::endl;
+ }
+
+ /* Create xpath evaluation context */
+ xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc);
+ if(xpathCtx == NULL) {
+ fprintf(stderr,"Error: unable to create new XPath context\n");
+ xmlFreeDoc(doc);
+ return(-1);
+ }
+
+ m_pedestalsName = getvalue(xpathCtx, "//pedestals" , 0);
+ //m_sumpedestalsName = getvalue(xpathCtx, "//sumpedestals", 0);
+ m_photopeakName = getvalue(xpathCtx, "//photopeak" , 0);
+ m_modulesmapName = getvalue(xpathCtx, "//modules" , 0);
+ m_channelsmapName = getvalue(xpathCtx, "//channels" , 0);
+ m_calibrationrootName = getvalue(xpathCtx, "//channelcalibration", 0);
+ m_threshold = atoi(getvalue(xpathCtx, "//adcthreshold", 0));
+ m_nx = atoi(getvalue(xpathCtx, "//nofcrystalsx", 0));
+ m_ny = atoi(getvalue(xpathCtx, "//nofcrystalsx", 0));
+ m_dx = atof(getvalue(xpathCtx, "//crystalpitchx", 0));
+ m_dy = atof(getvalue(xpathCtx, "//crystalpitchy", 0));
+
+ xmlXPathFreeContext(xpathCtx);
+ xmlFreeDoc(doc);
+ xmlCleanupParser();
+ xmlMemoryDump();
+ return 0;
+ };
+
+
+ double GetEnergy(int ch, int adc){
+ return (adc-gPedestals[ch])/(gPeak[ch]-gPedestals[ch])*m_peakscaling;
+ }
+ int GetRealPosition(int ipmt, float px,float py,float &cx,float &cy){
+ // iz CoG izracuna pozicijo na kristalu
+
+ int binx = m_crystalid[ipmt]->GetXaxis()->FindBin(px);
+ int biny = m_crystalid[ipmt]->GetYaxis()->FindBin(py);
+ int crystalid = m_crystalid[ipmt]->GetBinContent(binx,biny);
+ cx= ( crystalid % m_nx - m_nx/2. + 0.5 ) * m_dx;
+ cy= ( crystalid / m_nx - m_ny/2. + 0.5 ) * m_dy;
+
+
+ // razmazi pozicijo po povrsini kristala
+ double rndx = (m_rnd->Rndm()-0.5) * m_dx;
+ double rndy = (m_rnd->Rndm()-0.5) * m_dy;
+ cx+= rndx;
+ cy+= rndy;
+ return crystalid;
+ }
+
+ TVector3 GetGlobalPosition(int ipmt, double angle, float cx, float cy){
+
+ double phi = angle + module[ipmt].phi;
+ double &r = module[ipmt].r;
+ double sinphi = sin(phi);
+ double cosphi = cos(phi);
+
+ TVector3 rv( cosphi* r + sinphi * cx , sinphi* r - cosphi * cx, cy);
+ //if (m_debug) printf( "[%d] phi=%f %f global ( %f , %f , %f )\n" ,ipmt, module[ipmt].phi,angle, rv.x(), rv.y(), rv.z());
+ return rv;
+ }
+
+ Geometry(const char *fnameconfig, int debug){
+ m_debug = debug;
+ readxmlconfig(fnameconfig);
+
+ ReadModuleMap(m_modulesmapName);
+ ReadChannelMap(m_channelsmapName);
+
+ std::cout << "Reading ... " << m_calibrationrootName << std::endl;
+ m_calibration = new TFile(m_calibrationrootName);
+ for (int i=0; i<4;i++) {
+ char hn[256];
+ sprintf(hn,"pmt1%d_calib",i);
+ m_crystalid[i]= (TH1I *) m_calibration->Get(hn);
+ m_crystalid[i]->ls();
+ }
+
+ //readfile(m_sumpedestalsName ,gSumPedestals,4);
+ m_peakscaling=3000;
+
+ readfile(m_pedestalsName,gPedestals,4*16, 0);
+ readfile(m_photopeakName,gPeak,4*16, m_peakscaling);
+
+ m_rnd= new TRandom3();
+ };
+
+};
+
+
+class readdata {
+private:
+
+ TNtuple* gNtuple;
+
+
+ TH1F* m_Adc[MAXCH];
+ TH1F* m_AdcCut[MAXCH];
+ TH2F* m_Adc_vs_Nhits[MAXCH];
+ TH2F* m_Adc_vs_Sum[MAXCH];
+ TH2F* m_Adc_vs_Sum_Uncorected[MAXCH];
+ TH2F* m_Nhits;
+ TH1F* m_AdcSum[6];
+ TH2F* m_CenterOfGravity[6];
+ TH2F* m_ReconstructedPosition[6];
+ TH2F* m_CenterOfGravityforChAboveThreshold[6];
+ TH2F* m_GlobalPosition;
+ TH1F* m_AdcSumCluster[6];
+ TH2F* m_MaxAdc[4];
+ TH2F* m_SumAdc[4];
+
+ Geometry *m_geo; // pozicije modulov in kanalov
+ float gSum[6];
+ float gRawSum[6];
+ float gMax[6];
+ float gSumCluster[6];
+ float gNtdata[MAXCH*3];
+ int gNabove[5];
+ int m_neve; // number of events to process
+
+public:
+ PETProjDataMgr *m_projector; // projektor za sinograme
+
+ double m_rotation; // trenutna rotacija setupa
+ double m_threshold; //
+
+ double gData[MAXCH]; // korigirani ADC ji
+ double gAdc[MAXCH]; // raw ADC ji
+ int m_debug; // debug izpisi
+ int m_write; // ce je ta flag nastavljen, se v root file izpisuje ntuple
+ int m_coincidences; // stevilo koincidenc
+
+ //---------------------------------------------------
+ int Initialize(){
+
+ if (m_write) {
+ char varlist[1024], tmplist[1024];
+ sprintf( varlist,"sum0:sum1:sum2:sum3:nh0:nh1:nh2:nh3:c0:c1:c2:c3:px0:px1:px2:px3:py0:py1:py2:py3:rx0:rx1:rx2:rx3:ry0:ry1:ry2:ry3:max0:max1:max2:max3");
+ for (int i=0;i<MAXCH;i++) {
+ sprintf(tmplist,"%s",varlist);
+ sprintf(varlist,"%s:a%d",tmplist,i);
+ }
+ /*
+ energija .. (adc -pedestal) * scaling
+
+ sum0 .. sum3 vsota energij po vseh kanalih fotopomnozevalke
+ nh0 .. nh4 stevilo zadetkov na fotopomnozevalki (adc nad thresholdom)
+ c0 .. c4 vsota energij za kanale nad thresholdom
+ px0 .. px4 x koordinata COG na fotopomnozevalki
+ py0 .. py4 y koordinata COG na fotopomnozevalki
+ max0 .. max3 maximalna energija na kanalu na fotopomnozevalki
+ a0 .. a63 energija na i-tem kanalu
+ */
+ gNtuple = new TNtuple("nt","Pet RAW data",varlist);
+ printf ("#%s#\n",varlist);
+ }
+ m_Nhits = new TH2F("hnabove","Stevilo zadetkov nad nivojem diskriminacije",4,-0.5,3.5, 16,+0.5,16.5 ); // stevilo hitov nad thresholdom
+ m_GlobalPosition = new TH2F("globalxy","Rekonstruirana koordinata zadetka",200,-80,80,200,-80,80); // reconstruirana koordinata v globalnem sistemu
+
+
+ char hname[0xFF];
+ char hn[0xFF];
+ for (int i=0;i<MAXCH;i++){
+ sprintf(hname,"Raw ADC Ch. %d ;ADC;N",i);
+ sprintf(hn,"ach%d",i);
+ m_Adc[i] = new TH1F(hn,hname,4000,-0.5,3999.5); // osnovni adcji
+ sprintf(hname,"cADC za dogodke z manj kot x hitov na pmtju Ch. %d ;ADC;N",i);
+ sprintf(hn,"cutch%d",i);
+ m_AdcCut[i] = new TH1F(hn,hname,4000,-0.5,3999.5); // adcji za zadetke z manj kot 7 hiti na PMTju
+ sprintf(hname,"cADC proti stevilu zadetkov na pmtju Ch. %d ;ADC;N",i);
+ sprintf(hn,"singlech%d",i);
+ m_Adc_vs_Nhits[i] = new TH2F(hn,hname,200,0,4000, 17,-0.5,16.5); // adc proti stevilu zadetkov na PMTju
+ sprintf(hname,"cADC proti vsoti kanalov na pmtju, Ch. %d ;ADC;ADCsum",i);
+ sprintf(hn,"corrch%d",i);
+ m_Adc_vs_Sum[i] = new TH2F(hn,hname,200,0,4000, 200,0,12000 ); // raw adc proti vsoti kanalov na PMTju
+ sprintf(hname,"Raw ADC proti vsoti kanalov na pmtju, Ch. %d ;ADC;ADCsum",i);
+ sprintf(hn,"adcvssumch%d",i);
+ m_Adc_vs_Sum_Uncorected[i] = new TH2F(hn,hname,100,0,3500, 100,5000,12000 ); // adc proti vsoti kanalov na PMTju
+ }
+
+ for (int i=0;i<4;i++) {
+ sprintf(hname,"Vsota cADC na PMTju %d ;ADC;N",i);
+ sprintf(hn,"pmt%d",i);
+ m_AdcSum[i] = new TH1F(hn,hname,200,0,20000); // vsota adcjev na PMTju
+ sprintf(hname,"Vsota cADC na PMTju za dogodke z manj kot x zadetki %d ;ADC;N",i);
+ sprintf(hn,"clusterpmt%d",i);
+ m_AdcSumCluster[i] = new TH1F(hn,hname,200,0,20000); // vsota adcjev na pmtju, ki so nad threshholdom, za dogodke z manj kot 5 zadetki na PMTju
+
+ sprintf(hname,"Center naboja CoG PMT %d ;x;y",i);
+ sprintf(hn,"pmt1%d",i);
+ m_CenterOfGravity[i] = new TH2F(hn,hname,200,-0.25,3.25,200,-0.25,3.25); // center naboja na PMTju
+
+ sprintf(hname,"Center naboja CoG za zadetke nbad thresholdom PMT %d ;x;y",i);
+ sprintf(hn,"pmt2%d",i);
+ m_CenterOfGravityforChAboveThreshold[i] = new TH2F(hn,hname,200,0,3,200,0,3); // center naboja za zadetke nad thresholdom
+
+ sprintf(hname,"Rekonstruirana koordinata PMT %d ;x(mm);y(mm)",i);
+ sprintf(hn,"pmt3%d",i);
+ m_ReconstructedPosition[i] = new TH2F(hn,hname,200,-12,12,200,-12,12); // korigirana koordinata na pmtju
+
+ sprintf(hname,"SumAdc PMT %d ;SumAdc(a.u.);crystal ID",i);
+ sprintf(hn,"sumadc%d",i);
+ m_SumAdc[i] = new TH2F(hn,hname,200,0,12000,81,-0.5,80.5); // SumADC spekter za razlicne kristale
+ }
+
+ // store mapping
+ TH1I *hx= new TH1I("mapping", "ch mapping", 16, -0.5,15.5);
+ for (int j=0;j<16;j++) hx->Fill(j,m_geo->ch[j].idx);
+
+ return 0;
+ };
+
+ int FillHistograms(){
+ for (int i=0;i<MAXCH;i++) { // zanka cez vse elektronske kanale;
+ float x=gData[i];
+ int ipmt= i/16;
+ float s=gSum[ipmt];
+ // int idx= 2*(ipmt)+64;
+ // const float fSlope1=2;
+ // const float fSlope2=3;
+ // if ( ((s-gSumPedestals[ipmt] ) > fSlope1*(x-gPedestals[i])) && ((s-gSumPedestals[ipmt] ) < fSlope2*(x-gPedestals[ipmt])) ) {
+ if (gNabove[ipmt]<7){
+ m_Adc_vs_Sum[i]->Fill(x,s);
+ m_Adc_vs_Sum_Uncorected[i]->Fill(gAdc[i],gRawSum[ipmt]);
+ m_AdcCut[i]->Fill(x);
+ }
+ }
+
+ TVector3 r_coincidence[2];
+ int f_coincidence[2]={0,0};
+
+ for (int ipmt=0;ipmt<4;ipmt++){ // zanka preko pmtjev
+ int j2= (ipmt/2)*2+1-ipmt%2; // sosednja fotopomnozevalka
+
+ float posx[2]={0,0};
+ float posy[2]={0,0};
+ float sum[2]={0,0};
+ for (int ich=0;ich<16;ich++){ // zanka preko elektronskih kanalov na fotopomnozevalki
+ int ch= ich+ipmt*16;
+ if (gMax[ipmt]>m_threshold) {
+ posx[0]+= gData[ch]*m_geo->ch[ich].ix;
+ posy[0]+= gData[ch]*m_geo->ch[ich].iy;
+ sum[0] += gData[ch];
+
+ if (gData[ch]> 0.2*gMax[ipmt]){ // pri racunanju pozicije upostevaj le kanale, ki imajo vrednost vecjo od ratio*maksimalna na tisti fotopomnozevalki
+ posx[1]+= gData[ch]*m_geo->ch[ich].ix;
+ posy[1]+= gData[ch]*m_geo->ch[ich].iy;
+ sum[1] += gData[ch];
+ }
+ }
+ }
+
+ if (m_write) {
+ gNtdata[12+ipmt]=posx[0]/sum[0];
+ gNtdata[16+ipmt]=posy[0]/sum[0];
+ gNtdata[20+ipmt]=(sum[1]>0)?posx[1]/sum[1]:-10;
+ gNtdata[24+ipmt]=(sum[1]>0)?posy[1]/sum[1]:-10;
+ gNtdata[28+ipmt]=sum[0];
+ }
+ if (gMax[ipmt]>gMax[j2]){ // izberi fotopomnozevalko, ki ima vec zadetkov kot soseda TRG= (P1|P2) & (P3|P4)
+
+ if ( sum[0] > 0 ) {
+ float px=posx[0]/sum[0];
+ float py=posy[0]/sum[0];
+
+ m_CenterOfGravity[ipmt]->Fill(px,py);
+
+ float cx=0; //koordinata na pmtju
+ float cy=0;
+ int crystalID = m_geo->GetRealPosition(ipmt,px,py,cx,cy);
+ m_ReconstructedPosition[ipmt]->Fill(cx,cy);
+ m_SumAdc[ipmt]->Fill(gSum[ipmt], crystalID );
+ if (m_debug > 2) printf("ipmt=%d cx=%f cy=%f\n",ipmt,cx,cy);
+ r_coincidence[ipmt/2] = m_geo->GetGlobalPosition(ipmt, m_rotation, cx, cy) ;
+ f_coincidence[ipmt/2] =1;
+ m_GlobalPosition->Fill( r_coincidence[ipmt/2].x(),r_coincidence[ipmt/2].y());
+
+ }
+ if ( sum[1] > 0 ) {
+
+ float px=posx[1]/sum[1];
+ float py=posy[1]/sum[1];
+ m_CenterOfGravityforChAboveThreshold[ipmt]->Fill(px,py);
+ }
+
+ }
+ }
+
+ for (int i=0;i<MAXCH;i++) if (gData[i]>m_threshold) { m_Adc_vs_Nhits[i]->Fill(gData[i],gNabove[i/16]); }
+ for (int i=0;i<MAXCH;i++) if (gData[i]>m_threshold && gNabove[i/16]<5 ) { gSumCluster[i/16]+=gData[i]; }
+ for (int ipmt=0;ipmt<4;ipmt++) {
+ m_AdcSum[ipmt]->Fill(gSum[ipmt]);
+ m_AdcSumCluster[ipmt]->Fill(gSumCluster[ipmt]);
+ }
+
+ for (int i=0;i<4;i++) m_Nhits->Fill(i,gNabove[i]);
+
+ if (m_write) {
+ for (int i=0;i<4;i++) gNtdata[i]=gSum[i];
+ for (int i=0;i<4;i++) gNtdata[4+i]=gNabove[i];
+ for (int i=0;i<4;i++) gNtdata[8+i]=gSumCluster[i];
+ for (int i=0;i<MAXCH;i++) gNtdata[32+i]=gData[i];
+ gNtuple->Fill(gNtdata);
+ }
+ //
+ // coincidences
+ //
+ if (f_coincidence[0]&&f_coincidence[1]){
+ m_coincidences++;
+ if (m_debug > 1) {
+ printf( "Coincidence ( %f , %f , %f ) ( %f , %f , %f )\n" ,
+ r_coincidence[0].x(), r_coincidence[0].y(), r_coincidence[0].z(),
+ r_coincidence[1].x(), r_coincidence[1].y(), r_coincidence[1].z() );
+ }
+ m_projector->AddEvent( r_coincidence[0] , r_coincidence[1]);
+ }
+ return 0;
+ };
+
+
+int Process( const char *fnamein, int nevetoprocess){
+ gzFile fp=gzopen(fnamein,"r");
+ printf("Process....\n");
+ if (!fp) {
+ printf("File %s cannot be opened!\n", fnamein);
+ return 0;
+ } else {
+ printf("Processing file %s\n", fnamein);
+ }
+ const int bsize=20000;
+
+ unsigned int *buf = new unsigned int[bsize];
+ int nr=0;
+ int hdr[4];
+
+ int nread=0;
+ while (!gzeof(fp) ){
+ int nitems = gzread(fp,hdr, sizeof(int)*4 );
+ if (nitems !=4*sizeof(int)) break;
+ int recid = hdr[0];
+ int nb = hdr[1] - 4*sizeof(int);
+ if ( nb > bsize) {
+ nb=bsize;
+ printf ("Error bsize %d nb=%d\n",bsize,nb);
+ }
+ // int nint = nb/sizeof(int);
+ int n=hdr[3];
+ nitems=gzread(fp, buf, nb);
+ nr++;
+ if (nitems!=nb){ printf("Read Error nb=%d nitems=%d\n",nb,nitems); continue; }
+
+ switch (recid){
+ case RUNREC_ID:
+ runrec=(RUNREC *) (&buf[0]);
+ printf("RUNREC RECID = %u\n",hdr[0]);
+ printf("RUNREC Length = %u\n",hdr[1]);
+ printf("RUNREC File version = %u\n",hdr[2]);
+
+ printf("RUNREC Time = %u\n",hdr[3]);
+ printf("Number of events per step = %u\n",runrec->num_events);
+ printf("Number of channels measured = %u\n",runrec->num_channels);
+ printf("Pedestal = %u\n",runrec->pedestal);
+ printf("Scan type = %u :: 0 -> single data scan :: 1 -> XY position scan\n",runrec->xy);
+ printf("Number of steps in X = %d\n",runrec->nx);
+ printf("Start position X = %d\n",runrec->x0);
+ printf("Step size direction X = %d\n",runrec->dx);
+ printf("Number of steps in Y = %d\n",runrec->ny);
+ printf("Start position Y = %d\n",runrec->y0);
+ printf("Step size direction Y = %d\n",runrec->dy);
+
+ break;
+ case POSREC_ID:
+ posrec=(POSREC *) (&buf[0]);
+ printf("POSREC nx=%d , posx=%d setx=%d posy=%d sety=%d angle(deg)= %3.2f \n",n,posrec->mikro_pos_x,posrec->set_pos_x,posrec->mikro_pos_y,posrec->set_pos_y, posrec->set_pos_x*360./NSTEPS);
+ m_rotation = posrec->set_pos_x*2*TMath::Pi()/NSTEPS;
+ break;
+ case EVTREC_ID:
+ break;
+
+ }
+
+ if (recid!=EVTREC_ID) continue;
+ if (nread++==nevetoprocess) break;
+ int idx=1;
+ int neve=buf[0]/2;
+ //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ for (int nev=0;nev<neve;nev++){
+ int len=buf[idx];
+ int sb =buf[idx+1];
+ unsigned int *pbuf=&buf[idx+2];
+ if (sb!=0xffab) {
+ printf("0x%04x!0xffab len=%d\n",sb,len);
+ break;
+ }
+ // postavi na nic
+#define InitArrayWithAValue(arr,n,x) {for (int i=0;i<n;i++) arr[i]=x;}
+ InitArrayWithAValue( gSum , 4 , 0);
+ InitArrayWithAValue( gRawSum , 4 , 0);
+ InitArrayWithAValue( gMax , 4 , 0);
+ InitArrayWithAValue( gSumCluster, 4 , 0);
+ InitArrayWithAValue( gNabove , 4 , 0);
+ InitArrayWithAValue( gData , MAXCH, 0);
+ InitArrayWithAValue( gAdc , MAXCH, 0);
+ //------------------------------------------------------------
+ for (int i0=0;i0<len-2;i0+=2) {
+
+ int data0 = pbuf[i0];
+ int data1 = pbuf[i0+1];
+ int geo = (data1 >> 11) & 0x1f;
+ int ch = (data1&0x1f) | (geo<<5);
+ int dtype = (data1>>9)&0x3;
+ int adc = data0&0xfff;
+
+ switch (dtype) {
+ case 0x0:
+ if (ch<MAXCH) {
+ m_Adc[ch]->Fill(adc);
+
+ int ipmt = ch/16;
+
+ gAdc[ch]=adc;
+ gRawSum[ipmt]+=adc;
+
+ if (ch<64) gData[ch]=m_geo->GetEnergy(ch,adc); else gData[ch]=adc;
+
+ gSum[ipmt]+=gData[ch];
+ if (gData[ch] >gMax[ipmt] ) gMax[ipmt]= gData[ch];
+ if (gData[ch] >m_threshold ) gNabove[ipmt]++;
+ }
+ break;
+ case 0x10:
+ case 0x11:
+ case 0x01:
+ break;
+ }
+
+ };// for (int i0=0;i0<len-2;i0+=2)
+ //------------------------------------------------------------
+
+ idx+=len+1;
+ FillHistograms();
+ } // for (int nev=0;nev<neve;nev++)
+ //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ } // while (!gzeof(fp) )
+ gzclose(fp);
+ delete buf;
+ return nr;
+ }
+
+ readdata(const char *fnameconfig, std::vector<TString> &fnamein, const char *fnameout, int write2root=0, int nevetoread=-1, int debug=0 ){
+ m_debug = debug;
+ m_write = write2root;
+ m_neve = nevetoread;
+
+ m_geo = new Geometry(fnameconfig,m_debug);
+ m_threshold=m_geo->GetThreshold();
+
+ char rootname[0xFF]; sprintf(rootname,"%s.root",fnameout);
+ TFile *rootFile = new TFile(rootname,"RECREATE");
+ Initialize();
+ printf("Konec inicializacije ...\n");
+
+ m_projector=PETProjDataMgr::GetInstance();
+ m_projector->SetDebug(m_debug);
+ m_projector->SetRingDiameter(2*m_geo->module[0].r);
+ m_projector->SetFilename(fnameout);
+
+ m_coincidences=0;
+
+ int nr = 0;
+ for (unsigned int i=0;i< fnamein.size() ;i++) {
+ int neve = m_neve -nr;
+ if (m_neve == -1) nr += Process(fnamein[i].Data(), m_neve);
+ else if (neve>0) nr += Process(fnamein[i].Data(), neve);
+ }
+ printf("End...\nreaddata::Number of events read=%d\nreaddata::Number of coincidences=%d\n",nr,m_coincidences);
+
+ m_projector->WriteInterfile();
+ rootFile->Write();
+ rootFile->Close();
+ }
+
+ ~readdata(){
+ if (m_geo) delete m_geo;
+ };
+
+};
+
+#ifdef MAIN
+
+
+int main(int argc, char **argv){
+
+ std::cout << "Usage:" << argv[0] << " -c config.xml -i input.dat -o output -n nevents -w writentuple -d debuglevel "<< std::endl;
+ //-----------------------------------------------------------------
+ char configfile[256]="ini/config.xml";
+ char outputfile[256]="/tmp/petfmf";
+ int nevetoread =-1;
+ int writentuple =0;
+ int verbose =0;
+ char c;
+
+ std::vector<TString> inputfilelist;
+ while ((c=getopt (argc, argv, "c:i:o:n:w:d:")) != -1){
+ switch(c){
+ case 'c': sprintf(configfile,"%s",optarg); break;
+ case 'i': inputfilelist.push_back(TString(optarg)); break;
+ case 'o': sprintf(outputfile,"%s",optarg); break;
+ case 'n': nevetoread = atoi(optarg) ; break;
+ case 'w': writentuple = atoi(optarg); break;
+ case 'd': verbose = atoi(optarg); break;
+ default: abort();
+ }
+ }
+ std::cout << "Used: " << argv[0] << " -c " << configfile << " -o " << outputfile << " -n " << nevetoread << " -w " << writentuple << " -d " << verbose ;
+ for (unsigned int i=0;i< inputfilelist.size() ;i++) std::cout << " -i " << inputfilelist[i];
+ std::cout << std::endl;
+ //-----------------------------------------------------------------
+ new readdata(configfile, inputfilelist,outputfile, writentuple, nevetoread, verbose);
+
+ std::cout << "Output data files " << outputfile << std::endl;
+ return 0;
+}
+
+#endif
Index: praktikum/petrecofmf/FBP2D.sh
===================================================================
--- praktikum/petrecofmf/FBP2D.sh (nonexistent)
+++ praktikum/petrecofmf/FBP2D.sh (revision 113)
@@ -0,0 +1,58 @@
+#!/bin/bash
+#
+export STIRPATH=~rok/pet/reco/STIR2.1/STIR
+#
+export PATH=$PATH:$STIRPATH/bin
+EXE=FBP2D
+IMAGEVIEW=manip_image
+HVNAME=$2.hv
+
+TEMPFILE=`tempfile`
+rm -f $TEMPFILE
+echo -n "fbp2dparameters :=
+
+input file := $1
+output filename prefix := $2
+
+;input file := sino3D.hs
+;output filename prefix := output
+
+
+; output image parameters
+; zoom defaults to 1
+zoom := 1
+; image size defaults to whole FOV
+; ;;;;;;;; zakomentirano 8.1.2013 xy output image size (in pixels) := 100
+xy output image size (in pixels) := 200
+
+; can be used to call SSRB first
+; default means: call SSRB only if no axial compression is already present
+;num segments to combine with ssrb := -1
+
+; filter parameters, default to pure ramp
+alpha parameter for ramp filter := 0.5
+cut-off for ramp filter (in cycles) := 0.3
+
+; allow less padding. DO NOT USE
+; (unless you're sure that the object occupies only half the FOV)
+;Transaxial extension for FFT:=1
+
+; back projector that could be used (defaults to interpolating backprojector)
+; Back projector type:= some type
+
+; display data during processing for debugging purposes
+Display level := 1
+end :=
+" >> $TEMPFILE
+
+echo $EXE $TEMPFILE
+$EXE $TEMPFILE
+
+
+$IMAGEVIEW $HVNAME
+root src/sinoread.cxx+\(\"$HVNAME\",\"colz\"\)
+
+echo "za ponovni ogled slike"
+echo $IMAGEVIEW $HVNAME
+echo ali
+echo root src/sinoread.cxx+\(\"$HVNAME\",\"colz\"\)
/praktikum/petrecofmf/FBP2D.sh
Property changes:
Added: svn:executable
Index: praktikum/petrecofmf/xpath2.c
===================================================================
--- praktikum/petrecofmf/xpath2.c (nonexistent)
+++ praktikum/petrecofmf/xpath2.c (revision 113)
@@ -0,0 +1,139 @@
+/*
+
+xpath2 ini/config.xml '//anode'
+
+ * section: XPath
+ * synopsis: Load a document, locate subelements with XPath,
+ * usage: xpath2 <xml-file> <xpath-expr>
+ * test: xpath2 test3.xml '//discarded'
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <vector>
+
+#include <libxml/tree.h>
+#include <libxml/parser.h>
+#include <libxml/xpath.h>
+#include <libxml/xpathInternals.h>
+
+#if defined(LIBXML_XPATH_ENABLED) && defined(LIBXML_SAX1_ENABLED) && \
+ defined(LIBXML_OUTPUT_ENABLED)
+
+
+
+static int getvalue(xmlXPathContextPtr xpathCtx, const xmlChar * xpathExpr, std::vector <xmlChar *> &retval);
+static xmlChar * getvalue(xmlXPathContextPtr xpathCtx, const xmlChar * xpathExpr, int which=0);
+
+
+int
+main(int argc, char **argv) {
+ std::vector<xmlChar *> retval;
+ xmlDocPtr doc;
+ xmlXPathContextPtr xpathCtx;
+ int size;
+
+ xmlInitParser();
+ LIBXML_TEST_VERSION
+
+ /* Load XML document */
+ doc = xmlParseFile(argv[1]);
+ if (doc == NULL) {
+ fprintf(stderr, "Error: unable to parse file \"%s\"\n", argv[1]);
+ return(-1);
+ }
+
+ /* Create xpath evaluation context */
+ xpathCtx = xmlXPathNewContext(doc);
+ if(xpathCtx == NULL) {
+ fprintf(stderr,"Error: unable to create new XPath context\n");
+ xmlFreeDoc(doc);
+ return(-1);
+ }
+
+ for (int i=2; i<argc;i++){
+ size = getvalue(xpathCtx, BAD_CAST argv[i] , retval);
+ printf ("-->size = %d value %s\n",size, getvalue(xpathCtx, BAD_CAST argv[i]) );
+ }
+
+ xmlXPathFreeContext(xpathCtx);
+ xmlFreeDoc(doc);
+ xmlCleanupParser();
+ xmlMemoryDump();
+ return 0;
+}
+
+
+/**
+ * getvalue:
+ * @xpathCtx: the input xmlXPathContextPtr
+ * @xpathExpr: the xpath expression for evaluation.
+ *
+ * evaluates XPath expression
+ *
+ * Returns node content
+ */
+int getvalue(xmlXPathContextPtr xpathCtx, const xmlChar* xpathExpr, std::vector<xmlChar *> &retval) {
+
+ xmlXPathObjectPtr xpathObj;
+ xmlNodeSetPtr nodes;
+ xmlChar * ret=NULL;
+ int size;
+ int i;
+ assert(xpathExpr);
+ /* Evaluate xpath expression */
+ xpathObj = xmlXPathEvalExpression(xpathExpr, xpathCtx);
+ if(xpathObj == NULL) {
+ fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);
+ return retval.size();
+ }
+
+ nodes = xpathObj->nodesetval;
+ size = (nodes) ? nodes->nodeNr : 0;
+ for(i = 0; i< size; i++) {
+ ret = xmlNodeGetContent(nodes->nodeTab[i]);
+ if(ret) {
+ retval.push_back(ret);
+ printf("[%d] %s Return value: %s\n", i, xpathExpr, ret);
+ }
+ }
+
+ /* Cleanup of XPath data */
+ xmlXPathFreeObject(xpathObj);
+
+ return retval.size();
+}
+
+xmlChar * getvalue(xmlXPathContextPtr xpathCtx, const xmlChar* xpathExpr, int which) {
+
+ xmlXPathObjectPtr xpathObj;
+ xmlNodeSetPtr nodes;
+ xmlChar * ret=NULL;
+ int size;
+ int i;
+ assert(xpathExpr);
+ /* Evaluate xpath expression */
+ xpathObj = xmlXPathEvalExpression(xpathExpr, xpathCtx);
+ if(xpathObj == NULL) {
+ fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);
+ return NULL;
+ }
+
+ nodes = xpathObj->nodesetval;
+ size = (nodes) ? nodes->nodeNr : 0;
+ for(i = 0; i< size; i++) {
+ ret = xmlNodeGetContent(nodes->nodeTab[i]);
+ if(ret) {
+ printf("%s Return value: %s\n", xpathExpr, ret);
+ return ret;
+ }
+ }
+
+ /* Cleanup of XPath data */
+ xmlXPathFreeObject(xpathObj);
+
+ return NULL;
+}
+
+#endif
Index: praktikum/petrecofmf/fotovrh.cxx
===================================================================
--- praktikum/petrecofmf/fotovrh.cxx (nonexistent)
+++ praktikum/petrecofmf/fotovrh.cxx (revision 113)
@@ -0,0 +1,63 @@
+/*
+
+Skripta za kalibracijo visine fotovrhov
+Vhod: histogrami ADC vs sum
+
+Navodilo:
+1.klikni na GCUT toolbar
+2. poklikaj najvisji fotovrh v vsakem od histogramov . klikajh v vrstnem redu hoistogramov.
+3. zakljuci na zadnjem vrhu tako, da se enkrat kliknes na isto tocko
+4. ponovi vse za vse 4 fotopomnozevalke
+
+Izhod: Datoteka z visinami fotovrhov po kanalih
+
+ .x fotovrh.cxx("fmfpet.root")
+
+*/
+
+int fotovrh(char *fname, float max=100){
+
+char cmd[256];
+TFile *f= new TFile(fname);
+
+FILE *fp= fopen("fotovrh.dat","w");
+if (!fp) return 0;
+for (int ipmt=0;ipmt<4;ipmt++){
+ sprintf(cmd,"c%d",ipmt);
+ TCanvas *c= new TCanvas(cmd,cmd,1500,1000);
+
+ c->ToggleToolBar();
+ c->Divide(4,4);
+
+ for (int ch=0;ch<16;ch++){
+ c->cd(ch+1);
+ sprintf(cmd,"adcvssumch%d",ipmt*16+ch);
+ TH2 *h = (TH2 *) f->Get(cmd);
+ if (!h) {
+ cout << "Histogram not found " << cmd << endl;
+ continue;
+ } else {
+ if (max>0) h->SetMaximum(max);
+ h->Draw("colz");
+ }
+ }
+
+ c->WaitPrimitive("CUTG");
+ TCutG *mycutg = (TCutG*)gROOT->GetListOfSpecials()->FindObject("CUTG");
+ if (mycutg) {
+ sprintf(cmd,"cutg%d",ipmt);
+ mycutg->SetName(cmd);
+ mycutg->Print();
+ int n= mycutg->GetN();
+ double x,y;
+ for (int k=0;k<n-1;k++){
+ mycutg->GetPoint(k, x, y);
+ printf("%d %f\n",k+ipmt*16,x);
+ fprintf(fp,"%d %f\n",k+ipmt*16,x);
+ }
+ }
+
+}
+fclose(fp);
+return 0;
+}
Index: praktikum/petrecofmf/calibration.cxx
===================================================================
--- praktikum/petrecofmf/calibration.cxx (nonexistent)
+++ praktikum/petrecofmf/calibration.cxx (revision 113)
@@ -0,0 +1,77 @@
+/*
+
+Skripta za kalibracijo pozicije kanalov
+Vhod: histogrami utezenih pozicij
+
+Navodilo:
+1.klikni na GCUT toolbar
+2.poklikaj vrhove vseh kristalov po vrstnem redu zgoraj levo->zgoraj desno ->spodaj desno
+3. zakljuci na zadnjem vrhu tako, da se enkrat kliknes na isto tocko
+
+Izhod: Datoteka s kalibracijskimi histogrami, kjer je st. entrijev enako stevilki kanala
+
+ .x calibration.cxx("/tmp/pet.root","calibration.root")
+
+*/
+
+int calibration(char *fname, char* fout){
+
+char cmd[256];
+TFile *f= new TFile(fname);
+TFile *fnew= new TFile(fout,"RECREATE");
+for (int ii=0;ii<4;ii++){
+ sprintf(cmd,"c%d",ii);
+ TCanvas *c= new TCanvas(cmd,cmd,1500,1000);
+ sprintf(cmd,"pmt1%d",ii);
+ TH2 *h = (TH2 *) f->Get(cmd);
+ if (!h) {
+ cout << "Histogram not found " << cmd << endl;
+ continue;
+ } else {
+
+ }
+ c->ToggleToolBar();
+ c->Divide(2,1);
+ c->cd(1);
+ h->Draw("colz");
+ sprintf(cmd,"pmt1%d_calib",ii);
+ TH2F *h1= h->Clone(cmd);
+ h1->Reset();
+
+ c->WaitPrimitive("CUTG");
+ TCutG *mycutg = (TCutG*)gROOT->GetListOfSpecials()->FindObject("CUTG");
+ if (mycutg) {
+ sprintf(cmd,"cutg%d",ii);
+ mycutg->SetName(cmd);
+ mycutg->Print();
+ int n= mycutg->GetN();
+ double x,y;
+ TAxis *xaxis = h1->GetXaxis();
+ TAxis *yaxis = h1->GetYaxis();
+ for (int i=0; i < xaxis->GetNbins();i++){
+ double x0=yaxis->GetBinCenter(i+1);
+ for (int j=0; j < yaxis->GetNbins();j++){
+ double y0=yaxis->GetBinCenter(j+1);
+ double rmin=1e10;
+ int w=-1;
+ for (int k=0;k<n;k++){
+ mycutg->GetPoint(k, x, y);
+ double r2 = (x0-x)*(x0-x)+(y0-y)*(y0-y);
+ if (r2<rmin || k==0 ) {
+ w=k;
+ rmin=r2;
+ }
+
+ }
+ h1->Fill(x0,y0,w);
+ }
+ }
+ c->cd(2);
+ h1->Draw("colz");
+ h->Write();
+ mycutg->Write();
+ }
+}
+fnew->Write();
+return 0;
+}
Index: praktikum/petrecofmf/PETProjDataMgr.h
===================================================================
--- praktikum/petrecofmf/PETProjDataMgr.h (nonexistent)
+++ praktikum/petrecofmf/PETProjDataMgr.h (revision 113)
@@ -0,0 +1,120 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The GAMOS software is copyright of the Copyright Holders of *
+// * the GAMOS Collaboration. It is provided under the terms and *
+// * conditions of the GAMOS Software License, included in the file *
+// * LICENSE and available at http://fismed.ciemat.es/GAMOS/license .*
+// * These include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GAMOS collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the GAMOS Software license. *
+// ********************************************************************
+//
+
+#ifndef PETProjDataMgr_h
+#define PETProjDataMgr_h
+
+/*---------------------------------------------------------------------------
+ClassName: PETProjDataMgr
+Author: M. Canadas, P. Arce
+Changes: 01/11: creation
+
+-------------------------------------------------------------------------
+// Description
+
+ PET output (List-mode X Y Z) into projection data (sinograms).
+
+ Output data: Sinograms for PET image reconstruction. Interfile format, STIR compatible (.hs/.s),
+ STIR Software for Tomographic Image Reconstruction: http://stir.sourceforge.net/main.htm
+
+-----------------------------------------------------------------------*/
+
+#include <iostream>
+#include "TVector3.h"
+class TH2F;
+
+//#include <math.h>
+#include <map>
+#include <vector>
+#include <set>
+#include <string>
+
+class TH1;
+class TH2F;
+class TH3F;
+class TRandom3;
+
+ typedef unsigned short SINO_TYPE; //!!NOTE: Check "Write_sino3D" (.hv file) if SINO_TYPE changes !!
+
+//------------------------------------------------------------------------
+class PETProjDataMgr
+{
+ private:
+ PETProjDataMgr();
+ int m_Debug;
+ public:
+ ~PETProjDataMgr();
+
+ //! Get the only instance
+ static PETProjDataMgr* GetInstance();
+ void AddEvent(const TVector3& pos1, const TVector3& pos2);
+ //void ReadFile();
+ //PETOutput ReadEvent( G4bool& bEof );
+ void SetProjection(int axialplane, TH2F*h);
+ void WriteInterfile();
+ void SetNPlanes(int n){m_NOfPlanes=n;};
+ void SetNBin(int n){m_NOfBins=n;};
+ void SetNAng(int n){m_NOfAngles=n;};
+ void SetDebug(int n){m_Debug=n;};
+ void SetNumberOfChannels(int n){m_nch=n;};
+ void SetFilename( const char *fname){ m_Filename= TString(fname); };
+ void SetRingDiameter( double x){ m_RingDiameter = x; };
+ void SetAxialDistance( double x){ m_AxialDistance = x; };
+ TH2F *Phantom(int kaj);
+ int FwdProject(TH2F *img, TH2F *h=NULL);
+ int FwdProject(TH3F *img, TH3F *h=NULL);
+
+ int FwdProject(double x,double y, double z, int nmax, TH1 *h=NULL);
+ TVector3 Hits2Digits(const TVector3 &r);
+
+private:
+ static PETProjDataMgr* m_Instance;
+
+ double m_AxialDistance;
+ double m_RingDiameter;
+ int m_NOfPlanes;
+ int m_NOfBins;
+ int m_NOfAngles;
+ int m_MaxRingDifference;
+ int m_OutFormat;
+ int m_nch;
+ int m_TotalAxialPlanes;
+ TString m_Filename;
+ TRandom3 * m_rnd;
+ //G4bool bDumpCout;
+
+ SINO_TYPE ***m_projections;
+ SINO_TYPE *m_Buffer;
+
+ unsigned long int m_TotalCoincidences;
+ unsigned long int m_TotalProjectionCoincidences;
+
+
+ FILE *fp;
+};
+
+#endif
+