diff --git a/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h b/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h new file mode 100644 index 000000000..7623bc7e9 --- /dev/null +++ b/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h @@ -0,0 +1,202 @@ +#ifndef ETA_INTERPOLATION_ADAPTIVEBINS_H +#define ETA_INTERPOLATION_ADAPTIVEBINS_H + + +#include "tiffIO.h" +#include "etaInterpolationBase.h" + +class etaInterpolationAdaptiveBins : public etaInterpolationBase{ + private: + double calcDiff(double avg, double *hx, double *hy) { + double p_tot=0; + double diff=0; + double bsize=1./nSubPixels; + + for (int ipx=0; ipx=nSubPixels; ipx++) { + for (int ipy=0; ipy=nSubPixels; ipy++) { + p_tot=0; + for (int ibx=0; ibx=((ipx)*bsize) && hx[ibx+iby*nbeta]<((ipx+1)*bsize) && hy[ibx+iby*nbeta]>=((ipy)*bsize) && hy[ibx+iby*nbeta]<((ipy+1)*bsize)) { + p_tot+=heta[ibx+iby*nbeta]; + } + } + } + + diff+=(p_tot-avg)*(p_tot-avg); + + } + } + return diff; + + } + + void iterate(double *newhhx, double *newhhy) { + + double bsize=1./nSubPixels; + + double hy[nbeta]; //profile y + double hx[nbeta]; //profile x + double hix[nbeta]; //integral of projection x + double hiy[nbeta]; //integral of projection y + + double tot_eta_x=0; + double tot_eta_y=0; + for (int ipy=0; ipy=nSubPixels; ipy++) { + + + for (int ibx=0; ibx=((ipy)*bsize) && hhy[ibx+iby*nbeta]<((ipy+1)*bsize)) { + hx[ibx]+=heta[ibx+iby*nbeta]; + } + if ( hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<((ipy+1)*bsize)) { + hy[iby]+=heta[ibx+iby*nbeta]; + } + } + } + + + } + + hix[0]=hx[0]; + hiy[0]=hy[0]; + for (int ib=1; ib=((ipy)*bsize) && hhy[ibx+iby*nbeta]<((ipy+1)*bsize)) { + newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x); + } + if ( hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<((ipy+1)*bsize)) { + newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y); + } + } + } + } + + + } + + + public: + etaInterpolationAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){}; + + etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationBase(orig){}; + + virtual etaInterpolationAdaptiveBins* Clone() { + + return new etaInterpolationAdaptiveBins(this); + + }; + + + + virtual void prepareInterpolation(int &ok) + { + ok=1; + + + ///*Eta Distribution Rebinning*/// + double bsize=1./nSubPixels; //precision + // cout<<"nPixelsX = "< -ClassImp(EtaVEL); +// ClassImp(EtaVEL); -double Median(const TH1D * histo1) { +// double Median(const TH1D * histo1) { - int numBins = histo1->GetXaxis()->GetNbins(); - Double_t *x = new Double_t[numBins]; - Double_t* y = new Double_t[numBins]; - for (int i = 0; i < numBins; i++) { - x[i] = histo1->GetBinCenter(i); - y[i] = histo1->GetBinContent(i); - } - return TMath::Median(numBins, x, y); -} +// int numBins = histo1->GetXaxis()->GetNbins(); +// Double_t *x = new Double_t[numBins]; +// Double_t* y = new Double_t[numBins]; +// for (int i = 0; i < numBins; i++) { +// x[i] = histo1->GetBinCenter(i); +// y[i] = histo1->GetBinContent(i); +// } +// return TMath::Median(numBins, x, y); +// } double *EtaVEL::getPixelCorners(int x, int y){ diff --git a/slsDetectorCalibration/interpolations/etaVEL/interpolation_etaVEL.cpp b/slsDetectorCalibration/interpolations/etaVEL/interpolation_etaVEL.cpp deleted file mode 100644 index 1cdd3c718..000000000 --- a/slsDetectorCalibration/interpolations/etaVEL/interpolation_etaVEL.cpp +++ /dev/null @@ -1,134 +0,0 @@ -#include "interpolation_EtaVEL.h" -#include "TH2F.h" -#include "TCanvas.h" -#include "TROOT.h" -//#include "EtaVEL.h" -#include "EtaVEL.cpp" -/* -Zum erstellen der correction map ist createGainAndEtaFile(...) in EVELAlg.C der entry point. -Zum erstellen des HR images ist createImage(...) der entry point. -*/ -interpolation_EtaVEL::interpolation_EtaVEL(int nx, int ny, int ns, double etamin, double etamax, int p) : slsInterpolation(nx, ny, ns), newEta(NULL), heta(NULL), plot(p) { - newEta = new EtaVEL(nSubPixels,etamin,etamax,nPixelsX, nPixelsY); - heta= new TH2F("heta","heta",50*nSubPixels, etamin,etamax,50*nSubPixels, etamin,etamax); - heta->SetStats(kFALSE); -} - -interpolation_EtaVEL::~interpolation_EtaVEL() { - delete newEta; - delete heta; -} - - -void interpolation_EtaVEL::prepareInterpolation(int &ok, int maxit) { - int nit=0; - while ((newEta->converged != 1) && nit++Modified(); - gPad->Update(); - } - if (newEta->converged==1) ok=1; else ok=0; -} - -int interpolation_EtaVEL::addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay) { - Double_t sum, totquad, sDum[2][2]; - int corner =calcEta(cluster, etax, etay, sum, totquad, sDum); - //check if it's OK...should redo it every time? - //or should we fill a finer histogram and afterwards re-fill the newEta? - addToFlatField(etax, etay); - return corner; -} - -int interpolation_EtaVEL::addToFlatField(Double_t etax, Double_t etay) { - // newEta->fill(etaX,etaY); - heta->Fill(etax,etay); - return 0; -} - -void interpolation_EtaVEL::iterate() { - cout << " -------------- newEta refilled"<< endl; - for (int ibx=0; ibxGetNbinsX(); ibx++) { - for (int iby=0; ibyGetNbinsY(); iby++) { - newEta->fill(heta->GetXaxis()->GetBinCenter(ibx+1),heta->GetYaxis()->GetBinCenter(iby+1),heta->GetBinContent(ibx+1,iby+1)); - } - } - newEta->updatePixelPos(); - cout << " -------------- pixelPosition updated"<< endl; -} - -void interpolation_EtaVEL::DrawH() { - heta->Draw("col"); - (newEta->plotPixelBorder())->Draw(); -} - - -void interpolation_EtaVEL::getInterpolatedPosition(Int_t x, Int_t y, Double_t *cluster, Double_t &int_x, Double_t &int_y) { - - Double_t etax, etay, sum, totquad, sDum[2][2]; - - int corner =calcEta(cluster, etax, etay, sum, totquad, sDum); - - int bin = newEta->findBin(etax,etay); - if (bin<=0) { - int_x=-1; - int_y=-1; - return; - } - double subX = ((double)(newEta->getXBin(bin))+.5)/((double)newEta->getNPixels()); - double subY = ((double)(newEta->getYBin(bin))+.5)/((double)newEta->getNPixels()); - - double dX, dY; - switch (corner) { - case TOP_LEFT: - dX=-1.; - dY=+1.; - break; - case TOP_RIGHT: - dX=+1.; - dY=+1.; - break; - case BOTTOM_LEFT: - dX=-1.; - dY=-1.; - break; - case BOTTOM_RIGHT: - dX=+1.; - dY=-1.; - break; - default: - dX=0; - dY=0; - } - - int_x=((double)x)+ subX+0.5*dX; - int_y=((double)y)+ subY+0.5*dY; - - // cout << corner << " " << subX<< " " << subY << " " << dX << " " << dY << " " << int_x << " " << int_y << endl; - -}; - - -// void interpolation_EtaVEL::Streamer(TBuffer &b){newEta->Streamer(b);}; -void interpolation_EtaVEL::getInterpolatedBin(Double_t *cluster, Int_t &int_x, Int_t &int_y) { - - Double_t etax, etay, sum, totquad, sDum[2][2]; - - int corner =calcEta(cluster, etax, etay, sum, totquad, sDum); - - int bin = newEta->findBin(etax,etay); - if (bin<0) { - int_x=-1; - int_y=-1; - return; - } - int_x=newEta->getXBin(bin); - int_y=newEta->getYBin(bin); - - - -}; - diff --git a/slsDetectorCalibration/interpolations/interpolation_EtaVEL.h b/slsDetectorCalibration/interpolations/interpolation_EtaVEL.h deleted file mode 100644 index 93e2eaddd..000000000 --- a/slsDetectorCalibration/interpolations/interpolation_EtaVEL.h +++ /dev/null @@ -1,52 +0,0 @@ -#ifndef INTERPOLATION_ETAVEL_H -#define INTERPOLATION_ETAVEL_H - -#include -#include "EtaVEL.h" -#include "TH2F.h" -//#include "EtaVEL.cpp" -//class EtaVEL; - -class interpolation_EtaVEL: public slsInterpolation { - - public: - interpolation_EtaVEL(int nx=40, int ny=160, int ns=25, double etamin=-0.02, double etamax=1.02, int p=0); - ~interpolation_EtaVEL(); - - - //create eta distribution, eta rebinnining etc. - //returns flat field image - void prepareInterpolation(int &ok){prepareInterpolation(ok,10000);}; - void prepareInterpolation(int &ok, int maxit); - - //create interpolated image - //returns interpolated image - - //return position inside the pixel for the given photon - void getInterpolatedPosition(Int_t x, Int_t y, Double_t *data, Double_t &int_x, Double_t &int_y); - void getInterpolatedBin(Double_t *cluster, Int_t &int_x, Int_t &int_y); - - - - int addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay); - int addToFlatField(Double_t etax, Double_t etay); - int setPlot(int p=-1) {if (p>=0) plot=p; return plot;}; - int WriteH(){newEta->Write("newEta"); heta->Write("heta");}; - EtaVEL *setEta(EtaVEL *ev){if (ev) {delete newEta; newEta=ev;} return newEta;}; - TH2F *setEta(TH2F *ev){if (ev) {delete heta; heta=ev;} return heta;}; - void iterate(); - void DrawH(); - double getChiSq(){return newEta->getChiSq();}; - - - - protected: - EtaVEL *newEta; - TH2F *heta; - int plot; - - // ClassDefNV(interpolation_EtaVEL,1); - // #pragma link C++ class interpolation_EtaVEL-; -}; - -#endif diff --git a/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp b/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp index 6ce7380d5..341f27994 100644 --- a/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp +++ b/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp @@ -42,7 +42,7 @@ int main(int argc, char *argv[]) { int ok; int iprog=0; - moench03T1ReceiverData *decoder=new moench03T1ReceiverData(); + moench03T1ReceiverDataNew *decoder=new moench03T1ReceiverDataNew(); //moench03T1ZmqData *decoder=new moench03T1ZmqData(); singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200); // char tit[10000]; diff --git a/slsDetectorCalibration/moenchExecutables/readClusters.C b/slsDetectorCalibration/moenchExecutables/readClusters.C index f030e91eb..d75338508 100644 --- a/slsDetectorCalibration/moenchExecutables/readClusters.C +++ b/slsDetectorCalibration/moenchExecutables/readClusters.C @@ -13,23 +13,23 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { h2=new TH2F("h2",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5); //h2mult=new TH2F("h2mult",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5); // TH2F *hint=new TH2F("hint",fname,nx*ns, -0.5, nx-0.5, ny*ns, -0.5, ny-0.5); - TH1F *hf=new TH1F("hf","hf",1000,0,30000); + TH1F *hf=new TH1F("hf","hf",1000,0,3000000); //TH2F *hff=new TH2F("hff","hff",ns, -0.5, 0.5, ns, -0.5, +0.5); TH1F *hsp=new TH1F("hsp",fname,500,0,2000); TH1F *hsp1=new TH1F("hsp1",fname,500,0,2000); - TH1F *hsp2=new TH1F("hsp2",fname,500,0,2000); - TH1F *hsp3=new TH1F("hsp3",fname,500,0,2000); - hsp1->SetLineColor(2); - hsp2->SetLineColor(3); - hsp3->SetLineColor(4); + // TH1F *hsp2=new TH1F("hsp2",fname,500,0,1000); + // TH1F *hsp3=new TH1F("hsp3",fname,500,0,1000); + hsp1->SetLineColor(2); + // hsp2->SetLineColor(3); + // hsp3->SetLineColor(4); TCanvas *c=new TCanvas(); c->SetLogz(kTRUE); h2->Draw("colz"); TCanvas *c1=new TCanvas(); hsp->Draw(); hsp1->Draw("same"); - hsp2->Draw("same"); - hsp3->Draw("same"); + // hsp2->Draw("same"); + // hsp3->Draw("same"); TCanvas *c2=new TCanvas(); hf->Draw(); // hint->Draw("colz"); @@ -75,26 +75,34 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { //max at 340 //if (tot>200) { w=1; + + if (qt/tot>0.8 && qt/tot<1.2){ if (f0<0) f0=cl.iframe; hf->Fill(cl.iframe-f0); + if (qt>540) w++; + if (qt>820) w++; //(tot+3.5*phs)/phw; //} else //w=0; // if (w) { - // if (qt/tot>0.8 && qt/tot<1.2){ - if (cl.y<350) { - if (cl.y<100 || cl.y>300) { - hsp->Fill(qt); - hsp2->Fill(cl.get_data(0,0)); - } else { - hsp1->Fill(qt); - hsp3->Fill(cl.get_data(0,0)); - } + // if (cl.y<350) { + // if (cl.y<100 || cl.y>300) { + if (cl.x>150 && cl.x<250 && cl.y>200 && cl.y<250) + hsp1->Fill(qt); + else + hsp->Fill(qt); + + // hsp2->Fill(cl.get_data(0,0)); + // } else { + // hsp1->Fill(qt); + // hsp3->Fill(cl.get_data(0,0)); + // } //if (cl.x>160 && cl.x<260 && cl.y>30 && cl.y<80 && tot>0) //if (w==1) { - h2->Fill(cl.x,cl.y,w); - } + // if (w==2) + h2->Fill(cl.x,cl.y,w); + // } //} // hint->Fill(px+cl.x,py+cl.y,w); // if (cl.y<350) @@ -104,7 +112,7 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { //h2mult->Fill(cl.x,cl.y,w); - //} + } iph+=w; if (iph%100000==0) { c->Modified(); @@ -139,7 +147,7 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { return h2; - } else +} else cout << "could not open file " << fname << endl; return NULL;