diff --git a/slsDetectorCalibration/MovingStat.h b/slsDetectorCalibration/MovingStat.h index 7acba5706..801c5a64a 100755 --- a/slsDetectorCalibration/MovingStat.h +++ b/slsDetectorCalibration/MovingStat.h @@ -10,12 +10,12 @@ void SetN(int i) {n=i;}; int GetN() {return n;}; - void Calc(double x) { + inline void Calc(double x) { if (m_n 0) ? m_newM/m_n : 0.0; } @@ -71,12 +71,12 @@ return ( (m_n > 1) ? m_newM2/m_n : 0.0 ); } - double Variance() const + inline double Variance() const { return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 ); } - double StandardDeviation() const + inline double StandardDeviation() const { return ( (Variance() > 0) ? sqrt( Variance() ) : -1 ); } diff --git a/slsDetectorCalibration/moenchMakeTree.C b/slsDetectorCalibration/moenchMakeTree.C index 447050e1b..049bd18d9 100644 --- a/slsDetectorCalibration/moenchMakeTree.C +++ b/slsDetectorCalibration/moenchMakeTree.C @@ -17,10 +17,44 @@ #include "RunningStat.h" #include "MovingStat.h" #include "moench02ModuleData.h" +#include 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); + meC = me; sigC = sig; + tall->Fill(); + TThread::UnLock(); +} + void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign=1) { moench02ModuleData *decoder=new moench02ModuleData(); @@ -42,13 +76,8 @@ void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign int ix, iy, ir, ic; Double_t data[3][3]; - TTree* tall=new TTree(tname,tname); - tall->Branch("iFrame",&nf,"iframe/I"); - tall->Branch("x",&ix,"x/I"); - tall->Branch("y",&iy,"y/I"); - tall->Branch("data",data,"data[3][3]/D"); - tall->Branch("pedestal",&me,"pedestal/D"); - tall->Branch("rms",&sig,"rms/D"); + nPhotonsStat.Clear(); + nPhotonsStat.SetN(1000); for (ir=0; ir<160; ir++) { for (ic=0; ic<160; ic++) { @@ -59,7 +88,7 @@ void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign for (int irun=runmin; irunreadNextFrame(filebin))) { @@ -70,22 +99,22 @@ void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign dum=0; //no hit tot=0; - if (nf>1000) { + if (nf>nbg) { me=stat[iy][ix].Mean(); sig=stat[iy][ix].StandardDeviation(); val=sign*decoder->getChannelShort(buff,ix,iy)-me; dum=0; //no hit tot=0; - maxNei; + maxNei = 0; 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 = decoder->getChannelShort(buff,ix+ic,iy+ir)-stat[iy+ir][ix+ic].Mean(); - if (decoder->getChannelShort(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! + valNei = sign*decoder->getChannelShort(buff,ix+ic,iy+ir)-stat[iy+ir][ix+ic].Mean(); + if (sign*decoder->getChannelShort(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][ic]; + data[ir+1][ic+1] = valNei; maxNei = max(maxNei,valNei); } } @@ -94,24 +123,28 @@ void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign if (val<(me-3.*sig)) dum=2; //discard negative events! if(maxNei == val && dum == 1){ // this is an event and we are in the center - tall->Fill(); + storeEvent(nf,ix,iy,data,me,sig); nPhotons++; + //cout << "X: " << ix << " Y: " << iy << " val: " << val << " tot: " << tot << endl; + }else{ + dum = 3; //discard events (for pedestal) where sum of the neighbours is too large. } - if (tot>3.*sig) dum=3; //discard events (for pedestal) where sum of the neighbours is too large. + //if (tot>3.*sig && dum == 1){ dum=3;} + //cout << dum; } - if (nf<1000 || dum==0) + if (nfgetChannelShort(buff,ix,iy)); } } delete [] buff; - //cout << "=" ; - nPhotonsStat.Calc(nPhotons); + //cout << "="; cout.flush(); + if(nf>nbg) nPhotonsStat.Calc((double)nPhotons); nf++; } - cout << endl; - cout << " done. Avg. Photons/Frame: " << nPhotonsStat.Mean() << " sig: " << nPhotonsStat.StandardDeviation() << endl; + //cout << endl; + cout << "processed File " << fname << " done. Avg. Photons/Frame: " << nPhotonsStat.Mean() << " sig: " << nPhotonsStat.StandardDeviation() << endl; if (filebin.is_open()) filebin.close(); else @@ -120,15 +153,23 @@ void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign delete decoder; cout << "Read " << nf << " frames" << endl; - tall->Write(); - tall->GetCurrentFile()->Close(); + } +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); + if(argc != 6){ cout << "Usage: inFile outdir tname fileNrStart fileNrEnd" << endl; exit(-1); } + + int nThreads = 15; + TThread *threads[nThreads]; char *inFile = argv[1]; char *outDir = argv[2]; @@ -138,10 +179,37 @@ int main(int argc, char **argv){ TFile *f; char outfname[1000]; + char threadName[1000]; sprintf(outfname,"%s/%s.root",outDir,tName); f=new TFile(outfname,"RECREATE"); - moenchMakeTree(inFile,tName,start,end); + 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-1)*i; + t->runmax = start + (end-start)/(nThreads-1)*(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(); + } + + + //TThread::Ps(); + + for(int i = 0; i < nThreads; i++){ + threads[i]->Join(); + } + + + tall->Write(); + tall->GetCurrentFile()->Close(); } diff --git a/slsDetectorCalibration/slsDetectorData.h b/slsDetectorCalibration/slsDetectorData.h index 519014b7d..1c3d62a42 100644 --- a/slsDetectorCalibration/slsDetectorData.h +++ b/slsDetectorCalibration/slsDetectorData.h @@ -92,7 +92,7 @@ class slsDetectorData { */ - int getChannel(char *data, int ix, int iy=1) { + inline int getChannel(char *data, int ix, int iy=1) { return (*(data+dataMap[iy][ix]))^dataMask[iy][ix]; }; @@ -107,7 +107,7 @@ class slsDetectorData { */ - u_int16_t getChannelShort(char *data, int ix, int iy=1) { + inline u_int16_t getChannelShort(char *data, int ix, int iy=1) { u_int16_t m=dataMask[iy][ix], d=*(u_int16_t*)(data+dataMap[iy][ix]); // cout << ix << " " << iy << " " << dataMap[ix][iy] << endl; // return (*(dd+dataMap[ix][iy]))^((u_int16_t)dataMask[ix][iy]);