diff --git a/slsDetectorCalibration/moench02ModuleData.h b/slsDetectorCalibration/moench02ModuleData.h index ca2e753c4..7dd718f8f 100644 --- a/slsDetectorCalibration/moench02ModuleData.h +++ b/slsDetectorCalibration/moench02ModuleData.h @@ -109,8 +109,55 @@ class moench02ModuleData : public slsDetectorData { }; + void calculateCommonMode(int xmin=0, int xmax=160, int ymin=0, int ymax=160, double hc=0, double tc=0) { + + int scmin=xmin/40; + int scmax=(xmax-1)/40+1; + int isc=0, ix, iy; + + for (isc=scmin; isc0) + cmStat[isc].Calc(cmPed[isc]/nCm[isc]); + } + + + }; + + double getCommonMode(int ix, int iy) { + int isc=ix/40; + if (nCm[isc]>0) + return cmPed[isc]/nCm[isc]-cmStat[isc].Mean(); + else + return 0; + }; + + + char *readNextFrame(ifstream &filebin) { + if (oldbuff) delete [] oldbuff; @@ -223,10 +270,7 @@ class moench02ModuleData : public slsDetectorData { } else if (v<-nSigma*getPedestalRMS(ix,iy)) ret=NEGATIVE_PEDESTAL; } - - } - } } @@ -238,14 +282,19 @@ class moench02ModuleData : public slsDetectorData { } - + cx=ix; + cy=iy; return ret; }; - double getClusterElement(int ic, int ir){return cluster[ic+1][ir+1];}; + double getClusterElement(int ic, int ir, int cm=0){ + if (cm) cluster[ic+1][ir+1]-getCommonMode(cx,cy); + else return cluster[ic+1][ir+1]; + }; + double *getCluster(){return &cluster[0][0];}; @@ -271,8 +320,12 @@ class moench02ModuleData : public slsDetectorData { char *buff; char *oldbuff; int pedSub; + int cmSub; - + MovingStat cmStat[4]; //stores the common mode average per supercolumn + double cmPed[4]; // stores the common mode for this frame per supercolumn + double nCm[4]; // stores the number of pedestals to calculate the cm per supercolumn + int cx, cy; }; diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index bf9313830..c2ebd2efe 100644 --- a/slsDetectorCalibration/moenchReadData.C +++ b/slsDetectorCalibration/moenchReadData.C @@ -16,15 +16,35 @@ #include #include "moench02ModuleData.h" +//#include "MovingStat.h" + using namespace std; #define NC 160 #define NR 160 +/** + char *fformat, file name format + char *tit, title of the tree etc. + int runmin, minimum run number + int runmax, max run number + int nbins=1500, number of bins for spectrum hists + int hmin=-500, minimum for spectrum hists + int hmax=1000, maximum for spectrum hists + int sign=1, sign of the spectrum to find hits + double hc=0, readout correlation coefficient with previous pixel + double tc=0, time correlation coefficient with previous frame (case of bad reset) + int xmin=0, minimum x coordinate + int xmax=NC, maximum x coordinate + int ymin=0, minimum y coordinate + int ymax=NR, maximum y coordinate + int pedsub=1, enable pedestal subtraction + int cmsub=0 enable commonmode subtraction + +*/ - -THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, double tc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR, int pedsub=1) { +THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, double tc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR, int pedsub=1, int cmsub=0) { moench02ModuleData *decoder=new moench02ModuleData(); char *buff; @@ -32,8 +52,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb char fname[10000]; double oldval; int nf=0; - - + moench02ModuleData::eventType thisEvent=moench02ModuleData::PEDESTAL; // int iframe; @@ -80,7 +99,9 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb tall=decoder->initMoenchTree(tit, &iFrame, &x, &y, data, &ped, &sigma); } double tot, tl, tr, bl, br, v; + + int scmin=xmin/40, scmax=(xmax-1)/40+1; // 6% x-talk from previous pixel // 12% x-talk from previous frame @@ -92,20 +113,40 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb filebin.open((const char *)(fname), ios::in | ios::binary); while ((buff=decoder->readNextFrame(filebin))) { - + + + + if (nf>100) { + if (cmsub) { + for (int isc=scmin; isccalculateCommonMode(3+isc*40, 40*(isc-1)-3, 3, NR-3, hc, tc); + } + } + } + + + + for (ix=xmin-1; ix100) { - thisEvent= decoder->getEventType(ix, iy, hc, tc, 1); - } - + + + if (nf>100) { + thisEvent= decoder->getEventType(ix, iy, hc, tc, 1); + } + - - if (thisEvent==moench02ModuleData::PEDESTAL) { - decoder->addToPedestal( decoder->getChannelShort(ix, iy, hc, tc), ix, iy); - } + + + if (thisEvent==moench02ModuleData::PEDESTAL) { + if (cmsub && nf>1000) + decoder->addToPedestal( decoder->getChannelShort(ix, iy, hc, tc)- decoder->getCommonMode(ix,iy), ix, iy); + else + decoder->addToPedestal( decoder->getChannelShort(ix, iy, hc, tc), ix, iy); + + } /*********************************************************** @@ -118,14 +159,14 @@ Add here the function that you want to call: fill histos, make trees etc. tr=0; bl=0; br=0; - h1->Fill(decoder->getClusterElement(0,0), iy+NR*ix); + h1->Fill(decoder->getClusterElement(0,0, cmsub), iy+NR*ix); if (thisEvent==moench02ModuleData::PHOTON_MAX ) { for (ir=-1; ir<2; ir++) { for (ic=-1; ic<2; ic++) { - v=decoder->getClusterElement(ic,ir); + v=decoder->getClusterElement(ic,ir,cmsub); data[ic+1][ir+1]=v; tot+=v; if (ir<1) { @@ -154,34 +195,34 @@ Add here the function that you want to call: fill histos, make trees etc. if (bl>br && bl>tl && bl>tr) { h2->Fill(bl, iy+NR*ix); if (bl>0) { - hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,-1))/bl,iy+NR*ix); - hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(-1,0))/bl,iy+NR*ix); - hetaX->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(0,-1))/bl,iy+NR*(ix-1)); - hetaY->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(-1,0))/bl,(iy-1)+NR*ix); + hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/bl,iy+NR*ix); + hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/bl,iy+NR*ix); + hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/bl,iy+NR*(ix-1)); + hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/bl,(iy-1)+NR*ix); } } else if (br>bl && br>tl && br>tr) { h2->Fill(br, iy+NR*ix); if (br>0) { - hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,-1))/br,iy+NR*ix); - hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(1,0))/br,iy+NR*ix); - hetaX->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(0,-1))/br,iy+NR*(ix+1)); - hetaY->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(1,0))/br,iy-1+NR*ix); + hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/br,iy+NR*ix); + hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/br,iy+NR*ix); + hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/br,iy+NR*(ix+1)); + hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/br,iy-1+NR*ix); } } else if (tl>br && tl>bl && tl>tr) { h2->Fill(tl, iy+NR*ix); if (tl>0) { - hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,1))/tl,iy+NR*ix); - hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(-1,0))/tl,iy+NR*ix); - hetaX->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(0,1))/tl,iy+NR*(ix-1)); - hetaY->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(-1,0))/tl,iy+1+NR*ix); + hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tl,iy+NR*ix); + hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/tl,iy+NR*ix); + hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tl,iy+NR*(ix-1)); + hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/tl,iy+1+NR*ix); } } else if (tr>bl && tr>tl && tr>br) { h2->Fill(tr, iy+NR*ix); if (tr>0) { - hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,1))/tr,iy+NR*ix); - hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(1,0))/tr,iy+NR*ix); - hetaX->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(0,1))/tr,iy+NR*(ix+1)); - hetaY->Fill(1.-(decoder->getClusterElement(0,0)+decoder->getClusterElement(1,0))/tr,iy+1+NR*ix); + hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tr,iy+NR*ix); + hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/tr,iy+NR*ix); + hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tr,iy+NR*(ix+1)); + hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/tr,iy+1+NR*ix); } @@ -207,17 +248,17 @@ Add here the function that you want to call: fill histos, make trees etc. h1->Fill(decoder->getChannelShort(ix, iy), iy+NR*ix); } ////////////////////////////////////////////////////////// - } + cout << "=" ; nf++; - } + } cout << endl; if (filebin.is_open()) filebin.close(); else cout << "could not open file " << fname << endl; - } + } if (pedsub) tall->Write(tall->GetName(),TObject::kOverwrite); @@ -226,5 +267,5 @@ Add here the function that you want to call: fill histos, make trees etc. delete decoder; cout << "Read " << nf << " frames" << endl; return hs; -} + }