/petrecofmf/FBP2D.sh |
---|
File deleted |
Property changes: |
Deleted: svn:executable |
Index: petrecofmf/PETProjDataMgr.h |
=================================================================== |
--- petrecofmf/PETProjDataMgr.h (revision 112) |
+++ petrecofmf/PETProjDataMgr.h (nonexistent) |
@@ -1,120 +0,0 @@ |
-// |
-// ******************************************************************** |
-// * 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 |
- |
Index: petrecofmf/m16_reverse.map |
=================================================================== |
--- petrecofmf/m16_reverse.map (revision 112) |
+++ petrecofmf/m16_reverse.map (nonexistent) |
@@ -1,17 +0,0 @@ |
-#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: petrecofmf/Makefile |
=================================================================== |
--- petrecofmf/Makefile (revision 112) |
+++ petrecofmf/Makefile (nonexistent) |
@@ -1,60 +0,0 @@ |
- |
-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: petrecofmf/plot.cxx |
=================================================================== |
--- petrecofmf/plot.cxx (revision 112) |
+++ petrecofmf/plot.cxx (nonexistent) |
@@ -1,186 +0,0 @@ |
-//--------------------------------------------------- |
-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: petrecofmf/sumpedestals_zero.dat |
=================================================================== |
--- petrecofmf/sumpedestals_zero.dat (revision 112) |
+++ petrecofmf/sumpedestals_zero.dat (nonexistent) |
@@ -1,4 +0,0 @@ |
-0 6000 |
-1 5264.29 |
-2 5394.42 |
-3 5261.11 |
Index: petrecofmf/m16.map |
=================================================================== |
--- petrecofmf/m16.map (revision 112) |
+++ petrecofmf/m16.map (nonexistent) |
@@ -1,17 +0,0 @@ |
-#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: petrecofmf/sumpedestals.dat |
=================================================================== |
--- petrecofmf/sumpedestals.dat (revision 112) |
+++ petrecofmf/sumpedestals.dat (nonexistent) |
@@ -1,4 +0,0 @@ |
-0 5413.22 |
-1 5264.29 |
-2 5394.42 |
-3 5261.11 |
Index: petrecofmf/modules.map |
=================================================================== |
--- petrecofmf/modules.map (revision 112) |
+++ petrecofmf/modules.map (nonexistent) |
@@ -1,5 +0,0 @@ |
-#module ID (connector) r(mm) angle(deg) |
-0 61 0 |
-1 61 22.5 |
-2 61 180 |
-3 61 202.5 |
Index: petrecofmf/pedestals.cxx |
=================================================================== |
--- petrecofmf/pedestals.cxx (revision 112) |
+++ petrecofmf/pedestals.cxx (nonexistent) |
@@ -1,71 +0,0 @@ |
-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: petrecofmf/readdata.C |
=================================================================== |
--- petrecofmf/readdata.C (revision 112) |
+++ petrecofmf/readdata.C (nonexistent) |
@@ -1,793 +0,0 @@ |
-/* |
- |
-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: petrecofmf/xpath2.c |
=================================================================== |
--- petrecofmf/xpath2.c (revision 112) |
+++ petrecofmf/xpath2.c (nonexistent) |
@@ -1,139 +0,0 @@ |
-/* |
- |
-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: petrecofmf/fotovrh.cxx |
=================================================================== |
--- petrecofmf/fotovrh.cxx (revision 112) |
+++ petrecofmf/fotovrh.cxx (nonexistent) |
@@ -1,63 +0,0 @@ |
-/* |
- |
-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: petrecofmf/calibration.cxx |
=================================================================== |
--- petrecofmf/calibration.cxx (revision 112) |
+++ petrecofmf/calibration.cxx (nonexistent) |
@@ -1,77 +0,0 @@ |
-/* |
- |
-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: petrecofmf/pedestals_zero.dat |
=================================================================== |
--- petrecofmf/pedestals_zero.dat (revision 112) |
+++ petrecofmf/pedestals_zero.dat (nonexistent) |
@@ -1,64 +0,0 @@ |
-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: petrecofmf/pedestals.dat |
=================================================================== |
--- petrecofmf/pedestals.dat (revision 112) |
+++ petrecofmf/pedestals.dat (nonexistent) |
@@ -1,64 +0,0 @@ |
-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: petrecofmf/README |
=================================================================== |
--- petrecofmf/README (revision 112) |
+++ petrecofmf/README (nonexistent) |
@@ -1,2 +0,0 @@ |
-http://www-f9.ijs.si/wiki/Main/PetRekonstrukcija |
- |
Index: petrecofmf/fotovrh.dat |
=================================================================== |
--- petrecofmf/fotovrh.dat (revision 112) |
+++ petrecofmf/fotovrh.dat (nonexistent) |
@@ -1,64 +0,0 @@ |
-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: petrecofmf/calibration_all.root |
=================================================================== |
Cannot display: file marked as a binary type. |
svn:mime-type = application/octet-stream |
/petrecofmf/calibration_all.root |
---|
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: 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 |
+ |