moench multithread tree maker

git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@5 113b152e-814d-439b-b186-022a431db7b5

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<n) Add(x);
 	  else Push(x);
 	}
 
-	void Add(double x) {
+	inline void Add(double x) {
 	  m_n++;
 
             // See Knuth TAOCP vol 2, 3rd edition, page 232
@@ -38,7 +38,7 @@
 	}
 
 
-        void Push(double x)
+        inline void Push(double x)
         {
 	  if (m_n == 1)
             {
@@ -62,7 +62,7 @@
             return m_n;
         }
 
-        double Mean() const
+        inline double Mean() const
         {
             return (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 );
         }

slsDetectorCalibration/moenchMakeTree.C

@@ -17,10 +17,44 @@
 #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);
+  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);
+  nPhotonsStat.Clear();
-  tall->Branch("iFrame",&nf,"iframe/I");
+  nPhotonsStat.SetN(1000);
-  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");
 
   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; irun<runmax; irun++) {
     sprintf(fname,fformat,irun);
-    cout << "process file " << fname;
+    //cout << "process file " << fname; // cout.flush();
     filebin.open((const char *)(fname), ios::in | ios::binary);
    
     while ((buff=decoder->readNextFrame(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();
+		    valNei = sign*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!
+		    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 (nf<nbg || dum==0) 
 	      stat[iy][ix].Calc(sign*decoder->getChannelShort(buff,ix,iy));
 	   
 	  }
       }
       delete [] buff;
-      //cout << "=" ;
+      //cout << "="; cout.flush();
-      nPhotonsStat.Calc(nPhotons);
+      if(nf>nbg) nPhotonsStat.Calc((double)nPhotons);
       nf++;
     }
-    cout << endl;
+    //cout << endl;
-    cout << " done. Avg. Photons/Frame: " << nPhotonsStat.Mean() << " sig: " <<  nPhotonsStat.StandardDeviation() << 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();
 }
 

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]);