mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2026-01-31 03:34:57 +01:00
new separate receiver
This commit is contained in:
Maliakal Dhanya
244
slsReceiverSoftware/slsDetectorCalibration/moenchMakeTree.C
Normal file
244
slsReceiverSoftware/slsDetectorCalibration/moenchMakeTree.C
Normal file
@@ -0,0 +1,244 @@
|
||||
#include <TH1D.h>
|
||||
#include <TH2D.h>
|
||||
#include <TPad.h>
|
||||
#include <TDirectory.h>
|
||||
#include <TEntryList.h>
|
||||
#include <TFile.h>
|
||||
#include <TMath.h>
|
||||
#include <TTree.h>
|
||||
#include <TChain.h>
|
||||
#include <THStack.h>
|
||||
#include <TCanvas.h>
|
||||
#include <stdio.h>
|
||||
#include <deque>
|
||||
#include <list>
|
||||
#include <queue>
|
||||
#include <fstream>
|
||||
#include "RunningStat.h"
|
||||
#include "MovingStat.h"
|
||||
#include "moench02ModuleData.h"
|
||||
#include <TThread.h>
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
||||
//tree variables
|
||||
int xC,yC,iFrameC;
|
||||
double meC,sigC;
|
||||
Double_t dataC[3][3];
|
||||
TTree* tall;
|
||||
|
||||
typedef struct task_s{
|
||||
char *fformat;
|
||||
char *tname;
|
||||
int runmin;
|
||||
int runmax;
|
||||
int sign;
|
||||
} Task;
|
||||
|
||||
void setUpTree(char *tname){
|
||||
tall=new TTree(tname,tname);
|
||||
tall->Branch("iFrame",&iFrameC,"iframe/I");
|
||||
tall->Branch("x",&xC,"x/I");
|
||||
tall->Branch("y",&yC,"y/I");
|
||||
tall->Branch("data",dataC,"data[3][3]/D");
|
||||
tall->Branch("pedestal",&meC,"pedestal/D");
|
||||
tall->Branch("rms",&sigC,"rms/D");
|
||||
}
|
||||
|
||||
inline void storeEvent(int iF,int x,int y, Double_t data[][3], double me, double sig){
|
||||
TThread::Lock();
|
||||
xC = x; yC = y; iFrameC = iF;
|
||||
memcpy(dataC,data,sizeof(Double_t)*3*3);
|
||||
//cout << "X: " << x << " Y: " << y << endl;
|
||||
/* for(int i = 0; i < 3; i++){
|
||||
for(int j = 0; j < 3; j++){
|
||||
cout << "i: " << i << " j: " << j << " dataC " << dataC[i][j] << " data " << data[i][j] << endl;
|
||||
}
|
||||
}*/
|
||||
meC = me; sigC = sig;
|
||||
tall->Fill();
|
||||
TThread::UnLock();
|
||||
}
|
||||
|
||||
void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign=1) {
|
||||
double nThSigma = 3.;
|
||||
moench02ModuleData *decoder=new moench02ModuleData();
|
||||
char *buff;
|
||||
char fname[10000];
|
||||
|
||||
int nf=0;
|
||||
|
||||
int dum, nPhotons;
|
||||
double me, sig, tot, maxNei, val, valNei;
|
||||
|
||||
MovingStat stat[160][160];
|
||||
MovingStat nPhotonsStat;
|
||||
|
||||
ifstream filebin;
|
||||
|
||||
int nbg=20;
|
||||
|
||||
int ix, iy, ir, ic, mx,my;
|
||||
Double_t data[3][3];
|
||||
|
||||
nPhotonsStat.Clear();
|
||||
nPhotonsStat.SetN(1000);
|
||||
|
||||
for (ir=0; ir<160; ir++) {
|
||||
for (ic=0; ic<160; ic++) {
|
||||
stat[ir][ic].Clear();
|
||||
stat[ir][ic].SetN(nbg);
|
||||
}
|
||||
}
|
||||
|
||||
for (int irun=runmin; irun<runmax; irun++) {
|
||||
sprintf(fname,fformat,irun);
|
||||
//cout << "process file " << fname; // cout.flush();
|
||||
filebin.open((const char *)(fname), ios::in | ios::binary);
|
||||
|
||||
while ((buff=decoder->readNextFrame(filebin))) {
|
||||
nPhotons=0;
|
||||
//for (ix=0; ix<160; ix++){
|
||||
for (ix=0; ix<40; ix++){
|
||||
for (iy=0; iy<160; iy++) {
|
||||
|
||||
|
||||
|
||||
dum=0; //no hit
|
||||
// tot=0;
|
||||
|
||||
me=stat[iy][ix].Mean();
|
||||
sig=stat[iy][ix].StandardDeviation();
|
||||
val=sign*(decoder->getChannelShort(buff,ix,iy)-me);
|
||||
|
||||
|
||||
if (nf>nbg) {
|
||||
me=stat[iy][ix].Mean();
|
||||
sig=stat[iy][ix].StandardDeviation();
|
||||
val=sign*decoder->getChannel(buff,ix,iy)-me;
|
||||
|
||||
// dum=0; //no hit
|
||||
tot=0;
|
||||
maxNei = 0;
|
||||
|
||||
if (val>nThSigma*sig)//{ //hit check if neighbors are higher
|
||||
dum=1;
|
||||
for (ir=-1; ir<2; ir++){
|
||||
for (ic=-1; ic<2; ic++){
|
||||
if ((ix+ic)>=0 && (ix+ic)<160 && (iy+ir)>=0 && (iy+ir)<160) {
|
||||
valNei = sign*decoder->getChannel(buff,ix+ic,iy+ir)-stat[iy+ir][ix+ic].Mean();
|
||||
if (sign*decoder->getChannel(buff,ix+ic,iy+ir)>(stat[iy+ir][ix+ic].Mean()+3.*stat[iy+ir][ix+ic].StandardDeviation())) dum=1; //is a hit or neighbour is a hit!
|
||||
tot+=valNei;
|
||||
data[ir+1][ic+1] = valNei;
|
||||
if(valNei/stat[iy+ir][ix+ic].StandardDeviation() > maxNei){
|
||||
maxNei = valNei/stat[iy+ir][ix+ic].StandardDeviation();
|
||||
mx = ir;
|
||||
my = ic;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// }
|
||||
|
||||
if (val<(-nThSigma*sig)) dum=2; //discard negative events!
|
||||
|
||||
if(dum == 1 && mx == 0 && my == 0){ // this is an event and we are in the center
|
||||
storeEvent(nf,ix,iy,data,me,sig);
|
||||
nPhotons++;
|
||||
//cout << "BF X: " << ix << " Y: " << iy << " val: " << val << " tot: " << tot << " me: " << me << " sig: " << sig << endl;
|
||||
}
|
||||
|
||||
|
||||
//if (tot>9.*sig && dum == 1){ dum=3;}
|
||||
//cout << dum;
|
||||
}
|
||||
//if (ix==20 && iy==40)
|
||||
//cout << decoder->getChannelShort(buff,ix,iy)<< " " << val << " " << me << " " << sig << " " << dum << endl;
|
||||
|
||||
if ( dum==0) {
|
||||
|
||||
stat[iy][ix].Calc(decoder->getChannelShort(buff,ix,iy));
|
||||
}
|
||||
if (nf<nbg || dum==0)
|
||||
stat[iy][ix].Calc(sign*decoder->getChannel(buff,ix,iy));
|
||||
|
||||
}
|
||||
}
|
||||
delete [] buff;
|
||||
//cout << "="; cout.flush();
|
||||
if(nf>nbg) nPhotonsStat.Calc((double)nPhotons);
|
||||
nf++;
|
||||
}
|
||||
//cout << endl;
|
||||
cout << "processed File " << fname << " done. Avg. Photons/Frame: " << nPhotonsStat.Mean() << " sig: " << nPhotonsStat.StandardDeviation() << " " << runmax-irun << " files need processing" << endl;
|
||||
if (filebin.is_open())
|
||||
filebin.close();
|
||||
else
|
||||
cout << "could not open file " << fname << endl;
|
||||
}
|
||||
|
||||
delete decoder;
|
||||
cout << "Read " << nf << " frames" << endl;
|
||||
|
||||
|
||||
}
|
||||
|
||||
void *moenchMakeTreeTask(void *p){
|
||||
Task *t = (Task *)p;
|
||||
moenchMakeTree(t->fformat,t->tname,t->runmin,t->runmax,t->sign);
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
//to compile: g++ -DMYROOT -g `root-config --cflags --glibs` -o moenchMakeTree moenchMakeTree.C
|
||||
int main(int argc, char **argv){
|
||||
if(argc != 6){ cout << "Usage: inFile outdir tname fileNrStart fileNrEnd" << endl; exit(-1); }
|
||||
|
||||
int nThreads = 10;
|
||||
TThread *threads[nThreads];
|
||||
|
||||
char *inFile = argv[1];
|
||||
char *outDir = argv[2];
|
||||
char *tName = argv[3];
|
||||
int start = atoi(argv[4]);
|
||||
int end = atoi(argv[5]);
|
||||
|
||||
TFile *f;
|
||||
char outfname[1000];
|
||||
char threadName[1000];
|
||||
|
||||
sprintf(outfname,"%s/%s.root",outDir,tName);
|
||||
f=new TFile(outfname,"RECREATE");
|
||||
|
||||
cout << "outputfile: " << outfname << endl;
|
||||
setUpTree(tName);
|
||||
|
||||
for(int i = 0; i < nThreads; i++){
|
||||
sprintf(threadName,"t%i",i);
|
||||
Task *t = (Task *)malloc(sizeof(Task));
|
||||
t->fformat = inFile;
|
||||
t->tname = tName;
|
||||
t->sign = 1.;
|
||||
t->runmin = start + (end-start)/(nThreads)*i;
|
||||
t->runmax = start + (end-start)/(nThreads)*(i+1);
|
||||
if(i == nThreads - 1) t->runmax = end;
|
||||
cout << "start thread " << i << " start: " << t->runmin << " end " << t->runmax << endl;
|
||||
threads[i] = new TThread(threadName, moenchMakeTreeTask, t);
|
||||
threads[i]->Run();
|
||||
//moenchMakeTreeTask(t);
|
||||
}
|
||||
|
||||
|
||||
//TThread::Ps();
|
||||
|
||||
for(int i = 0; i < nThreads; i++){
|
||||
threads[i]->Join();
|
||||
}
|
||||
|
||||
|
||||
tall->Write();
|
||||
tall->Print();
|
||||
f->Close();
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user