Subversion Repositories f9daq

Compare Revisions

No changes between revisions

Ignore whitespace Rev 112 → Rev 113

/petrecofmf/modules.map
File deleted
/petrecofmf/pedestals.cxx
File deleted
/petrecofmf/readdata.C
File deleted
/petrecofmf/xpath2.c
File deleted
/petrecofmf/calibration.cxx
File deleted
/petrecofmf/fotovrh.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
-application/octet-stream
\ No newline at end of property
/petrecofmf/config_zero.xml
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
Property changes:
Deleted: svn:mime-type
-application/xml
\ No newline at end of property
/petrecofmf/config.xml
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
Property changes:
Deleted: svn:mime-type
-application/xml
\ No newline at end of property
/petrecofmf/PETProjDataMgr.C
File deleted
/petrecofmf/pedestals_new.dat
File deleted
/petrecofmf/ph511.dat
File deleted
/petrecofmf/sinoread.cxx
File deleted
/petrecofmf/FBP2D.sh
File deleted
Property changes:
Deleted: svn:executable
/petrecofmf/PETProjDataMgr.h
File deleted
/petrecofmf/m16_reverse.map
File deleted
/petrecofmf/Makefile
File deleted
/petrecofmf/sumpedestals_zero.dat
File deleted
/petrecofmf/plot.cxx
File deleted
/petrecofmf/m16.map
File deleted
/petrecofmf/sumpedestals.dat
File deleted
/praktikum/petrecofmf/Makefile
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
/praktikum/petrecofmf/README
0,0 → 1,2
http://www-f9.ijs.si/wiki/Main/PetRekonstrukcija
 
/praktikum/petrecofmf/sumpedestals_zero.dat
0,0 → 1,4
0 6000
1 5264.29
2 5394.42
3 5261.11
/praktikum/petrecofmf/m16.map
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
/praktikum/petrecofmf/sumpedestals.dat
0,0 → 1,4
0 5413.22
1 5264.29
2 5394.42
3 5261.11
/praktikum/petrecofmf/modules.map
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
/praktikum/petrecofmf/pedestals_zero.dat
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
/praktikum/petrecofmf/pedestals.dat
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
/praktikum/petrecofmf/fotovrh.dat
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
/praktikum/petrecofmf/calibration_all.root
Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream
Property changes:
Added: svn:mime-type
+application/octet-stream
\ No newline at end of property
/praktikum/petrecofmf/config_zero.xml
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
Property changes:
Added: svn:mime-type
+application/xml
\ No newline at end of property
/praktikum/petrecofmf/config.xml
Cannot display: file marked as a binary type.
svn:mime-type = application/xml
Property changes:
Added: svn:mime-type
+application/xml
\ No newline at end of property
/praktikum/petrecofmf/pedestals_new.dat
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
/praktikum/petrecofmf/ph511.dat
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
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
/praktikum/petrecofmf/m16_reverse.map
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
/praktikum/petrecofmf/plot.cxx
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;
}
//---------------------------------------------------
/praktikum/petrecofmf/PETProjDataMgr.C
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;
}
/praktikum/petrecofmf/pedestals.cxx
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;
}
/praktikum/petrecofmf/sinoread.cxx
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);
}
/praktikum/petrecofmf/readdata.C
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
/praktikum/petrecofmf/FBP2D.sh
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\"\)
Property changes:
Added: svn:executable
/praktikum/petrecofmf/xpath2.c
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
/praktikum/petrecofmf/fotovrh.cxx
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;
}
/praktikum/petrecofmf/calibration.cxx
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;
}
/praktikum/petrecofmf/PETProjDataMgr.h
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