Added Moench03 common mode and photon finder/pedestal funcs

slsDetectorCalibration/moench03CommonMode.h

@@ -0,0 +1,45 @@
+#ifndef MOENCH03COMMONMODE_H
+#define MOENCH03COMMONMODE_H
+
+#include "commonModeSubtraction.h"
+
+class moench03CommonMode : public commonModeSubtraction {
+  /** @short class to calculate the common mode noise for moench02 i.e. on 4 supercolumns separately */
+ public:
+  /** constructor -  initalizes a commonModeSubtraction with 4 different regions of interest 
+      \param nn number of samples for the moving average
+  */
+  moench03CommonMode(int nn=1000) : commonModeSubtraction(nn,32){} ;
+    
+
+    /** add value to common mode as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
+	\param val value to add to the common mode
+	\param ix pixel coordinate in the x direction
+	\param iy pixel coordinate in the y direction
+    */
+    virtual void addToCommonMode(double val, int ix=0, int iy=0) {
+      // (void) iy;  
+      int isc=ix/25+(iy/200)*16;
+      if (isc>=0 && isc<nROI) {
+	cmPed[isc]+=val; 
+	nCm[isc]++;
+      }
+    };
+     /**returns common mode value as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
+	\param ix pixel coordinate in the x direction
+	\param iy pixel coordinate in the y direction
+	\returns common mode value
+    */
+    virtual double getCommonMode(int ix=0, int iy=0) {
+      (void) iy;
+      int isc=ix/25+(iy/200)*16;
+      if (isc>=0 && isc<nROI) { 
+	if (nCm[isc]>0) return cmPed[isc]/nCm[isc]-cmStat[isc].Mean(); 
+      }
+      return 0;
+    };
+
+};
+
+
+#endif

slsDetectorCalibration/moench03CtbData.h

