diff --git a/slsDetectorCalibration/moench02ModuleData.h b/slsDetectorCalibration/moench02ModuleData.h index 7dc2d1bdf..ca2e753c4 100644 --- a/slsDetectorCalibration/moench02ModuleData.h +++ b/slsDetectorCalibration/moench02ModuleData.h @@ -8,6 +8,7 @@ #include #include +#include using namespace std; @@ -110,7 +111,7 @@ class moench02ModuleData : public slsDetectorData { char *readNextFrame(ifstream &filebin) { - if (oldbuff) + if (oldbuff) delete [] oldbuff; oldbuff=buff; @@ -159,7 +160,7 @@ class moench02ModuleData : public slsDetectorData { if (np>0) cout << "Too few packets for this frame! "<< fnum << " " << pnum << endl; delete [] data; - buff=NULL; + data=NULL; } buff=data; return data; @@ -183,7 +184,7 @@ class moench02ModuleData : public slsDetectorData { if (ix>0 && hc!=0) val-=hc*slsDetectorData::getChannelShort(buff,ix-1,iy); - if (oldbuff && tc!=0) + if (oldbuff && tc!=0) val+=tc*slsDetectorData::getChannelShort(oldbuff,ix,iy); @@ -248,7 +249,19 @@ class moench02ModuleData : public slsDetectorData { double *getCluster(){return &cluster[0][0];}; + TTree *initMoenchTree(char *tname, Int_t *iframe, Int_t * x, Int_t* y, Double_t data[3][3], Double_t *ped, Double_t *sigma) { + TTree* tall=new TTree(tname,tname); + tall->Branch("iFrame",iframe,"iframe/I"); + tall->Branch("x",x,"x/I"); + tall->Branch("y",y,"y/I"); + tall->Branch("data",data,"data[3][3]/D"); + tall->Branch("pedestal",ped,"pedestal/D"); + tall->Branch("rms",sigma,"rms/D"); + return tall; + }; + + private: diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index 078d5931a..4a16ef9ce 100644 --- a/slsDetectorCalibration/moenchReadData.C +++ b/slsDetectorCalibration/moenchReadData.C @@ -23,7 +23,8 @@ using namespace std; -TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, int nsigma=5, double hc=0, double tc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR) { + +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) { moench02ModuleData *decoder=new moench02ModuleData(); char *buff; @@ -35,21 +36,52 @@ TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int moench02ModuleData::eventType thisEvent=moench02ModuleData::PEDESTAL; + // int iframe; + // double *data, ped, sigma; - TH2F *h2=NULL; - h2=new TH2F("h2",fformat,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5); + // data=decoder->getCluster(); - int val, dum; - double me=0, sig, tot, vv; + TH2F *h2; + TH2F *h3; + TH2F *hetaX; + TH2F *hetaY; - MovingStat stat[160][160]; + THStack *hs=new THStack("hs",fformat); + + TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5); + hs->Add(h1); + + if (pedsub) { + + h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5); + 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); + 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(hetaY); + } + + ifstream filebin; - int nbg=1000; int ix=20, iy=20, ir, ic; + + Int_t iFrame, x, y; + Double_t data[3][3],ped,sigma; + + TTree *tall; + if (pedsub) { + cout << "Subtracting pedestal!!!! " << endl; + tall=decoder->initMoenchTree(tit, &iFrame, &x, &y, data, &ped, &sigma); + } + double tot, tl, tr, bl, br, v; + + // 6% x-talk from previous pixel // 12% x-talk from previous frame @@ -63,14 +95,16 @@ TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int for (ix=xmin-1; ix100) { thisEvent= decoder->getEventType(ix, iy, hc, tc, 1); } + + if (thisEvent==moench02ModuleData::PEDESTAL) { - decoder->addToPedestal(decoder->getClusterElement(0,0), ix, iy); + decoder->addToPedestal( decoder->getChannelShort(ix, iy, hc, tc), ix, iy); } @@ -78,6 +112,101 @@ TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int Add here the function that you want to call: fill histos, make trees etc. ************************************************************/ // getClusterElement to access quickly the photon and the neighbours + if (nf>1000) { + tot=0; + tl=0; + tr=0; + bl=0; + br=0; + h1->Fill(decoder->getClusterElement(0,0), 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); + data[ic+1][ir+1]=v; + tot+=v; + if (ir<1) { + + if (ic<1) { + bl+=v; + + } + if (ic>-1) { + br+=v; + } + } + + if (ir>-1) { + if (ic<1) { + tl+=v; + } + if (ic>-1) { + tr+=v; + } + } + } + } + + + 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); + } + } 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); + } + } 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); + } + } 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); + } + + + + } + + h3->Fill(tot, iy+NR*ix); + iFrame=decoder->getFrameNumber(buff); + ped=decoder->getPedestal(ix,iy); + sigma=decoder->getPedestalRMS(ix,iy); + x=ix; + y=iy; + + tall->Fill(); + + } + + + + } + } else { + //cout << ix << " " << iy << endl; + h1->Fill(decoder->getChannelShort(ix, iy), iy+NR*ix); + } + ////////////////////////////////////////////////////////// } cout << "=" ; @@ -89,13 +218,13 @@ Add here the function that you want to call: fill histos, make trees etc. else cout << "could not open file " << fname << endl; } - - + if (pedsub) + tall->Write(); delete decoder; cout << "Read " << nf << " frames" << endl; - return h2; + return hs; }