mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-04-22 22:40:02 +02:00
added moench03 T CPU data structure
This commit is contained in:
parent
ed7d6b2995
commit
37ba8f9a71
@ -470,6 +470,13 @@ TF1* energyCalibration::fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar)
|
||||
|
||||
|
||||
|
||||
TF1* energyCalibration::fitSpectrumPixel(TH1 *h1, Double_t *mypar, Double_t *emypar) {
|
||||
initFitFunction(fspixel,h1);
|
||||
return fitFunction(fspixel, h1, mypar, emypar);
|
||||
}
|
||||
|
||||
|
||||
|
||||
TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff) {
|
||||
|
||||
TGraphErrors *gr;
|
||||
|
||||
@ -317,6 +317,16 @@ class energyCalibration {
|
||||
TF1 *fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar);
|
||||
|
||||
|
||||
/**
|
||||
fits histogram with the spectrum
|
||||
\param h1 1d-histogram to be fitted
|
||||
\param mypar pointer to fit parameters array
|
||||
\param emypar pointer to fit parameter errors
|
||||
\returns the fitted function - can be used e.g. to get the Chi2 or similar
|
||||
*/
|
||||
TF1 *fitSpectrumPixel(TH1 *h1, Double_t *mypar, Double_t *emypar);
|
||||
|
||||
|
||||
/**
|
||||
calculates gain and offset for the set of inflection points
|
||||
\param nscan number of energy scans
|
||||
|
||||
678
slsDetectorCalibration/etaVEL/EVELAlg.C
Normal file
678
slsDetectorCalibration/etaVEL/EVELAlg.C
Normal file
@ -0,0 +1,678 @@
|
||||
#include <TH1D.h>
|
||||
#include <TH2D.h>
|
||||
#include <TPad.h>
|
||||
#include <TDirectory.h>
|
||||
#include <TEntryList.h>
|
||||
#include <TFile.h>
|
||||
#include <TMath.h>
|
||||
#include <TTree.h>
|
||||
#include <TChain.h>
|
||||
#include <THStack.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TF1.h>
|
||||
#include <TLegend.h>
|
||||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <deque>
|
||||
#include <list>
|
||||
#include <queue>
|
||||
#include <fstream>
|
||||
|
||||
#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.
|
||||
*/
|
||||
int etabins = 25;
|
||||
int nEtas = 25;
|
||||
Double_t dum[3][3];
|
||||
Int_t x,y,f,q;
|
||||
|
||||
int counter[5];
|
||||
int remoteCounter[5];
|
||||
|
||||
//TH2D *sum = new TH2D("sum","sum",3,-0.1,2.1,3,-0.1,2.1);
|
||||
//TH2F *subPos = new TH2F("subPos","subPos", 100, -1.,1. ,100, -1.,1.);
|
||||
TH2D *subPosAEta = new TH2D("subPosAEta","subPosAEta", 50, -.5,1.5 ,50, -.5,1.5);
|
||||
TH2D *subPosBEta = new TH2D("subPosBEta","subPosBEta", 50, -.5,1.5 ,50, -.5,1.5);
|
||||
|
||||
|
||||
|
||||
TH1D *cE = new TH1D("clusterEnergy","clusterEnergy",400, 0.,4000.);
|
||||
//TH1D *cES = new TH1D("clusterEnergyS","clusterEnergyS",400, 0.,4000.);
|
||||
|
||||
|
||||
TH2D *cES3vs2 = new TH2D("clusterEnergy3vs2","clusterEnergy3vs2",800, 0.,8000.,600,0.,6000.);
|
||||
TH2D *cES3vs2S = new TH2D("clusterEnergy3vs2S","clusterEnergy3vs2S",800, 0.,8000.,600,0.,6000.);
|
||||
|
||||
double th = 0.99;
|
||||
double sigmas = 1.0;
|
||||
|
||||
TH2D *imgRLR = new TH2D("imgRLR","imgRLR",160,0.0,160.0 ,160 ,0.0,160.0);
|
||||
TH2D *imgLR = new TH2D("imgLR","imgLR",160*2,0.0,160.0 ,160*2 ,0.0,160.0);
|
||||
|
||||
TH2D *clusHist= new TH2D("clusHist","clusHist",3,-0.5,2.5,3,-0.5,2.5);
|
||||
TH2D *clusHistC= new TH2D("clusHistC","clusHistC",3,-0.5,2.5,3,-0.5,2.5);
|
||||
|
||||
int **imgArray;
|
||||
|
||||
int findShape(Double_t cluster[3][3], double sDum[2][2]){
|
||||
int corner = -1;
|
||||
|
||||
double sum = cluster[0][0] + cluster[1][0] + cluster[2][0] + cluster[0][1] + cluster[1][1] + cluster[2][1] + cluster[0][2] + cluster[1][2] + cluster[2][2];
|
||||
|
||||
double sumTL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL
|
||||
double sumTR = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL
|
||||
double sumBL = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR
|
||||
double sumBR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR
|
||||
double sumMax = 0;
|
||||
|
||||
|
||||
//double **sDum = subCluster;
|
||||
Double_t ssDum[2][2];
|
||||
|
||||
// if(sumTL >= sumMax){
|
||||
sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0];
|
||||
sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1];
|
||||
|
||||
ssDum[0][0] = cluster[0][0]; ssDum[1][0] = cluster[0][1];
|
||||
ssDum[0][1] = cluster[1][0]; ssDum[1][1] = cluster[1][1];
|
||||
|
||||
corner = 2;
|
||||
sumMax=sumTL;
|
||||
// }
|
||||
|
||||
if(sumTR >= sumMax){
|
||||
sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0];
|
||||
sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1];
|
||||
|
||||
ssDum[0][0] = cluster[2][0]; ssDum[1][0] = cluster[2][1];
|
||||
ssDum[0][1] = cluster[1][0]; ssDum[1][1] = cluster[1][1];
|
||||
|
||||
corner = 0;
|
||||
sumMax=sumTR;
|
||||
}
|
||||
|
||||
if(sumBL >= sumMax){
|
||||
sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1];
|
||||
sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2];
|
||||
|
||||
ssDum[0][0] = cluster[0][2]; ssDum[1][0] = cluster[0][1];
|
||||
ssDum[0][1] = cluster[1][2]; ssDum[1][1] = cluster[1][1];
|
||||
|
||||
corner = 3;
|
||||
sumMax=sumBL;
|
||||
}
|
||||
|
||||
if(sumBR >= sumMax){
|
||||
sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1];
|
||||
sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2];
|
||||
|
||||
ssDum[0][0] = cluster[2][2]; ssDum[1][0] = cluster[2][1];
|
||||
ssDum[0][1] = cluster[1][2]; ssDum[1][1] = cluster[1][1];
|
||||
|
||||
corner = 1;
|
||||
sumMax=sumBR;
|
||||
}
|
||||
|
||||
switch(corner){
|
||||
case 0:
|
||||
cES3vs2->Fill(sum,sumTR); break;
|
||||
case 1:
|
||||
cES3vs2->Fill(sum,sumBR); break;
|
||||
case 2:
|
||||
cES3vs2->Fill(sum,sumTL); break;
|
||||
case 3:
|
||||
cES3vs2->Fill(sum,sumBL); break;
|
||||
}
|
||||
|
||||
counter[corner]++;
|
||||
remoteCounter[q]++;
|
||||
|
||||
// cout << "local corner is: " << corner << " remote corner is: " << q << endl;
|
||||
|
||||
return corner;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
int placePhoton( TH2D *img, double subCluster[2][2], int cX, int cY, int corner, double *sX, double *sY, double *scX, double *scY){
|
||||
double tot = subCluster[0][0] + subCluster[0][1] + subCluster[1][0] + subCluster[1][1];
|
||||
double t = subCluster[1][0] + subCluster[1][1];
|
||||
double r = subCluster[0][1] + subCluster[1][1];
|
||||
|
||||
double xHitC = r/tot;
|
||||
double yHitC = t/tot;
|
||||
|
||||
imgRLR->Fill(cX,cY);
|
||||
|
||||
cE->Fill(tot);
|
||||
|
||||
double dX, dY;
|
||||
|
||||
//before looking at annas code
|
||||
/* if(corner == 0){ dX=-1.; dY=-1.; }
|
||||
if(corner == 1){ dX=-1.; dY=+1.; }
|
||||
if(corner == 2){ dX=+1.; dY=-1.; }
|
||||
if(corner == 3){ dX=+1.; dY=+1.; }*/
|
||||
|
||||
if(corner == 0){ dX=-1.; dY=+1.; } //top left
|
||||
if(corner == 1){ dX=+1.; dY=+1.; } //top right
|
||||
if(corner == 2){ dX=-1.; dY=-1.; } //bottom left
|
||||
if(corner == 3){ dX=+1.; dY=-1.; } //bottom right
|
||||
|
||||
imgLR->Fill(cX+0.25*dX,cY+0.25*dY);
|
||||
|
||||
double posX = ((double)cX) + 0.5*dX + xHitC;
|
||||
double posY = ((double)cY) + 0.5*dY + yHitC;
|
||||
|
||||
subPosBEta->Fill(xHitC ,yHitC);
|
||||
if(img){
|
||||
img->Fill(posX,posY);
|
||||
}
|
||||
|
||||
if(xHitC < 0.02&& yHitC < 0.02){
|
||||
|
||||
cES3vs2S->Fill(dum[0][0]+dum[0][1]+dum[0][2]+dum[1][0]+dum[1][1]+dum[1][2]+dum[2][0]+dum[2][1]+dum[2][2],subCluster[0][0]+subCluster[0][1]+subCluster[1][0]+subCluster[1][1]);
|
||||
}
|
||||
|
||||
|
||||
if(sX && sY && scX && scY){
|
||||
*sX = xHitC; //0.5 + 0.5*dX + xHitC;
|
||||
*sY = yHitC; //0.5 + 0.5*dY + yHitC;
|
||||
*scX = ((double)cX) + 0.5*dX;
|
||||
*scY = ((double)cY) + 0.5*dY;
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void placePhotonCorr(TH2D *img, EtaVEL *e,double sX, double sY, double scX, double scY){
|
||||
int bin = e->findBin(sX,sY);
|
||||
if(bin <= 0) return;
|
||||
double subX = ((double)(e->getXBin(bin))+.5)/((double)e->getNPixels());
|
||||
double subY = ((double)(e->getYBin(bin))+.5)/((double)e->getNPixels());
|
||||
|
||||
if(img!=NULL){
|
||||
img->Fill(scX+ subX , scY+ subY);
|
||||
}
|
||||
subPosAEta->Fill(subX,subY);
|
||||
|
||||
int iscx = scX;
|
||||
int iscy = scY;
|
||||
if(iscx >=nx || iscx<0 || iscy >=ny || iscy<0) return;
|
||||
//cout << iscx*e->getNPixels()+e->getXBin(bin) << " " << iscy*e->getNPixels()+e->getXBin(bin) << endl;
|
||||
if(img==NULL) return;
|
||||
imgArray[iscx*e->getNPixels()+e->getXBin(bin)][iscy*e->getNPixels()+e->getYBin(bin)]++;
|
||||
}
|
||||
|
||||
void gainCorrection(Double_t corrected[3][3], TH2D *gainMap){
|
||||
|
||||
for(int xx = 0; xx < 3; xx++)
|
||||
for(int yy = 0; yy < 3; yy++){
|
||||
if(gainMap && gainMap->GetBinContent(x+xx+2,y+yy+2) != 0){
|
||||
corrected[xx][yy] = dum[xx][yy] / gainMap->GetBinContent(x+xx+2,y+yy+2);
|
||||
clusHistC->Fill(xx,yy,corrected[xx][yy]);
|
||||
}
|
||||
else
|
||||
corrected[xx][yy] = dum[xx][yy];
|
||||
|
||||
clusHist->Fill(xx,yy,dum[xx][yy]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
EtaVEL *plotEtaDensity(TChain* tree2, TEntryList *el, EtaVEL *oldEta = NULL, TH2D **img = NULL, TH2D *gainMap=NULL, int nPixels=25) {
|
||||
|
||||
|
||||
|
||||
EtaVEL *newEta = new EtaVEL(25,-0.02,1.02);
|
||||
|
||||
Long64_t listEntries=el->GetN();
|
||||
Long64_t treeEntry;
|
||||
Long64_t chainEntry;
|
||||
|
||||
Int_t treenum=0;
|
||||
tree2->SetEntryList(el);
|
||||
|
||||
double gainCorrC[3][3];
|
||||
double subCluster[2][2];
|
||||
double sX, sY, scX, scY;
|
||||
|
||||
cout << "Events: " << listEntries << endl;
|
||||
if(oldEta == NULL){ cout << "Old Eta is NULL " << endl; }
|
||||
for(int i = 0; i<4; i++){ counter[i] = 0; remoteCounter[i] = 0; }
|
||||
|
||||
for (Long64_t il =0; il<listEntries;il++) {
|
||||
treeEntry = el->GetEntryAndTree(il,treenum);
|
||||
|
||||
chainEntry = treeEntry+tree2->GetTreeOffset()[treenum];
|
||||
if (tree2->GetEntry(chainEntry)) {
|
||||
|
||||
gainCorrection(gainCorrC,gainMap);
|
||||
//cout << gainCorrC[1][1] << endl;
|
||||
|
||||
//finds corner
|
||||
int corner = findShape(gainCorrC,subCluster);
|
||||
|
||||
int validEvent;
|
||||
|
||||
|
||||
if(img){
|
||||
validEvent = placePhoton(img[0],subCluster,x,y, corner, &sX, &sY, &scX, &scY);
|
||||
}else{
|
||||
//calc etaX, etaY
|
||||
validEvent = placePhoton(NULL,subCluster,x,y, corner, &sX, &sY, &scX, &scY);
|
||||
}
|
||||
|
||||
//fill etavel
|
||||
newEta->fill(sX,sY);
|
||||
|
||||
|
||||
|
||||
|
||||
if(oldEta && img && img[1]){
|
||||
placePhotonCorr(img[1],oldEta, sX,sY, scX, scY);
|
||||
}else{
|
||||
placePhotonCorr(NULL,newEta,sX,sY,scX,scY);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
//cout << il << endl;
|
||||
int ssize = 500000;
|
||||
if(il % ssize == 0 && il != 0 && oldEta==NULL){
|
||||
|
||||
cout << " -------------- "<< endl;
|
||||
newEta->updatePixelPos();
|
||||
|
||||
|
||||
//newEta->resolveSelfIntersect();
|
||||
char tit[1000];
|
||||
/* TFile *ff = new TFile("/scratch/Spider.root","UPDATE");
|
||||
sprintf(tit,"subPosAEta%i",newEta->getIt()); subPosAEta->SetName(tit);
|
||||
subPosAEta->Write(); subPosAEta->Reset();
|
||||
sprintf(tit,"subPosBEta%i",newEta->getIt()); subPosBEta->SetName(tit);
|
||||
subPosBEta->Write(); subPosBEta->Reset();
|
||||
sprintf(tit,"Eta%i",newEta->getIt()); newEta->Write(tit);
|
||||
ff->Close(); */
|
||||
//il = 0;
|
||||
}
|
||||
|
||||
if(il % ssize == ssize-1){
|
||||
double prog = (double)il/(double)listEntries*100.;
|
||||
cout << prog << "%" << endl;
|
||||
//if(prog > 19.) return newEta;
|
||||
if(newEta->converged == 1){ cout << "converged ... " << endl; return newEta; }
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
cout << "local corners: " ;
|
||||
for(int i = 0; i<4; i++) cout << i << ": " << counter[i] << " || " ;
|
||||
cout << endl;
|
||||
|
||||
//cout << "remote corners: " ;
|
||||
//for(int i = 0; i<4; i++) cout << i << ": " << remoteCounter[i] << " || " ;
|
||||
//cout << endl;
|
||||
|
||||
return newEta;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
TChain *openTree(char *tname, char *fname,double lEc, double hEc, double rms=5., char *chainName=">>thischan"){
|
||||
TChain *tree2;
|
||||
// TH1D **etaDI;
|
||||
char cut[1000];
|
||||
|
||||
tree2=new TChain(tname);
|
||||
tree2->Add(fname);
|
||||
tree2->Print();
|
||||
|
||||
//sprintf(cut,"(x<=40) && (data[%d][%d]>%f*rms) && Sum$(data)<%f && Sum$(data)>%f",1,1,rms, hEc, lEc);
|
||||
// sprintf(cut,"(x<=40) && (data[%d][%d]>%f*rms)",1,1,rms);// && Sum$(data)<%f && Sum$(data)>%f",1,1,rms, hEc, lEc);
|
||||
sprintf(cut,"(x<=40) && Sum$(data)<%f && Sum$(data)>%f", hEc, lEc);
|
||||
// sprintf(cut,"");
|
||||
cout << cut << endl;
|
||||
|
||||
tree2->Draw(chainName, cut, "entrylist");
|
||||
|
||||
|
||||
tree2->SetBranchAddress("iFrame",&f);
|
||||
tree2->SetBranchAddress("x",&x);
|
||||
tree2->SetBranchAddress("y",&y);
|
||||
tree2->SetBranchAddress("data",dum);
|
||||
//tree2->SetBranchAddress("q",&q);
|
||||
|
||||
cout << "openTree : end" << endl;
|
||||
return tree2;
|
||||
}
|
||||
|
||||
EtaVEL *etaDensity(char *tname, char *fname, double lEc = 1000, double hEc=3000, TH2D *gainMap=NULL, int nPixels=25) {
|
||||
/** open tree and make selection */
|
||||
TChain *tree2 = openTree(tname,fname,lEc,hEc);
|
||||
TEntryList *elist = (TEntryList*)gDirectory->Get("thischan");
|
||||
if(elist == NULL) { cout << "could not open tree " << endl; return NULL; }
|
||||
|
||||
EtaVEL *etaDen = plotEtaDensity(tree2,elist,NULL,NULL,gainMap,nPixels);
|
||||
|
||||
|
||||
//etaDen->Draw("colz");
|
||||
cout << "done" << endl;
|
||||
|
||||
return etaDen;
|
||||
}
|
||||
|
||||
void interpolate(char *tname, char *fname, EtaVEL *etaDI, double lEc = 1000, double hEc=3000, TH2D *gainMap=NULL) {
|
||||
|
||||
TChain *tree2 = openTree(tname,fname,lEc,hEc,5.,">>intChain");
|
||||
TEntryList *elist = (TEntryList*)gDirectory->Get("intChain");
|
||||
if(elist == NULL) { cout << "could not open tree " << endl; return; }
|
||||
|
||||
double nPixels = (double)etaDI->getNPixels();
|
||||
|
||||
TH2D **img = new TH2D*[3];
|
||||
img[0] = new TH2D("img","img",nPixels*160,0.0,160.0 ,nPixels*160 ,0.0,160.0);
|
||||
img[1] = new TH2D("imgE","imgE",nPixels*160,0.0,160.0 ,nPixels*160 ,0.0,160.0);
|
||||
|
||||
int inPixels = etaDI->getNPixels();
|
||||
|
||||
imgArray = new int*[inPixels*160];
|
||||
for(int i = 0; i < inPixels*160; i++){
|
||||
imgArray[i] = new int[inPixels*160];
|
||||
for(int j = 0; j < inPixels*160; j++){
|
||||
imgArray[i][j] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
cout << "starting" << endl;
|
||||
plotEtaDensity(tree2,elist, etaDI,img,gainMap);
|
||||
|
||||
//img->Draw("colz");
|
||||
}
|
||||
|
||||
|
||||
TH2D *createGainMap(char *tname, char *fname, double lEc = 0,double hEc=10000){
|
||||
char name[100];
|
||||
TH1D *avgSpec3 = new TH1D("avgSpec3", "avgSpec3",hEc/20,0,hEc);
|
||||
TH1D ***specs3 = new TH1D**[160];
|
||||
TH1D ***specs1 = new TH1D**[160];
|
||||
for(int xx = 0; xx < 160; xx++){
|
||||
specs3[xx] = new TH1D*[160];
|
||||
specs1[xx] = new TH1D*[160];
|
||||
for(int yy = 0; yy < 160; yy++){
|
||||
sprintf(name,"S3x%iy%i",xx,yy);
|
||||
specs3[xx][yy] = new TH1D(name,name,hEc/20,0,hEc);
|
||||
sprintf(name,"S1x%iy%i",xx,yy);
|
||||
specs1[xx][yy] = new TH1D(name,name,hEc/20,0,hEc);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
TChain *tree2 = openTree(tname,fname,0,hEc,5.,">>gainChan");
|
||||
TEntryList *elist = (TEntryList*)gDirectory->Get("gainChan");
|
||||
if(elist == NULL) { cout << "could not open tree " << endl; return NULL; }
|
||||
|
||||
Long64_t listEntries=elist->GetN();
|
||||
Long64_t treeEntry;
|
||||
Long64_t chainEntry;
|
||||
|
||||
Int_t treenum=0;
|
||||
tree2->SetEntryList(elist);
|
||||
|
||||
cout << "Events: " << listEntries << endl;
|
||||
for(int i = 0; i<4; i++) counter[i] = 0;
|
||||
for (Long64_t il =0; il<listEntries;il++) {
|
||||
treeEntry = elist->GetEntryAndTree(il,treenum);
|
||||
chainEntry = treeEntry+tree2->GetTreeOffset()[treenum];
|
||||
|
||||
if (tree2->GetEntry(chainEntry)) {
|
||||
double sum = 0;
|
||||
for(int xx = 0; xx < 3; xx++)
|
||||
for(int yy = 0; yy < 3; yy++)
|
||||
sum += dum[xx][yy];
|
||||
specs3[x][y]->Fill(sum);
|
||||
specs1[x][y]->Fill(dum[1][1]);
|
||||
avgSpec3->Fill(sum);
|
||||
}
|
||||
}
|
||||
|
||||
TH2D *gainMap3 = new TH2D("gainMap3","gainMap3",160,-0.5,160.-0.5,160,-.5,160.-.5);
|
||||
TH2D *gainMap1 = new TH2D("gainMap1","gainMap1",160,-0.5,160.-0.5,160,-.5,160.-.5);
|
||||
for(int xx = 0; xx < 160; xx++){
|
||||
for(int yy = 0; yy < 160; yy++){
|
||||
TF1 *gf3 = new TF1("gf3","gaus", lEc, hEc);
|
||||
specs3[xx][yy]->Fit(gf3,"Q");
|
||||
double e3 = gf3->GetParameter(1);
|
||||
gainMap3->Fill(xx,yy,e3);
|
||||
|
||||
TF1 *gf1 = new TF1("gf1","gaus", lEc, hEc);
|
||||
specs1[xx][yy]->Fit(gf1,"Q");
|
||||
double e1 = gf1->GetParameter(1);
|
||||
gainMap1->Fill(xx,yy,e1);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
return gainMap3;
|
||||
}
|
||||
|
||||
void writeMatlab2DHisto(int xx, int yy,char *outFileName){
|
||||
ofstream outFile;
|
||||
outFile.open (outFileName);
|
||||
|
||||
cout << "create matlab file with " << xx << " xbins and " << yy << " ybins" << endl;
|
||||
|
||||
for(int y = 0; y < yy; y++){
|
||||
for(int x = 0; x < xx; x++){
|
||||
outFile << imgArray[x][y] << "\t";
|
||||
}
|
||||
outFile << endl;
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
}
|
||||
|
||||
//COMPLETE STUFF
|
||||
|
||||
void createImage(char *tdir, char *tname, char *ftname, char *ifname = NULL, int useGM=0, double lEth=-1., double hEth=-1.){
|
||||
imgRLR->Reset();
|
||||
imgLR->Reset();
|
||||
|
||||
char fname[1000];
|
||||
char inFName[1000];
|
||||
char outFName[1000];
|
||||
char moutFName[1000];
|
||||
if(ifname == NULL){
|
||||
sprintf(fname,"%s/%s_*.root",tdir,tname);
|
||||
}else{
|
||||
sprintf(fname,"%s",ifname);
|
||||
}
|
||||
|
||||
if(useGM) sprintf(inFName,"%s/%s-PlotsWGMVEL.root",tdir,ftname);
|
||||
else sprintf(inFName,"%s/%s-PlotsVEL.root",tdir,ftname);
|
||||
|
||||
sprintf(outFName,"%s/%s-ImgVEL.root",tdir,tname);
|
||||
sprintf(moutFName,"%s/%s-ImgVEL.mf",tdir,tname);
|
||||
|
||||
TFile *inFile = new TFile(inFName,"READ");
|
||||
|
||||
cout << "Image Tree File Name: " << fname << endl;
|
||||
cout << "Eta File Name: " << inFName << endl;
|
||||
cout << "Out File Name: " << outFName << endl;
|
||||
cout << "Matlab Out File Name: " << moutFName << endl;
|
||||
|
||||
TH2D *gm = NULL;
|
||||
if(useGM){
|
||||
cout << "Load gain map" << endl;
|
||||
gm = (TH2D *)gDirectory->Get("gainMap");
|
||||
if(gm == NULL){ cout << "can not find gainMap in file" << endl; return; }
|
||||
}
|
||||
|
||||
cout << "Load eta" << endl;
|
||||
EtaVEL *ee = (EtaVEL *)gDirectory->Get("etaDist");
|
||||
|
||||
cout << "Select Energy BW" << endl;
|
||||
TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
|
||||
if(spec == NULL){ cout << "can not find avgSpec3" << endl; return; }
|
||||
|
||||
TF1 *gf3 = new TF1("gf3","gaus", 0, 10000);
|
||||
spec->Fit(gf3,"Q");
|
||||
double avgE = gf3->GetParameter(1);
|
||||
double sigE = gf3->GetParameter(2);
|
||||
cout << "avgE: " << avgE << " sigE: " << sigE << endl;
|
||||
cout << endl;
|
||||
|
||||
if(lEth == -1.) lEth = avgE-5.*sigE;
|
||||
if(hEth == -1.) hEth = avgE+5.*sigE;
|
||||
cout << lEth << " < E < " << hEth << " (eV)" << endl;
|
||||
|
||||
cout << "start with interpolation" << endl;
|
||||
interpolate( tname, fname, ee,lEth,hEth ,gm);
|
||||
|
||||
|
||||
TH2D *img = (TH2D *)gDirectory->Get("img");
|
||||
if(img == NULL){ cout << "could not find 2d-histogram: img " << endl; return; }
|
||||
|
||||
|
||||
TH2D *imgE = (TH2D *)gDirectory->Get("imgE");
|
||||
if(imgE == NULL){ cout << "could not find 2d-histogram: imgE " << endl; return; }
|
||||
|
||||
|
||||
//TH2D *imgEOM = (TH2D *)gDirectory->Get("imgEOM");
|
||||
//if(imgEOM == NULL){ cout << "could not find 2d-histogram: imgEOM " << endl; return; }
|
||||
|
||||
TFile *outFile = new TFile(outFName,"UPDATE");
|
||||
imgLR->Write();
|
||||
imgRLR->Write();
|
||||
imgE->Write();
|
||||
//imgEOM->Write();
|
||||
img->Write();
|
||||
outFile->Close();
|
||||
inFile->Close();
|
||||
cout << "writing matlab file: " << moutFName << endl;
|
||||
writeMatlab2DHisto(160*ee->getNPixels(),160*ee->getNPixels(),moutFName);
|
||||
cout << "Done : " << outFName << endl;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
\par tdir input tree directory
|
||||
\par tname input tree name
|
||||
\par ifname input file name if different than tdir/tname_*.root
|
||||
\par useGM use gain map
|
||||
\par maxExpEinEv spectrum maximum
|
||||
\par nPixels sub-pixels bins
|
||||
\par lEth low threshold
|
||||
\par hEth high threshold
|
||||
|
||||
*/
|
||||
|
||||
|
||||
EtaVEL *createGainAndEtaFile(char *tdir, char *tname, char *ifname=NULL, int useGM=0, double maxExpEinEv=25000., int nPixels =25, double lEth=-1., double hEth=-1.){
|
||||
char fname[1000];
|
||||
char outFName[1000];
|
||||
|
||||
|
||||
if(ifname == NULL){
|
||||
sprintf(fname,"%s/%s_*.root",tdir,tname);
|
||||
}else{
|
||||
sprintf(fname,"%s",ifname);
|
||||
}
|
||||
|
||||
if(useGM) sprintf(outFName,"%s/%s-PlotsWGVEL.root",tdir,tname);
|
||||
else sprintf(outFName,"%s/%s-PlotsVEL.root",tdir,tname);
|
||||
|
||||
|
||||
cout << "Tree File Name: " << fname << endl;
|
||||
cout << "Output File Name: " << outFName << endl;
|
||||
|
||||
/** creates gain map and 3x3 spectrum */
|
||||
cout << "Creating gain map: " << endl;
|
||||
TH2D *gm = createGainMap(tname,fname,0,maxExpEinEv/10.);
|
||||
gm->SetName("gainMap");
|
||||
|
||||
|
||||
/** gets average 3x3 spectrum and fits it with a gaus */
|
||||
TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
|
||||
if(spec == NULL){ cout << "can not find avgSpec3" << endl; return NULL; }
|
||||
TF1 *gf3 = new TF1("gf3","gaus", 0, maxExpEinEv/10.);
|
||||
spec->Fit(gf3,"Q");
|
||||
double avgE = gf3->GetParameter(1);
|
||||
double sigE = gf3->GetParameter(2);
|
||||
cout << "avgE: " << avgE << " sigE: " << sigE << endl;
|
||||
cout << endl;
|
||||
|
||||
|
||||
/** sets high and low threshold if not given by the user */
|
||||
if(lEth == -1.) lEth = avgE-5.*sigE;
|
||||
if(hEth == -1.) hEth = avgE+5.*sigE;
|
||||
cout << lEth << " < E < " << hEth << " (eV)" << endl;
|
||||
|
||||
|
||||
|
||||
|
||||
cout << "calculating eta stuff" << endl;
|
||||
|
||||
EtaVEL *newEta;
|
||||
if(useGM) newEta = etaDensity(tname,fname,lEth,hEth,gm,nPixels);
|
||||
else newEta = etaDensity(tname,fname,lEth,hEth,NULL,nPixels);
|
||||
|
||||
cout << "writing to file " << outFName << endl;
|
||||
|
||||
TFile *outFile = new TFile(outFName,"UPDATE");
|
||||
|
||||
newEta->Write("etaDist");
|
||||
|
||||
gm->Write();
|
||||
spec->Write();
|
||||
subPosAEta->Write();
|
||||
cES3vs2->Write();
|
||||
|
||||
outFile->Close();
|
||||
cout << "Done : " << outFName << endl;
|
||||
return newEta;
|
||||
}
|
||||
|
||||
void exportSpec(char *tdir, char *tname){
|
||||
char tfname[1000];
|
||||
char ofname[1000];
|
||||
char cleanName[1000];
|
||||
|
||||
for(int p = 0; p < strlen(tname);p++){
|
||||
cleanName[p+1] = '\0';
|
||||
cleanName[p] = tname[p];
|
||||
|
||||
if(tname[p] == '-') cleanName[p] = '_';
|
||||
}
|
||||
|
||||
sprintf(tfname,"%s/%s-PlotsVEL.root",tdir,tname);
|
||||
sprintf(ofname,"%s/%s_SpecVEL.m",tdir,cleanName);
|
||||
TFile *tf = new TFile(tfname);
|
||||
TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
|
||||
|
||||
ofstream outFile;
|
||||
outFile.open (ofname);
|
||||
|
||||
if(outFile.fail()){
|
||||
cout << "Could not open file : " << ofname << endl;
|
||||
return;
|
||||
}
|
||||
|
||||
cout << "create matlab file with with spec " << ofname << endl;
|
||||
|
||||
|
||||
outFile << cleanName << " = [ " << endl;
|
||||
for(int i = 0; i < spec->GetNbinsX(); i++){
|
||||
outFile << i << " " << spec->GetBinCenter(i) << " " << spec->GetBinContent(i) << " ; " << endl;
|
||||
}
|
||||
|
||||
outFile << " ] ; " << endl;
|
||||
|
||||
outFile.close();
|
||||
}
|
||||
679
slsDetectorCalibration/etaVEL/EtaVEL.cpp
Normal file
679
slsDetectorCalibration/etaVEL/EtaVEL.cpp
Normal file
@ -0,0 +1,679 @@
|
||||
#include "EtaVEL.h"
|
||||
#include <iomanip>
|
||||
|
||||
|
||||
ClassImp(EtaVEL);
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
|
||||
double *EtaVEL::getPixelCorners(int x, int y){
|
||||
double tlX,tlY,trX,trY,blX,blY,brX,brY;
|
||||
tlX = xPPos[getCorner(x,y+1)];
|
||||
tlY = yPPos[getCorner(x,y+1)];
|
||||
trX = xPPos[getCorner(x+1,y+1)];
|
||||
trY = yPPos[getCorner(x+1,y+1)];
|
||||
blX = xPPos[getCorner(x,y)];
|
||||
blY = yPPos[getCorner(x,y)];
|
||||
brX = xPPos[getCorner(x+1,y)];
|
||||
brY = yPPos[getCorner(x+1,y)];
|
||||
|
||||
//cout << "gPC: TL: " << getCorner(x,y+1) << " TR: " << getCorner(x+1,y+1) << " BL " << getCorner(x,y) << " BR " << getCorner(x+1,y) << endl;
|
||||
|
||||
double *c = new double[8];
|
||||
c[0] = tlX; c[1] = trX; c[2] = brX; c[3] = blX;
|
||||
c[4] = tlY; c[5] = trY; c[6] = brY; c[7] = blY;
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
int EtaVEL::findBin(double xx, double yy){
|
||||
|
||||
double tlX,tlY,trX,trY,blX,blY,brX,brY;
|
||||
/********Added by anna ******/
|
||||
// if (xx<min) xx=min+1E-6;
|
||||
// if (xx>max) xx=max-1E-6;
|
||||
// if (yy<min) yy=min+1E-6;
|
||||
// if (yy>max) yy=max-1E-6;
|
||||
/**************/
|
||||
|
||||
|
||||
int bin = -1;
|
||||
for(int x = 0; x < nPixels; x++){
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
double *c = getPixelCorners(x,y);
|
||||
tlX = c[0]; trX = c[1]; brX = c[2]; blX = c[3];
|
||||
tlY = c[4]; trY = c[5]; brY = c[6]; blY = c[7];
|
||||
|
||||
///if(y == 0){
|
||||
// cout << "x: " << x << " blY " << blY << " brY " << brY << endl;
|
||||
//}
|
||||
|
||||
int out = 0;
|
||||
|
||||
double tb = 0;
|
||||
double bb = 0;
|
||||
double lb = 0;
|
||||
double rb = 0;
|
||||
|
||||
if((trX-tlX)>0.)
|
||||
tb = (trY - tlY)/(trX-tlX);
|
||||
|
||||
if((brX-blX)>0.)
|
||||
bb = (brY - blY)/(brX-blX);
|
||||
|
||||
if((tlY-blY)>0.)
|
||||
lb = (tlX - blX)/(tlY-blY);
|
||||
|
||||
if((trY-brY)>0.)
|
||||
rb = (trX - brX)/(trY-brY);
|
||||
|
||||
double ty = tlY + tb * (xx - tlX);
|
||||
double by = blY + bb * (xx - blX);
|
||||
|
||||
double lx = blX + lb * (yy - blY);
|
||||
double rx = brX + rb * (yy - brY);
|
||||
|
||||
|
||||
|
||||
|
||||
if(yy >= ty) out++;
|
||||
if(yy < by) out++;
|
||||
if(xx < lx) out++;
|
||||
if(xx >= rx) out++;
|
||||
|
||||
//cout << "ty " << ty << endl;
|
||||
//cout << "by " << by << endl;
|
||||
//cout << "lx " << lx << endl;
|
||||
//cout << "rx " << rx << endl;
|
||||
|
||||
//double dist = (xx - xPPos[getBin(x,y)]) * (xx - xPPos[getBin(x,y)]) + (yy - yPPos[getBin(x,y)]) * (yy - yPPos[getBin(x,y)]);
|
||||
//cout << "x " << x << " y " << y << " out " << out << " ty " << ty << endl;
|
||||
//cout << "tl " << tlX << "/" << tlY << " tr " << trX << "/" << trY << endl;
|
||||
//cout << "bl " << blX << "/" << blY << " br " << brX << "/" << brY << endl;
|
||||
|
||||
//cout << " tb " << tb << endl;
|
||||
|
||||
|
||||
delete[] c;
|
||||
if(out == 0){ return getBin(x,y); }
|
||||
}
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
void EtaVEL::createLogEntry(){
|
||||
if(it >= nIterations){
|
||||
cerr << "log full" << endl;
|
||||
}
|
||||
log[it].itN = it;
|
||||
log[it].xPos = new double[nPixels*nPixels+1];
|
||||
log[it].yPos = new double[nPixels*nPixels+1];
|
||||
log[it].binCont = new double[nPixels*nPixels+1];
|
||||
for(int x = 0; x < nPixels; x++)
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
log[it].xPos[getBin(x,y)] = xPPos[getBin(x,y)];
|
||||
log[it].yPos[getBin(x,y)] = yPPos[getBin(x,y)];
|
||||
log[it].binCont[getBin(x,y)] = binCont[getBin(x,y)];
|
||||
}
|
||||
it++;
|
||||
}
|
||||
|
||||
void EtaVEL::updatePixelCorner(){
|
||||
double w = 20;
|
||||
int rows = (nPixels+1)*(nPixels+1) + 4 + 4 * 4;//(4*(nPixels+1))-4;
|
||||
int cols = (nPixels+1)*(nPixels+1);
|
||||
|
||||
double *rVx = new double[rows];
|
||||
double *rVy = new double[rows];
|
||||
|
||||
double *posMat = new double[rows*cols];
|
||||
for(int i = 0 ; i < rows*cols; i++) posMat[i] = 0;
|
||||
int boundaryPoint = 0;
|
||||
|
||||
cout << "linear sys stuff" << endl;
|
||||
|
||||
double minELength = 100000000000000; int minX=-1, minY=-1;
|
||||
|
||||
for(int y = 0; y < nPixels+1; y++){
|
||||
for(int x = 0; x < nPixels+1; x++){
|
||||
double bx = 0, by = 0;
|
||||
|
||||
//boundary conditions
|
||||
|
||||
if((x == 0 && y % 5 == 0) ||
|
||||
(x == nPixels && y % 5 == 0) ||
|
||||
(y == 0 && x % 5 == 0) ||
|
||||
(y == nPixels && x % 5 == 0)){
|
||||
|
||||
bx = xPPos[getCorner(x,y)];
|
||||
//cout << "bP " << boundaryPoint << " bx " << bx << endl;
|
||||
by = yPPos[getCorner(x,y)];
|
||||
rVx[(nPixels+1)*(nPixels+1) + boundaryPoint] = bx*w;
|
||||
rVy[(nPixels+1)*(nPixels+1) + boundaryPoint] = by*w;
|
||||
posMat[(nPixels+1)*(nPixels+1)*cols + boundaryPoint * cols + getCorner(x,y)-1] = w;
|
||||
boundaryPoint++;
|
||||
}
|
||||
|
||||
double tot = 4 - (x == 0) - (y == 0) - (x == nPixels) - (y == nPixels);
|
||||
//cout << "totW: " << tot << endl;
|
||||
//tot = 4.;
|
||||
double eLength = 0;
|
||||
if(x != 0) eLength += edgeL[getEdgeX(x-1,y)];
|
||||
if(y != 0) eLength += edgeL[getEdgeY(x,y-1)];
|
||||
if(x != nPixels) eLength += edgeL[getEdgeX(x,y)];
|
||||
if(y != nPixels) eLength += edgeL[getEdgeY(x,y)];
|
||||
|
||||
/*cout << "Corner X:" <<x << " Y: " << y ;
|
||||
cout << " C# " << getCorner(x,y);
|
||||
cout << " eXl " << getEdgeX(x-1,y) << "(C# " << getCorner(x-1,y) << " ) ";
|
||||
cout << " eXr " << getEdgeX(x,y) << "(C# " << getCorner(x+1,y) << " ) ";
|
||||
cout << " eYb " << getEdgeY(x,y-1) << "(C# " << getCorner(x,y-1) << " ) ";
|
||||
cout << " eYt " << getEdgeY(x,y) << "(C# " << getCorner(x,y+1) << " ) " << endl; */
|
||||
//" totW: " << tot << " totE: " << eLength << endl;
|
||||
|
||||
if(eLength < minELength & tot == 4){
|
||||
minELength = eLength;
|
||||
minX = x; minY = y;
|
||||
}
|
||||
|
||||
|
||||
//matrixes updated
|
||||
if(x != 0) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x-1,y)-1] = -edgeL[getEdgeX(x-1,y)]/eLength;
|
||||
if(y != 0) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x,y-1)-1] = -edgeL[getEdgeY(x,y-1)]/eLength;;
|
||||
if(x != nPixels) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x+1,y)-1] = -edgeL[getEdgeX(x,y)]/eLength;;
|
||||
if(y != nPixels) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x,y+1)-1] = -edgeL[getEdgeY(x,y)]/eLength;;
|
||||
|
||||
posMat[y*(nPixels+1)*cols+x*cols+getCorner(x,y)-1] = 1.;
|
||||
rVx[getCorner(x,y)-1] = 0.;
|
||||
rVy[getCorner(x,y)-1] = 0.;
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
cout << "Min Corner X: " <<minX << " Y: " << minY << " C# " << getCorner(minX,minY) << " length " << minELength << endl;
|
||||
|
||||
TMatrixD *k = new TMatrixD(rows,cols);
|
||||
TVectorD *fx = new TVectorD(rows,rVx);
|
||||
TVectorD *fy = new TVectorD(rows,rVy);
|
||||
// f->Print();
|
||||
k->SetMatrixArray(posMat);
|
||||
// k->Print();
|
||||
|
||||
|
||||
//solve linear system
|
||||
|
||||
Bool_t ok;
|
||||
TDecompSVD *s = new TDecompSVD(*k);
|
||||
s->Solve(*fx);
|
||||
s->Solve(*fy);
|
||||
|
||||
double *fxA = fx->GetMatrixArray();
|
||||
double *fyA = fy->GetMatrixArray();
|
||||
|
||||
|
||||
for(int y = 0; y < nPixels+1; y++){
|
||||
for(int x = 0; x < nPixels+1; x++){
|
||||
//do not update boundaries
|
||||
|
||||
if(!(x == 0 ||
|
||||
x == nPixels||
|
||||
y == 0 ||
|
||||
y == nPixels)){
|
||||
xPPos[getCorner(x,y)] = fxA[getCorner(x,y)-1];
|
||||
yPPos[getCorner(x,y)] = fyA[getCorner(x,y)-1];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void EtaVEL::updatePixelPos(){
|
||||
double xMov, yMov, d1Mov, d2Mov;
|
||||
createLogEntry();
|
||||
double *chMap = getChangeMap();
|
||||
int ch =0;
|
||||
|
||||
cout << "update edge lengths" << endl;
|
||||
for(int x = 0; x < nPixels; x++)
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
|
||||
|
||||
/*cout << "Pixel X:" <<x << " Y: " << y << " P# " << getBin(x,y) << " eXb " << getEdgeX(x,y);
|
||||
cout << " eXt " << getEdgeX(x,y+1) << " eYl " << getEdgeY(x,y) << " eYr " << getEdgeY(x+1,y) << endl;
|
||||
*/
|
||||
|
||||
edgeL[getEdgeX(x,y)] *= chMap[getBin(x,y)];
|
||||
edgeL[getEdgeX(x,y+1)] *= chMap[getBin(x,y)];
|
||||
edgeL[getEdgeY(x,y)] *= chMap[getBin(x,y)];
|
||||
edgeL[getEdgeY(x+1,y)] *= chMap[getBin(x,y)];
|
||||
|
||||
//cout << "Pixel x: " << x << " y: " << y << " Ch: " << chMap[getBin(x,y)] << " counts: " << binCont[getBin(x,y)] << endl;
|
||||
//cout << "BE " << getEdgeX(x,y) << endl;
|
||||
//cout << "TE " << getEdgeX(x,y+1) << endl;
|
||||
//cout << "LE " << getEdgeY(x,y) << endl;
|
||||
//cout << "RE " << getEdgeY(x+1,y) << endl;
|
||||
binCont[getBin(x,y)] = 0;
|
||||
}
|
||||
|
||||
updatePixelCorner();
|
||||
|
||||
//double *pSize = getSizeMap();
|
||||
double totEdgeLength = 0;
|
||||
for(int e = 1; e < 2*nPixels*(nPixels+1)+1; e++){
|
||||
totEdgeLength += edgeL[e];
|
||||
}
|
||||
cout << "tot edge Length: " << totEdgeLength << endl;
|
||||
|
||||
totCont = 0.;
|
||||
|
||||
}
|
||||
|
||||
double *EtaVEL::getSizeMap(){
|
||||
double tlX,tlY,trX,trY,blX,blY,brX,brY;
|
||||
double *szMap = new double[nPixels*nPixels+1];
|
||||
for(int x = 1; x < nPixels-1; x++)
|
||||
for(int y = 1; y < nPixels-1; y++){
|
||||
double *c = getPixelCorners(x,y);
|
||||
tlX = c[0]; trX = c[1]; brX = c[2]; blX = c[3];
|
||||
tlY = c[4]; trY = c[5]; brY = c[6]; blY = c[7];
|
||||
|
||||
//double area = dtl * dtr / 2. + dtr * dbr / 2. + dbr * dbl / 2. + dbl * dtl / 2.;
|
||||
|
||||
//http://en.wikipedia.org/wiki/Shoelace_formula
|
||||
double sl1 = tlX * trY + trX * brY + brX * blY + blX * tlY;
|
||||
double sl2 = tlY * trX + trY * brX + brY * blX + blY * tlX;
|
||||
double area = 1./2. * (- sl1 + sl2);
|
||||
if(area < 0.){
|
||||
cout << "negative area: X " << x << " Y " << y << " area " << endl;
|
||||
edgeL[getEdgeX(x,y)] *= 2.;
|
||||
edgeL[getEdgeX(x,y+1)] *= 2.;
|
||||
edgeL[getEdgeY(x,y)] *= 2.;
|
||||
edgeL[getEdgeY(x+1,y)] *= 2.;
|
||||
|
||||
}
|
||||
szMap[getBin(x,y)] = area / (max - min) / (max - min) * nPixels * nPixels;
|
||||
delete[] c;
|
||||
|
||||
}
|
||||
return szMap;
|
||||
}
|
||||
|
||||
double *EtaVEL::getChangeMap(){
|
||||
double *chMap = new double[nPixels*nPixels+1];
|
||||
double avg = totCont/(double)(nPixels*nPixels);
|
||||
// TH1D *hmed=getCounts();
|
||||
// double med = Median(hmed);
|
||||
// delete hmed;
|
||||
double acc = TMath::Sqrt(avg);
|
||||
cout << "totC: " << totCont << " avg " << avg << " acc: " << acc << endl;//<< " med " << med
|
||||
double totOffAcc = 0.;
|
||||
int totInRange03s = 0;
|
||||
int totInRange07s = 0;
|
||||
int totInRange12s = 0;
|
||||
int totInRange20s = 0;
|
||||
int totInRange25s = 0;
|
||||
double dd;
|
||||
int totInBins = 0;
|
||||
|
||||
//double
|
||||
chi_sq=0;
|
||||
|
||||
int maxC = 0, maxX=-1, maxY=-1;
|
||||
double minC = 1000000000000000, minX, minY;
|
||||
|
||||
for(int x = 0; x < nPixels; x++){
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
totInBins += binCont[getBin(x,y)];
|
||||
double r = (double)binCont[getBin(x,y)];
|
||||
if(r > 0. & totCont > 0.){
|
||||
dd=sqrt(r/avg);
|
||||
/**Added by Anna */
|
||||
if (dd>2.) dd=1.5;
|
||||
if (dd<0.5) dd=0.75;
|
||||
chMap[getBin(x,y)] = dd;
|
||||
/** */
|
||||
//if( chMap[getBin(x,y)] < 1.){ chMap[getBin(x,y)] = 1/1.2; }
|
||||
//if( chMap[getBin(x,y)] > 1.){ chMap[getBin(x,y)] = 1.2; }
|
||||
//if( chMap[getBin(x,y)] < 1/1.2){ chMap[getBin(x,y)] = 1/1.2; }
|
||||
//if( chMap[getBin(x,y)] > 1.2){ chMap[getBin(x,y)] = 1.2; }
|
||||
}else if(totCont > 0.){
|
||||
chMap[getBin(x,y)] =0.5; //1/1.2;
|
||||
}else{
|
||||
chMap[getBin(x,y)] = 1.;
|
||||
}
|
||||
|
||||
//if(r < avg + 2*acc && r > avg - 2*acc){ totInRange++;}// chMap[getBin(x,y)] = 1.; }
|
||||
|
||||
/** Commente away by Anna
|
||||
if(converged == 0 && r < med+20*acc){ chMap[getBin(x,y)] = 1.; }
|
||||
if(converged == 2 && r < med+20*acc && r > med-03*acc){ chMap[getBin(x,y)] = 1.; }
|
||||
if(r < med+03*acc){ totInRange03s++; }
|
||||
if(r < med+07*acc){ totInRange07s++; }
|
||||
if(r < med+12*acc){ totInRange12s++; }
|
||||
if(r < med+20*acc){ totInRange20s++; }
|
||||
if(r < med+25*acc){ totInRange25s++; }
|
||||
*/
|
||||
|
||||
//cout << "x " << x << " y " << y << " r " << r << " ch " << chMap[getBin(x,y)] << endl;
|
||||
// if(r - avg > acc){ totOffAcc += r-avg;}
|
||||
//if(r - avg < -acc){ totOffAcc += avg-r;}
|
||||
totOffAcc += (avg-r)*(avg-r);
|
||||
chi_sq+=(avg-r)*(avg-r)/r;
|
||||
//cout << " x " << x << " y " << y << " bC " << binCont[x*nPixels+y] << " r " << r << endl;
|
||||
|
||||
if(r > maxC){ maxC = r; maxX = x; maxY = y; }
|
||||
if(r < minC){minC = r; minX = x; minY = y; }
|
||||
|
||||
}
|
||||
}
|
||||
// cout << "totInBins " << totInBins << " zero Bin " << binCont[0] << endl;
|
||||
cout << "AvgOffAcc: " << sqrt(totOffAcc/(double)(nPixels*nPixels)) << endl;
|
||||
cout << "***********Reduced Chi Square: " << chi_sq/((double)(nPixels*nPixels)) << endl;
|
||||
// cout << "totInRange03 (<" << med+03*acc << "): " << totInRange03s << endl;
|
||||
// cout << "totInRange07 (<" << med+07*acc << "): " << totInRange07s << endl;
|
||||
// cout << "totInRange12 (<" << med+12*acc << "): " << totInRange12s << endl;
|
||||
// cout << "totInRange20 (<" << med+20*acc << "): " << totInRange20s << endl;
|
||||
// cout << "totInRange25 (<" << med+25*acc << "): " << totInRange25s << endl;
|
||||
double maxSig = (maxC - avg)*(maxC - avg) / avg;//acc;
|
||||
double minSig = (avg - minC)*(avg - minC) / avg;//acc;
|
||||
cout << "Max Pixel X: " << maxX << " Y: " << maxY << " P# " << getBin(maxX,maxY) << " count: " << maxC << " sig: "<< maxSig << endl;
|
||||
cout << "Min Pixel X: " << minX << " Y: " << minY << " P# " << getBin(minX,minY) << " count: " << minC << " sig: "<< minSig << endl;
|
||||
|
||||
// if(maxSig <= 25){ converged = 2; cout << "reached first converstion step!!!" << endl; }
|
||||
//if(minSig <= 7 && converged == 2) { converged = 1; }
|
||||
if (chi_sq<3) converged=2;
|
||||
if (chi_sq<1) converged=1;
|
||||
cout << "Conversion step "<< converged << endl;
|
||||
return chMap;
|
||||
}
|
||||
|
||||
TH2D *EtaVEL::getContent(int it, int changeType){
|
||||
TH2D *cont = new TH2D("cont","cont",nPixels,min,max,nPixels,min,max);
|
||||
double *chMap = NULL;
|
||||
if(changeType ==1) chMap = getChangeMap();
|
||||
double *szMap = getSizeMap();
|
||||
for(int x = 0; x < nPixels; x++)
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
if(changeType ==2 ){
|
||||
cont->SetBinContent(x+1,y+1,szMap[getBin(x,y)]);
|
||||
}
|
||||
if(changeType ==1 ){
|
||||
cont->SetBinContent(x+1,y+1,chMap[getBin(x,y)]);
|
||||
}
|
||||
if(changeType ==0 ){
|
||||
if(it == -1){
|
||||
cont->SetBinContent(x+1,y+1,binCont[getBin(x,y)]);
|
||||
//cout << "x " << x << " y " << y << " cont " << binCont[getBin(x,y)] << endl;
|
||||
}
|
||||
else{cont->SetBinContent(x+1,y+1,log[it].binCont[getBin(x,y)]);}
|
||||
}
|
||||
}
|
||||
return cont;
|
||||
}
|
||||
|
||||
TH1D *EtaVEL::getCounts(){
|
||||
TH1D *ch = new TH1D("ch","ch",500,0,totCont/(nPixels*nPixels)*4);
|
||||
for(int x = 0; x < nPixels; x++)
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
ch->Fill(binCont[getBin(x,y)]);
|
||||
}
|
||||
return ch;
|
||||
|
||||
}
|
||||
|
||||
void EtaVEL::printGrid(){
|
||||
|
||||
double *colSum = new double[nPixels+1];
|
||||
double *rowSum = new double[nPixels+1];
|
||||
|
||||
for(int i = 0; i < nPixels+1; i++){
|
||||
colSum[i] = 0.;
|
||||
rowSum[i] = 0.;
|
||||
for(int j = 0; j < nPixels; j++){
|
||||
rowSum[i] += edgeL[getEdgeX(j,i)];
|
||||
colSum[i] += edgeL[getEdgeY(i,j)];
|
||||
}
|
||||
}
|
||||
|
||||
cout << endl;
|
||||
|
||||
cout.precision(3); cout << fixed;
|
||||
cout << " ";
|
||||
for(int x = 0; x < nPixels+1; x++){
|
||||
cout << setw(2) << x << " (" << colSum[x] << ") ";
|
||||
}
|
||||
cout << endl;
|
||||
for(int y = 0; y < nPixels+1; y++){
|
||||
cout << setw(2) << y << " ";
|
||||
for(int x = 0; x < nPixels+1; x++){
|
||||
cout << "(" << xPPos[getCorner(x,y)] << "/" << yPPos[getCorner(x,y)] << ") " ;
|
||||
if(x < nPixels) cout << " -- " << edgeL[getEdgeX(x,y)]/rowSum[y]*(max-min) << " -- ";
|
||||
}
|
||||
cout << " | " << rowSum[y] << endl;
|
||||
|
||||
if(y < nPixels){
|
||||
cout << " ";
|
||||
for(int x = 0; x < nPixels+1; x++){
|
||||
cout << edgeL[getEdgeY(x,y)]/colSum[x]*(max-min) << " ";
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
}
|
||||
delete[] colSum;
|
||||
delete[] rowSum;
|
||||
|
||||
}
|
||||
|
||||
TMultiGraph *EtaVEL::plotPixelBorder(int plotCenters){
|
||||
TMultiGraph *mg = new TMultiGraph();
|
||||
double cx[5], cy[5];
|
||||
for(int x = 0; x < nPixels; x++)
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
double *c = getPixelCorners(x,y);
|
||||
cx[0]=c[0]; cx[1]=c[1]; cx[2]=c[2]; cx[3]=c[3]; cx[4]=c[0];
|
||||
cy[0]=c[4]; cy[1]=c[5]; cy[2]=c[6]; cy[3]=c[7]; cy[4]=c[4];
|
||||
|
||||
|
||||
TGraph *g = new TGraph(5,cx,cy);
|
||||
mg->Add(g);
|
||||
if(plotCenters){
|
||||
g = new TGraph(1,&(xPPos[getBin(x,y)]),&(yPPos[getBin(x,y)]));
|
||||
mg->Add(g);
|
||||
}
|
||||
delete[] c;
|
||||
}
|
||||
return mg;
|
||||
}
|
||||
|
||||
TMultiGraph *EtaVEL::plotLog(int stepSize, int maxIt){
|
||||
int mIt;
|
||||
TMultiGraph *mg = new TMultiGraph();
|
||||
double **xposl = new double*[nPixels*nPixels+1];
|
||||
double **yposl = new double*[nPixels*nPixels+1];
|
||||
if(maxIt==-1){ mIt = it; } else{ mIt = maxIt; };
|
||||
cout << "mIt " << mIt << " steps " << mIt/stepSize << endl;
|
||||
for(int x = 0; x < nPixels; x++){
|
||||
for(int y = 0; y < nPixels; y++){
|
||||
xposl[getBin(x,y)] = new double[mIt/stepSize];
|
||||
yposl[getBin(x,y)] = new double[mIt/stepSize];
|
||||
for(int i = 0; i < mIt/stepSize; i++){
|
||||
xposl[getBin(x,y)][i] = log[i*stepSize].xPos[getBin(x,y)];
|
||||
yposl[getBin(x,y)][i] = log[i*stepSize].yPos[getBin(x,y)];
|
||||
}
|
||||
TGraph *g = new TGraph(mIt/stepSize,xposl[getBin(x,y)],yposl[getBin(x,y)]);
|
||||
g->SetLineColor((x*y % 9) + 1);
|
||||
|
||||
if(x == 0) g->SetLineColor(2);
|
||||
if(y == 0) g->SetLineColor(3);
|
||||
if(x == nPixels-1) g->SetLineColor(4);
|
||||
if(y == nPixels-1) g->SetLineColor(5);
|
||||
mg->Add(g);
|
||||
}
|
||||
}
|
||||
return mg;
|
||||
}
|
||||
|
||||
void EtaVEL::serialize(ostream &o){
|
||||
// b.WriteVersion(EtaVEL::IsA());
|
||||
char del = '|';
|
||||
o << min << del;
|
||||
o << max << del;
|
||||
o << ds << del;
|
||||
o << nPixels << del;
|
||||
o << it << del;
|
||||
o << totCont << del;
|
||||
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
|
||||
o << xPPos[i] << del;
|
||||
o << yPPos[i] << del;
|
||||
}
|
||||
for(int i = 0; i < nPixels*nPixels+1; i++){
|
||||
o << binCont[i] << del;
|
||||
}
|
||||
|
||||
for(int i = 0; i < it; i++){
|
||||
o << log[i].itN << del;
|
||||
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
|
||||
o << log[i].xPos[j] << del;
|
||||
o << log[i].yPos[j] << del;
|
||||
}
|
||||
for(int j = 0; j < nPixels*nPixels+1; j++){
|
||||
o << log[i].binCont[j] << del;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void EtaVEL::deserialize(istream &is){
|
||||
delete[] xPPos;
|
||||
delete[] yPPos;
|
||||
delete[] binCont;
|
||||
|
||||
char del;
|
||||
|
||||
is >> min >> del;
|
||||
is >> max >> del;
|
||||
is >> ds >> del;
|
||||
is >> nPixels >> del;
|
||||
is >> it >> del;
|
||||
is >> totCont >> del;
|
||||
|
||||
xPPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
yPPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
binCont = new double[nPixels*nPixels+1];
|
||||
|
||||
cout << "d";
|
||||
|
||||
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
|
||||
is >> xPPos[i] >> del;
|
||||
is >> yPPos[i] >> del;
|
||||
}
|
||||
|
||||
cout << "d";
|
||||
|
||||
for(int i = 0; i < nPixels*nPixels+1; i++){
|
||||
is >> binCont[i] >> del;
|
||||
}
|
||||
|
||||
cout << "d";
|
||||
|
||||
for(int i = 0; i < it; i++){
|
||||
is >> log[i].itN >> del;
|
||||
log[i].xPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
log[i].yPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
log[i].binCont = new double[nPixels*nPixels+1];
|
||||
|
||||
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
|
||||
is >> log[i].xPos[j] >> del;
|
||||
is >> log[i].yPos[j] >> del;
|
||||
}
|
||||
for(int j = 0; j < nPixels*nPixels+1; j++){
|
||||
is >> log[i].binCont[j] >> del;
|
||||
}
|
||||
cout << "d";
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
void EtaVEL::Streamer(TBuffer &b){
|
||||
if (b.IsReading()) {
|
||||
Version_t v = b.ReadVersion();
|
||||
|
||||
delete[] xPPos;
|
||||
delete[] yPPos;
|
||||
delete[] binCont;
|
||||
|
||||
b >> min;
|
||||
b >> max;
|
||||
b >> ds;
|
||||
b >> nPixels;
|
||||
b >> it;
|
||||
b >> totCont;
|
||||
|
||||
xPPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
yPPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
binCont = new double[nPixels*nPixels+1];
|
||||
|
||||
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
|
||||
b >> xPPos[i];
|
||||
b >> yPPos[i];
|
||||
}
|
||||
for(int i = 0; i < nPixels*nPixels+1; i++){
|
||||
b >> binCont[i];
|
||||
}
|
||||
|
||||
for(int i = 0; i < it; i++){
|
||||
b >> log[i].itN;
|
||||
log[i].xPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
log[i].yPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
log[i].binCont = new double[nPixels*nPixels+1];
|
||||
|
||||
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
|
||||
b >> log[i].xPos[j];
|
||||
b >> log[i].yPos[j];
|
||||
}
|
||||
for(int j = 0; j < nPixels*nPixels+1; j++){
|
||||
b >> log[i].binCont[j];
|
||||
}
|
||||
}
|
||||
|
||||
} else {
|
||||
b.WriteVersion(EtaVEL::IsA());
|
||||
b << min;
|
||||
b << max;
|
||||
b << ds;
|
||||
b << nPixels;
|
||||
b << it;
|
||||
b << totCont;
|
||||
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
|
||||
b << xPPos[i];
|
||||
b << yPPos[i];
|
||||
}
|
||||
for(int i = 0; i < nPixels*nPixels+1; i++){
|
||||
b << binCont[i];
|
||||
}
|
||||
|
||||
for(int i = 0; i < it; i++){
|
||||
b << log[i].itN;
|
||||
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
|
||||
b << log[i].xPos[j];
|
||||
b << log[i].yPos[j];
|
||||
}
|
||||
for(int j = 0; j < nPixels*nPixels+1; j++){
|
||||
b << log[i].binCont[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
164
slsDetectorCalibration/etaVEL/EtaVEL.h
Normal file
164
slsDetectorCalibration/etaVEL/EtaVEL.h
Normal file
@ -0,0 +1,164 @@
|
||||
#include <iostream>
|
||||
#include <TGraph.h>
|
||||
#include <TAxis.h>
|
||||
#include <TMultiGraph.h>
|
||||
#include <TH2D.h>
|
||||
#include <TMath.h>
|
||||
#include <TObject.h>
|
||||
#include <TBuffer.h>
|
||||
|
||||
#include <TMatrixD.h>
|
||||
|
||||
#include <TDecompSVD.h>
|
||||
//#include <TDecompQRH.h>
|
||||
|
||||
|
||||
#include <TH1.h>
|
||||
#include <TMath.h>
|
||||
#include <vector>
|
||||
|
||||
#include <ostream>
|
||||
#include <istream>
|
||||
|
||||
using namespace std;
|
||||
|
||||
#ifndef ETAVPS
|
||||
#define ETAVPS
|
||||
|
||||
typedef struct {
|
||||
int itN;
|
||||
double *xPos;
|
||||
double *yPos;
|
||||
double *binCont;
|
||||
} itLog;
|
||||
|
||||
|
||||
|
||||
class EtaVEL : public TObject{
|
||||
|
||||
public:
|
||||
EtaVEL(int numberOfPixels = 25, double minn=0., double maxx=1., int nnx=160, int nny=160) : nPixels(numberOfPixels), min(minn), max(maxx), converged(0), nx(nnx), ny(nny), chi_sq(0){
|
||||
//acc = 0.02;
|
||||
ds = 0.005;
|
||||
|
||||
init();
|
||||
}
|
||||
void init(){
|
||||
double pOffset = (max-min)/(double)nPixels;
|
||||
xPPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
yPPos = new double[(nPixels+1)*(nPixels+1)+1];
|
||||
binCont = new double[nPixels*nPixels+1];
|
||||
totCont = 0.;
|
||||
edgeL = new double[2*nPixels*(nPixels+1)+1];
|
||||
|
||||
for(int ii = 0; ii < 2*nPixels*(nPixels+1)+1; ii++){
|
||||
edgeL[ii] = 1.0;
|
||||
//cout << "ii " << ii << endl;
|
||||
}
|
||||
|
||||
for(int x = 0; x < nPixels+1; x++){
|
||||
for(int y = 0; y < nPixels+1; y++){
|
||||
xPPos[getCorner(x,y)] = min + (double)x * pOffset;
|
||||
yPPos[getCorner(x,y)] = min + (double)y * pOffset;
|
||||
|
||||
if(x < nPixels && y < nPixels) binCont[getBin(x,y)] = 0;
|
||||
}
|
||||
}
|
||||
// edgeL[1] = 3.0;
|
||||
updatePixelCorner();
|
||||
it = 0;
|
||||
|
||||
log = new itLog[nIterations];
|
||||
}
|
||||
|
||||
void fill(double x, double y, double amount = 1.){
|
||||
totCont+=amount;
|
||||
int bin = findBin(x,y);
|
||||
if(bin < 0) {
|
||||
//cout << "can not find bin x: " << x << " y: " << y << endl;
|
||||
totCont-=amount;
|
||||
}
|
||||
binCont[bin]+=amount;
|
||||
|
||||
}
|
||||
|
||||
int getBin(int x, int y){
|
||||
if(x < 0 || x >= nPixels || y < 0 || y >= nPixels){
|
||||
//cout << "getBin: out of bounds : x " << x << " y " << y << endl;
|
||||
return 0;
|
||||
}
|
||||
return y*nPixels+x+1;
|
||||
}
|
||||
|
||||
int getXBin(int bin){
|
||||
return (bin-1)%nPixels;
|
||||
}
|
||||
|
||||
int getYBin(int bin){
|
||||
return (bin-1)/nPixels;
|
||||
}
|
||||
|
||||
int getCorner(int x, int y){
|
||||
return y*(nPixels+1)+x+1;
|
||||
}
|
||||
|
||||
int getEdgeX(int x,int row){
|
||||
int ret = row*nPixels+x+1;
|
||||
//cout << "| edge X x " << x << " row " << row << ": "<< ret << " | ";
|
||||
return ret;
|
||||
}
|
||||
|
||||
int getEdgeY(int col, int y){
|
||||
int ret = nPixels*(nPixels+1)+col*nPixels+y+1;
|
||||
//cout << "| edge Y col " << col << " y " << y << ": "<< ret << " | ";
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
int getIt(){ return it; };
|
||||
|
||||
int getNPixels(){ return nPixels; }
|
||||
double *getXPPos(){ return xPPos; }
|
||||
double *getYPPos(){ return yPPos; }
|
||||
|
||||
void updatePixelCorner();
|
||||
double *getPixelCorners(int x, int y);
|
||||
int findBin(double xx, double yy);
|
||||
void createLogEntry();
|
||||
|
||||
void updatePixelPos();
|
||||
double *getSizeMap();
|
||||
double *getChangeMap();
|
||||
TH2D *getContent(int it=-1, int changeType = 0);
|
||||
TMultiGraph *plotPixelBorder(int plotCenters=0);
|
||||
TMultiGraph *plotLog(int stepSize=1, int maxIt=-1);
|
||||
void printGrid();
|
||||
TH1D *getCounts();
|
||||
|
||||
void serialize(ostream &o);
|
||||
void deserialize(istream &is);
|
||||
|
||||
int converged ;
|
||||
double getChiSq(){return chi_sq;};
|
||||
|
||||
private:
|
||||
itLog *log;
|
||||
int it;
|
||||
const static int nIterations =10000;
|
||||
int nx, ny;
|
||||
int nPixels;
|
||||
double *xPPos;
|
||||
double *yPPos;
|
||||
double *binCont;
|
||||
double totCont;
|
||||
double *edgeL;
|
||||
// double acc;
|
||||
double ds;
|
||||
double min,max;
|
||||
double chi_sq;
|
||||
|
||||
ClassDefNV(EtaVEL,1);
|
||||
#pragma link C++ class EtaVEL-;
|
||||
};
|
||||
|
||||
#endif
|
||||
3109
slsDetectorCalibration/etaVEL/interpolationMacro.C
Normal file
3109
slsDetectorCalibration/etaVEL/interpolationMacro.C
Normal file
File diff suppressed because it is too large
Load Diff
52
slsDetectorCalibration/etaVEL/interpolation_EtaVEL.h
Normal file
52
slsDetectorCalibration/etaVEL/interpolation_EtaVEL.h
Normal file
@ -0,0 +1,52 @@
|
||||
#ifndef INTERPOLATION_ETAVEL_H
|
||||
#define INTERPOLATION_ETAVEL_H
|
||||
|
||||
#include <slsInterpolation.h>
|
||||
#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
|
||||
134
slsDetectorCalibration/etaVEL/interpolation_etaVEL.cpp
Normal file
134
slsDetectorCalibration/etaVEL/interpolation_etaVEL.cpp
Normal file
@ -0,0 +1,134 @@
|
||||
#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++<maxit) {
|
||||
cout << " -------------- new step "<< nit << endl;
|
||||
iterate();
|
||||
}
|
||||
if (plot) {
|
||||
Draw();
|
||||
gPad->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; ibx<heta->GetNbinsX(); ibx++) {
|
||||
for (int iby=0; iby<heta->GetNbinsY(); 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);
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
224
slsDetectorCalibration/etaVEL/slsInterpolation.h
Normal file
224
slsDetectorCalibration/etaVEL/slsInterpolation.h
Normal file
@ -0,0 +1,224 @@
|
||||
#ifndef SLS_INTERPOLATION_H
|
||||
#define SLS_INTERPOLATION_H
|
||||
|
||||
#include <TObject.h>
|
||||
#include <TTree.h>
|
||||
#include <TH2F.h>
|
||||
|
||||
#ifndef DEF_QUAD
|
||||
#define DEF_QUAD
|
||||
enum quadrant {
|
||||
TOP_LEFT=0,
|
||||
TOP_RIGHT=1,
|
||||
BOTTOM_LEFT=2,
|
||||
BOTTOM_RIGHT=3,
|
||||
UNDEFINED_QUADRANT=-1
|
||||
};
|
||||
#endif
|
||||
|
||||
class slsInterpolation : public TObject{
|
||||
|
||||
public:
|
||||
slsInterpolation(int nx=40, int ny=160, int ns=25) :nPixelsX(nx), nPixelsY(ny), nSubPixels(ns) {hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);};
|
||||
|
||||
//create eta distribution, eta rebinnining etc.
|
||||
//returns flat field image
|
||||
virtual void prepareInterpolation(int &ok)=0;
|
||||
|
||||
//create interpolated image
|
||||
//returns interpolated image
|
||||
virtual TH2F *getInterpolatedImage(){return hint;};
|
||||
//return position inside the pixel for the given photon
|
||||
virtual void getInterpolatedPosition(Int_t x, Int_t y, Double_t *data, Double_t &int_x, Double_t &int_y)=0;
|
||||
|
||||
TH2F *addToImage(Double_t int_x, Double_t int_y){hint->Fill(int_x, int_y); return hint;};
|
||||
|
||||
|
||||
|
||||
virtual int addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay)=0;
|
||||
virtual int addToFlatField(Double_t etax, Double_t etay)=0;
|
||||
|
||||
//virtual void Streamer(TBuffer &b);
|
||||
|
||||
|
||||
static int calcQuad(Double_t *cl, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]){
|
||||
|
||||
int corner = UNDEFINED_QUADRANT;
|
||||
Double_t *cluster[3];
|
||||
cluster[0]=cl;
|
||||
cluster[1]=cl+3;
|
||||
cluster[2]=cl+6;
|
||||
|
||||
sum = cluster[0][0] + cluster[1][0] + cluster[2][0] + cluster[0][1] + cluster[1][1] + cluster[2][1] + cluster[0][2] + cluster[1][2] + cluster[2][2];
|
||||
|
||||
double sumBL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL
|
||||
double sumTL = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL
|
||||
double sumBR = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR
|
||||
double sumTR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR
|
||||
double sumMax = 0;
|
||||
double t, r;
|
||||
|
||||
// if(sumTL >= sumMax){
|
||||
sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0];
|
||||
sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1];
|
||||
corner = BOTTOM_LEFT;
|
||||
sumMax=sumBL;
|
||||
// }
|
||||
|
||||
if(sumTL >= sumMax){
|
||||
sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0];
|
||||
sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1];
|
||||
|
||||
corner = TOP_LEFT;
|
||||
sumMax=sumTL;
|
||||
}
|
||||
|
||||
if(sumBR >= sumMax){
|
||||
sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1];
|
||||
sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2];
|
||||
|
||||
corner = BOTTOM_RIGHT;
|
||||
sumMax=sumBR;
|
||||
}
|
||||
|
||||
if(sumTR >= sumMax){
|
||||
sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1];
|
||||
sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2];
|
||||
|
||||
corner = TOP_RIGHT;
|
||||
sumMax=sumTR;
|
||||
}
|
||||
|
||||
totquad=sumMax;
|
||||
|
||||
return corner;
|
||||
|
||||
}
|
||||
|
||||
static int calcEta(Double_t totquad, Double_t sDum[2][2], Double_t &etax, Double_t &etay){
|
||||
Double_t t,r;
|
||||
|
||||
if (totquad>0) {
|
||||
t = sDum[1][0] + sDum[1][1];
|
||||
r = sDum[0][1] + sDum[1][1];
|
||||
etax=r/totquad;
|
||||
etay=t/totquad;
|
||||
}
|
||||
return 0;
|
||||
|
||||
}
|
||||
|
||||
|
||||
static int calcEta(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]) {
|
||||
int corner = calcQuad(cl,sum,totquad,sDum);
|
||||
calcEta(totquad, sDum, etax, etay);
|
||||
|
||||
return corner;
|
||||
}
|
||||
|
||||
|
||||
static int calcEtaL(Double_t totquad, int corner, Double_t sDum[2][2], Double_t &etax, Double_t &etay){
|
||||
Double_t t,r, toth, totv;
|
||||
if (totquad>0) {
|
||||
switch(corner) {
|
||||
case TOP_LEFT:
|
||||
t = sDum[1][1] ;
|
||||
r = sDum[0][1] ;
|
||||
toth=sDum[1][1]+sDum[1][0];
|
||||
totv=sDum[0][1]+sDum[1][1];
|
||||
break;
|
||||
case TOP_RIGHT:
|
||||
t = sDum[1][0] ;
|
||||
r = sDum[0][1] ;
|
||||
toth=sDum[0][0]+t;
|
||||
totv=sDum[0][0]+r;
|
||||
break;
|
||||
case BOTTOM_LEFT:
|
||||
r = sDum[1][1] ;
|
||||
t = sDum[1][1] ;
|
||||
toth=sDum[1][0]+t;
|
||||
totv=sDum[0][1]+r;
|
||||
break;
|
||||
case BOTTOM_RIGHT:
|
||||
t = sDum[1][0] ;
|
||||
r = sDum[1][1] ;
|
||||
toth=sDum[1][1]+t;
|
||||
totv=sDum[0][1]+r;
|
||||
break;
|
||||
default:
|
||||
etax=-1;
|
||||
etay=-1;
|
||||
return 0;
|
||||
}
|
||||
etax=r/totv;
|
||||
etay=t/toth;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int calcEtaL(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]) {
|
||||
int corner = calcQuad(cl,sum,totquad,sDum);
|
||||
calcEtaL(totquad, corner, sDum, etax, etay);
|
||||
|
||||
return corner;
|
||||
}
|
||||
|
||||
|
||||
|
||||
static int calcEtaC3(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]){
|
||||
|
||||
int corner = calcQuad(cl,sum,totquad,sDum);
|
||||
calcEta(sum, sDum, etax, etay);
|
||||
return corner;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
static int calcEta3(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum) {
|
||||
Double_t l,r,t,b;
|
||||
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
|
||||
if (sum>0) {
|
||||
l=cl[0]+cl[3]+cl[6];
|
||||
r=cl[2]+cl[5]+cl[8];
|
||||
b=cl[0]+cl[1]+cl[2];
|
||||
t=cl[6]+cl[7]+cl[8];
|
||||
etax=(-l+r)/sum;
|
||||
etay=(-b+t)/sum;
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
|
||||
static int calcEta3X(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum) {
|
||||
Double_t l,r,t,b;
|
||||
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
|
||||
if (sum>0) {
|
||||
l=cl[3];
|
||||
r=cl[5];
|
||||
b=cl[1];
|
||||
t=cl[7];
|
||||
etax=(-l+r)/sum;
|
||||
etay=(-b+t)/sum;
|
||||
}
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
protected:
|
||||
int nPixelsX, nPixelsY;
|
||||
int nSubPixels;
|
||||
TH2F *hint;
|
||||
|
||||
|
||||
// ClassDefNV(slsInterpolation,1);
|
||||
// #pragma link C++ class slsInterpolation-;
|
||||
};
|
||||
|
||||
#endif
|
||||
154
slsDetectorCalibration/moench03TCtbData.h
Normal file
154
slsDetectorCalibration/moench03TCtbData.h
Normal file
@ -0,0 +1,154 @@
|
||||
#ifndef MOENCH03TCTBDATA_H
|
||||
#define MOENCH03TCTBDATA_H
|
||||
#include "slsDetectorData.h"
|
||||
|
||||
|
||||
|
||||
class moench03CtbData : public slsDetectorData<uint16_t> {
|
||||
|
||||
private:
|
||||
|
||||
int iframe;
|
||||
int nadc;
|
||||
int sc_width;
|
||||
int sc_height;
|
||||
public:
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
Implements the slsReceiverData structure for the moench02 prototype read out by a module i.e. using the slsReceiver
|
||||
(160x160 pixels, 40 packets 1286 large etc.)
|
||||
\param c crosstalk parameter for the output buffer
|
||||
|
||||
*/
|
||||
|
||||
|
||||
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]={300,325,350,375,300,325,350,375, \
|
||||
200,225,250,275,200,225,250,275,\
|
||||
100,125,150,175,100,125,150,175,\
|
||||
0,25,50,75,0,25,50,75};
|
||||
int row, col;
|
||||
|
||||
int isample;
|
||||
int iadc;
|
||||
int ix, iy;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
for (iadc=0; iadc<nadc; iadc++) {
|
||||
for (int i=0; i<sc_width*sc_height; i++) {
|
||||
col=adc_nr[iadc]+(i%sc_width);
|
||||
if (iadc<16) {
|
||||
row=199-i/sc_width;
|
||||
} else {
|
||||
row=200+i/sc_width;
|
||||
}
|
||||
dataMap[row][col]=(nadc*i+iadc)*2;
|
||||
if (dataMap[row][col]<0 || dataMap[row][col]>=2*400*400)
|
||||
cout << "Error: pointer " << dataMap[row][col] << " out of range "<< endl;
|
||||
|
||||
}
|
||||
}
|
||||
for (int i=0; i<nx*ny; i++) {
|
||||
isample=i/nadc;
|
||||
iadc=i%nadc;
|
||||
ix=isample%sc_width;
|
||||
iy=isample/sc_width;
|
||||
if (iadc<(nadc/2)) {
|
||||
xmap[i]=adc_nr[iadc]+ix;
|
||||
ymap[i]=ny/2-1-iy;
|
||||
} else {
|
||||
xmap[i]=adc_nr[iadc]+ix;
|
||||
ymap[i]=ny/2+iy;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
iframe=0;
|
||||
// cout << "data struct created" << endl;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
|
||||
Returns the frame number for the given dataset. Purely virtual func.
|
||||
\param buff pointer to the dataset
|
||||
\returns frame number
|
||||
|
||||
*/
|
||||
|
||||
|
||||
virtual int getFrameNumber(char *buff){(void)buff; return iframe;};
|
||||
|
||||
/**
|
||||
|
||||
Returns the packet number for the given dataset. purely virtual func
|
||||
\param buff pointer to the dataset
|
||||
\returns packet number number
|
||||
|
||||
|
||||
virtual int getPacketNumber(char *buff)=0;
|
||||
|
||||
*/
|
||||
|
||||
/**
|
||||
|
||||
Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func
|
||||
\param data pointer to the memory to be analyzed
|
||||
\param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot
|
||||
\param dsize size of the memory slot to be analyzed
|
||||
\returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found
|
||||
|
||||
*/
|
||||
virtual char *findNextFrame(char *data, int &ndata, int dsize){ndata=dsize; setDataSize(dsize); return data;};
|
||||
|
||||
|
||||
/**
|
||||
|
||||
Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
|
||||
\param filebin input file stream (binary)
|
||||
\returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete
|
||||
|
||||
*/
|
||||
virtual char *readNextFrame(ifstream &filebin){
|
||||
// int afifo_length=0;
|
||||
uint16_t *afifo_cont;
|
||||
int ib=0;
|
||||
if (filebin.is_open()) {
|
||||
afifo_cont=new uint16_t[dataSize/2];
|
||||
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;
|
||||
return NULL;
|
||||
}
|
||||
}
|
||||
return NULL;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
@ -228,7 +228,7 @@ class slsDetectorData {
|
||||
\returns data for the selected channel, with inversion if required as double
|
||||
|
||||
*/
|
||||
virtual double getValue(char *data, int ix, int iy=0) {return (double)getChannel(data, ix, iy);};
|
||||
virtual double getValue(char *data, int ix, int iy=0) {/* cout << " x "<< ix << " y"<< iy << " val " << getChannel(data, ix, iy)<< endl;*/return (double)getChannel(data, ix, iy);};
|
||||
|
||||
|
||||
/**
|
||||
|
||||
@ -97,7 +97,7 @@ public:
|
||||
break;
|
||||
} else {
|
||||
//cprintf(BG_RED, "Too many packets for this frame! fnum:%d, pnum:%d np:%d\n",fnum,pnum,np);
|
||||
cout << "Too many packets for this frame! "<< fnum << " " << pnum << endl;cprintf(BG_RED,"Exiting\n");exit(-1);
|
||||
cout << "Too many packets for this frame! "<< fnum << " " << pnum << endl;//cprintf(BG_RED,"Exiting\n");exit(-1);
|
||||
retval=NULL;
|
||||
}
|
||||
}
|
||||
@ -105,7 +105,7 @@ public:
|
||||
if (np<nPackets) {
|
||||
if (np>0){
|
||||
//cprintf(BG_RED, "Too few packets for this frame! fnum:%d, pnum:%d np:%d\n",fnum,pnum,np);
|
||||
cout << "Too few packets for this frame! "<< fnum << " " << pnum << " " << np <<endl;cprintf(BG_RED,"Exiting\n");exit(-1);
|
||||
cout << "Too few packets for this frame! "<< fnum << " " << pnum << " " << np <<endl;//cprintf(BG_RED,"Exiting\n");exit(-1);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user