@@ -29,7 +29,10 @@ class moench03CtbData : public slsDetectorData<uint16_t> {
   moench03CtbData(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*2*32, NULL, NULL) , nadc(32), sc_width(25), sc_height(200) {
 
     
-    int adc_nr[32]={200,225,250,275,300,325,350,375,0,25,50,75,100,125,150,175,175,150,125,100,75,50,25,0,375,350,325,300,275,250,225,200};
+    int adc_nr[32]={200,225,250,275,300,325,350,375,\
+		    0,25,50,75,100,125,150,175,\
+		    175,150,125,100,75,50,25,0,\
+		    375,350,325,300,275,250,225,200};
     int row, col;
 
     int isample;
@@ -126,18 +129,24 @@ class moench03CtbData : public slsDetectorData<uint16_t> {
   */
     virtual char *readNextFrame(ifstream &filebin){
       //	int afifo_length=0;  
- 	uint16_t *afifo_cont; 
+      uint16_t *afifo_cont; 
-    if (filebin.is_open()) {
+      int ib=0;
+      if (filebin.is_open()) {
 	afifo_cont=new uint16_t[dataSize/2];
- 	if (filebin.read((char*)afifo_cont,dataSize)) {
+ 	while (filebin.read(((char*)afifo_cont)+ib,2)) {
+	  ib+=2;
+	  if (ib==dataSize) break;
+	}
+	if (ib>0) {
 	  iframe++;
+	  cout << ib << "-" << endl;
 	  return (char*)afifo_cont;
 	} else {
-	delete [] afifo_cont;
+	  delete [] afifo_cont;
-	return NULL;
+	  return NULL;
 	}
-    }     
+      }     
-    return NULL;
+      return NULL;
     };
 
 

slsDetectorCalibration/moench03ReadData.C

@@ -15,7 +15,7 @@
 //#include <queue>
 #include <fstream>
 #include "moench03CtbData.h"
-#include "moenchCommonMode.h"
+#include "moench03CommonMode.h"
 #define MYROOT1
 #include "singlePhotonDetector.h"
 
@@ -28,22 +28,26 @@ using namespace std;
 
 
 #define MY_DEBUG 1
+
 #ifdef MY_DEBUG
 #include <TCanvas.h>
 #endif
 
 
 
-TH2F *readImage(ifstream &filebin, int dsize=5000*32) {
+TH2F *readImage(ifstream &filebin,   TH2F *h2=NULL, TH2F *hped=NULL) {
   moench03CtbData *decoder=new moench03CtbData();
   char *buff=decoder->readNextFrame(filebin);
-  TH2F *h2=NULL;
+  
 
 //   TH1F *h1=new TH1F("h1","",400*400,0,400*400);
 //   int ip=0;
   if (buff) {
-    h2=new TH2F("h2","",400,0,400,400,0,400);
+    if (h2==NULL) {
-    h2->SetStats(kFALSE);
+      h2=new TH2F("h2","",400,0,400,400,0,400);
+      h2->SetStats(kFALSE);
+    }
+    cout << "." << endl;
     for (int ix=0; ix<400; ix++) {
       for (int iy=0; iy<400; iy++) {
 	//	cout <<  decoder->getDataSize() << " " << decoder->getValue(buff,ix,iy)<< endl;
@@ -51,23 +55,37 @@ TH2F *readImage(ifstream &filebin, int dsize=5000*32) {
 	//	h1->SetBinContent(++ip,decoder->getValue(buff,ix,iy));
       }
     }
+    if (hped) h2->Add(hped,-1);
+    return h2;
   }
-  return h2;
+  return NULL;
 }
 
 
 TH2F *readImage(char *fname, int iframe=0, TH2F *hped=NULL) {
   ifstream filebin;
   filebin.open((const char *)(fname), ios::in | ios::binary);
-  TH2F *h2=readImage(filebin);
+  TH2F *h2=new TH2F("h2","",400,0,400,400,0,400);
-  for (int i=0; i<iframe; i++) {
+  TH2F *hh;
-    if (h2) delete h2; 
+  hh=readImage(filebin, h2, hped );
-    else break;
+  if (hh==NULL) {
-    h2=readImage(filebin);
+
-    cout << "="<< endl;
+    delete h2;
+    return NULL;
   }
-  filebin.close();
+  for (int i=0; i<iframe; i++) {
-  if (h2!=NULL && hped!=NULL)
+    if (hh==NULL) break;
+    hh=readImage(filebin, h2, hped );
+    if (hh)
+      cout << "="<< endl;
+    else {
+      delete h2;
+      return NULL;
+    }
+  }
+  if (filebin.is_open())
+    filebin.close();
+  if (hped!=NULL)
     h2->Add(hped,-1);
 
   return h2;
@@ -96,6 +114,7 @@ TH2F *calcPedestal(char *fformat, int runmin, int runmax){
 	}
       }
       delete [] buff;
+      cout << "="<< endl;
       ii++;
     }
     if (filebin.is_open())
@@ -140,13 +159,22 @@ Loops over data file to find single photons, fills the tree (and writes it to fi
 THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1) {
   
   moench03CtbData *decoder=new moench03CtbData();
-   moenchCommonMode *cmSub=NULL;
+  cout << "decoder allocated " << endl;
-//   if (cmsub)
+
-//     cmSub=new moenchCommonMode();
+   moench03CommonMode *cmSub=NULL;
+   if (cmsub) {
+     cmSub=new moench03CommonMode(100);
+     cout << "common mode allocated " << endl;
+
+   } else {
+     
+     cout << "non allocating common mode  " << endl;
+   }
    int iev=0;
   int nph=0;
 
-  singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub);
+  singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 10, 100);
+  cout << "filter allocated " << endl;
 
   char *buff;
   char fname[10000];
@@ -159,29 +187,36 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
 
   // data=decoder->getCluster();
 
-  TH2F *h2;
-  TH2F *h3;
-  TH2F *hetaX;
-  TH2F *hetaY;
 
   THStack *hs=new THStack("hs",fformat);
 
+  cout << "hstack allocated " << endl;
 
 
   TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
   hs->Add(h1);
+  cout << "h1 allocated " << endl;
 
+  TH2F *h2;
+  TH2F *h3;
   if (hitfinder) {
     h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
+  cout << "h2 allocated " << endl;
     h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
-    hetaX=new TH2F("hetaX",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
+  cout << "h3 allocated " << endl;
-    hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
+    // hetaX=new TH2F("hetaX",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
+    //  hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
     hs->Add(h2);
     hs->Add(h3);
-    hs->Add(hetaX);
+    // hs->Add(hetaX);
-    hs->Add(hetaY);
+    // hs->Add(hetaY);
   }
-
+  if (hs->GetHists()) {
+    for (int i=0; i<3; i++)
+      if (hs->GetHists()->At(1)) cout << i << " " ; 
+    cout << " histos allocated " << endl;
+  } else
+    cout << "no hists in stack " << endl;
   
   
   ifstream filebin;
@@ -199,6 +234,8 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
   
 
 #ifdef MY_DEBUG
+
+  cout << "debug mode " << endl;
   
   TCanvas *myC;
   TH2F *he;
@@ -207,21 +244,22 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
   TCanvas *cH3;
 
   if (hitfinder) {
-    myC=new TCanvas();
+    myC=new TCanvas("myc");
     he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
     he->SetStats(kFALSE);
     he->Draw("colz");
-    cH1=new TCanvas();
+    cH1=new TCanvas("ch1");
     cH1->SetLogz();
     h1->Draw("colz");
-    cH2=new TCanvas();
+    cH2=new TCanvas("ch2");
     cH2->SetLogz();
     h2->Draw("colz");
-    cH3=new TCanvas();
+    cH3=new TCanvas("ch3");
     cH3->SetLogz();
     h3->Draw("colz");
   }
 #endif
+
   filter->newDataSet();
 
 
@@ -255,12 +293,12 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
 
 #ifdef MY_DEBUG
 	  if (hitfinder) {
-	    if (iev%1000==0)
+	    //  if (iev%10==0)
 	      he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
 	  }
 #endif	  
 	  
-	  if (nf>1000) {
+	  //  if (nf>10) {
 	    h1->Fill(filter->getClusterTotal(1), iy+NR*ix);
 	    if (hitfinder) {
 	      
@@ -280,21 +318,22 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
 	    }
 	    
 
-	  }
+	    // }
 	}
 	  //////////////////////////////////////////////////////////
     
 #ifdef MY_DEBUG
-      if (iev%1000==0) {
+      cout << iev << " " << h1->GetEntries() << " " << h2->GetEntries() << endl;
-	myC->Modified();
+//       if (iev%10==0) {
-	myC->Update();
+// 	myC->Modified();
-	cH1->Modified();
+// 	myC->Update();
-	cH1->Update();
+// 	cH1->Modified();
-	cH2->Modified();
+// 	cH1->Update();
-	cH2->Update();
+// 	cH2->Modified();
-	cH3->Modified();
+// 	cH2->Update();
-	cH3->Update();
+// 	cH3->Modified();
-      }
+// 	cH3->Update();
+//       }
       iev++;
 #endif		    
       nf++;
@@ -311,7 +350,18 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
   if (hitfinder)
     tall->Write(tall->GetName(),TObject::kOverwrite);
   
+   //////////////////////////////////////////////////////////
     
+#ifdef MY_DEBUG  
+	myC->Modified();
+	myC->Update();
+	cH1->Modified();
+	cH1->Update();
+	cH2->Modified();
+	cH2->Update();
+	cH3->Modified();
+	cH3->Update();
+#endif	
   
   delete decoder;
   cout << "Read " << nf << " frames" << endl;

slsDetectorCalibration/readMoench03Data.C

@@ -0,0 +1,143 @@
+#include "moench03ReadData.C"
+
+
+
+void readMoench03Data(char *tit, int sign=1){
+
+
+
+
+
+  char fname[1000];
+  TFile *fout;
+  THStack *hs;
+
+  sprintf(fname,"/scratch/roberto/photons.root");
+  fout=new TFile(fname,"RECREATE");
+
+  sprintf(fname,"/scratch/roberto/run_%%d.raw");
+
+  hs=moench03ReadData(fname,"photons",136,1135,1500,-500,2500,1,0.,1,399,1,399, 0,1); 
+//   cout << "returned" << endl;
+//   hs->SetName(tit);
+//   hs->SetTitle(tit);
+//   cout << "name/title set" << endl;
+  
+  
+//   if (hs->GetHists()) {
+//     for (int i=0; i<3; i++) {
+//       if (hs->GetHists()->At(i)) {
+// 	cout << i << " " ;
+// 	(TH2F*)(hs->GetHists()->At(i))->Write();
+//       }
+//     } 
+//     cout << " histos written " << endl;
+//   } else
+//     cout << "no hists in stack " << endl;
+  
+
+//   fout->Close();
+
+
+
+}
+
+
+
+// void raedNoiseDataN(char *tit, int sign=1){
+
+
+
+//   char fname[1000];
+//   char f[1000];
+//   TFile *fout;
+//   THStack *hs2N;
+
+//   sprintf(fname,"/data/moench_xbox_20140116/noise_%s.root",tit);
+//   fout=new TFile(fname,"RECREATE");
+
+//   sprintf(fname,"/data/moench_xbox_20140116/noise_%s_f00000%%04d000_0.raw",tit);
+
+//   hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0,0); 
+//   hs2N->SetName(tit);
+//   hs2N->SetTitle(tit);
+//   (TH2F*)(hs2N->GetHists()->At(0))->Write();
+
+//   // (TH2F*)(hs2N->GetHists()->At(1))->Write();
+//   // (TH2F*)(hs2N->GetHists()->At(2))->Write();
+//   // (TH2F*)(hs2N->GetHists()->At(3))->Write();
+//   // (TH2F*)(hs2N->GetHists()->At(4))->Write();
+
+
+//   fout->Close();
+
+
+
+// }
+
+
+
+// void g4() {
+
+//   raedNoiseData("cds_g4_low_gain");
+//   raedNoiseData("cds_g4_sto1_only");
+//   raedNoiseData("cds_g4_no sto");
+  
+
+
+// }
+
+// void no_cds() {
+
+//   raedNoiseData("cds_disable_low_gain",-1);
+//   raedNoiseData("cds_disable_sto1_only",-1);
+//   raedNoiseData("cds_disable_no sto",-1);
+  
+
+
+// }
+
+// void all_gains() {
+
+//   raedNoiseData("cds_g2");
+//   raedNoiseData("cds_g2HC");
+//   raedNoiseData("cds_g1_2");
+//   raedNoiseData("cds_g2_3");
+  
+
+
+// }
+
+// void all_low_gains() {
+
+//   raedNoiseData("cds_g2_low_gain");
+//   raedNoiseData("cds_g2HC_low_gain");
+//   raedNoiseData("cds_g1_2_low_gain");
+//   raedNoiseData("cds_g2_3_low_gain");
+// }
+
+// /*
+// clkdivider data
+// /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv17_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 12:40 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv25_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 13:26 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv35_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 14:09 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv50_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 14:54 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv70_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 16:42 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv110_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 17:27 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv170_f000000010000_0.raw
+// */
+
+
+// /* oversampled data
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 18:12 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_0.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 18:47 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_1.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 19:22 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_2.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 20:02 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_3.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 20:41 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_4.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 21:16 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_5.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 21:56 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_6.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 22:35 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_7.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 23:11 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_8.raw
+// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 23:50 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_9.raw
+// */
+

slsDetectorCalibration/singlePhotonDetector.h

@@ -109,7 +109,7 @@ class singlePhotonDetector {
 	\param cm commonModeSubtraction algorithm to be used (NULL unsets) 
 	\returns pointer to the actual common mode subtraction algorithm
     */
-    commonModeSubtraction setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
+    commonModeSubtraction *setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
 
 
     /**