mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-06-21 00:58:01 +02:00
added patternGenerator and slsDteectorCalibration directory in order to compile the ctbGui
This commit is contained in:
678
slsDetectorCalibration/interpolations/etaVEL/EVELAlg.C
Normal file
678
slsDetectorCalibration/interpolations/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/interpolations/etaVEL/EtaVEL.cpp
Normal file
679
slsDetectorCalibration/interpolations/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/interpolations/etaVEL/EtaVEL.h
Normal file
164
slsDetectorCalibration/interpolations/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
|
||||
393
slsDetectorCalibration/interpolations/etaVEL/EtaVELTr.py
Normal file
393
slsDetectorCalibration/interpolations/etaVEL/EtaVELTr.py
Normal file
@ -0,0 +1,393 @@
|
||||
import numpy as np
|
||||
import math
|
||||
|
||||
maxf = 2
|
||||
minf = 0.5
|
||||
|
||||
class EtaVELTr:
|
||||
def __init__(self, numberOfPixels = 25, minn=0., maxx=1.):
|
||||
self.nPixels = numberOfPixels
|
||||
self.nCorners = self.nPixels + 1
|
||||
self.minEta = minn
|
||||
self.maxEta = maxx
|
||||
#self.corners = []
|
||||
self.edgesX = []
|
||||
self.edgesY = []
|
||||
self.edgesE = []
|
||||
self.edgesF = []
|
||||
self.edgesG = []
|
||||
self.edgesH = []
|
||||
self.counts = []
|
||||
self.pOffset = (self.maxEta-self.minEta)/self.nPixels
|
||||
|
||||
self.cPosX = []
|
||||
self.cPosY = []
|
||||
self.zPosX = []
|
||||
self.zPosY = []
|
||||
|
||||
self.sqSums = []
|
||||
self.cIteration = 0
|
||||
self.bSqSum = 0
|
||||
|
||||
self.initGrid()
|
||||
#self.calculatePixelCorners()
|
||||
self.update()
|
||||
|
||||
def initGrid(self):
|
||||
dd = 1 / math.sqrt(2)
|
||||
|
||||
#self.corners = [ [self.minEta + x * pOffset, self.minEta + y * pOffset] for y in range(self.nPixels)] for x in range(self.nPixels)
|
||||
self.cPosX = [ [self.minEta + x * self.pOffset for x in range(self.nCorners)] for y in range(self.nCorners) ]
|
||||
self.cPosY = [ [self.minEta + y * self.pOffset for x in range(self.nCorners)] for y in range(self.nCorners) ]
|
||||
self.counts = [ [ [0,0,0,0] for x in range(self.nPixels) ] for y in range(self.nPixels) ]
|
||||
self.edgesX = [ [ 1 for x in range(self.nCorners) ] for y in range(self.nCorners + 1) ]
|
||||
self.edgesY = [ [ 1 for x in range(self.nCorners+1) ] for y in range(self.nCorners) ]
|
||||
self.edgesE = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
|
||||
self.edgesF = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
|
||||
self.edgesG = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
|
||||
self.edgesH = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
|
||||
|
||||
|
||||
|
||||
def update(self):
|
||||
self.normalizeEdgeLengths()
|
||||
self.calculateEdgeLengths2()
|
||||
self.calculatePixelCorners2()
|
||||
conv = False
|
||||
out = 0
|
||||
outList = []
|
||||
tot = self.nPixels*self.nPixels*4
|
||||
sqSum = 0
|
||||
avg = self.getAvgCounts()
|
||||
avgPS = self.getAvgCounts() + 1*math.sqrt(self.getAvgCounts())
|
||||
avgMS = self.getAvgCounts() - 1*math.sqrt(self.getAvgCounts())
|
||||
for y in range(self.nPixels):
|
||||
for x in range(self.nPixels):
|
||||
for t in range(4):
|
||||
sqSum += (avg -self.counts[y][x][t]) * (avg -self.counts[y][x][t])
|
||||
if self.counts[y][x][t] > avgPS or self.counts[y][x][t] < avgMS:
|
||||
|
||||
out += 1
|
||||
outList.append([y,x,t,self.counts[y][x][t]])
|
||||
|
||||
outList = sorted(outList,key=lambda t: abs(self.counts[t[0]][t[1]][t[2]]/self.getAvgCounts()))
|
||||
self.counts = [ [ [0,0,0,0] for x in range(self.nPixels) ] for y in range(self.nPixels) ]
|
||||
print("There are {} of {} triangles out of 1 std ({} %)".format(out,tot,out/tot*100))
|
||||
print("Total Sq Err: {}".format(sqSum))
|
||||
|
||||
self.sqSums.append(sqSum)
|
||||
|
||||
if len(self.sqSums) > 2 and np.diff(self.sqSums)[-1] > 0:
|
||||
print("converged after {} steps: sqSums {} diff(sqSums) {}".format(self.cIteration,self.sqSums,np.diff(self.sqSums)))
|
||||
conv = True
|
||||
self.bSqSum = self.sqSums[-2]
|
||||
|
||||
self.cIteration += 1
|
||||
return [conv,outList,sqSum]
|
||||
|
||||
def normalizeEdgeLengths(self):
|
||||
sumL = 0
|
||||
sumL += sum(map(sum,zip(*self.edgesX)))
|
||||
sumL += sum(map(sum,zip(*self.edgesY)))
|
||||
sumL += sum(map(sum,zip(*self.edgesE)))
|
||||
sumL += sum(map(sum,zip(*self.edgesF)))
|
||||
sumL += sum(map(sum,zip(*self.edgesG)))
|
||||
sumL += sum(map(sum,zip(*self.edgesH)))
|
||||
avgL = sumL/(4*self.nPixels*self.nPixels+2*self.nCorners*(self.nCorners+1))
|
||||
print("total Sum is {} avg: {}".format(sumL,avgL))
|
||||
self.edgesX = [ [ x/avgL for x in y ] for y in self.edgesX ]
|
||||
self.edgesY = [ [ x/avgL for x in y ] for y in self.edgesY ]
|
||||
self.edgesE = [ [ x/avgL for x in y ] for y in self.edgesE ]
|
||||
self.edgesF = [ [ x/avgL for x in y ] for y in self.edgesF ]
|
||||
self.edgesG = [ [ x/avgL for x in y ] for y in self.edgesG ]
|
||||
self.edgesH = [ [ x/avgL for x in y ] for y in self.edgesH ]
|
||||
|
||||
def _shapeF(self,f):
|
||||
f = (f - 1) * 0.6 + 1
|
||||
if f > maxf:
|
||||
return maxf
|
||||
if f < minf:
|
||||
return minf
|
||||
return f
|
||||
|
||||
def calculateEdgeLengths2(self):
|
||||
if self.getTotalCounts() == 0:
|
||||
return
|
||||
avg = self.getAvgCounts()
|
||||
|
||||
for y in range(self.nPixels):
|
||||
for x in range(self.nPixels):
|
||||
for t in range(4):
|
||||
pc = self.counts[y][x][t]
|
||||
if pc == 0:
|
||||
f = maxf
|
||||
else:
|
||||
f = math.sqrt(avg/pc)
|
||||
if pc > avg-math.sqrt(avg) and pc < avg+math.sqrt(avg):
|
||||
f = 1.
|
||||
sf = self._shapeF(f)
|
||||
if t == 0:
|
||||
self.edgesX[y][x] = self.edgesX[y][x] / sf
|
||||
self.edgesE[y][x] = self.edgesE[y][x] / sf
|
||||
self.edgesF[y][x] = self.edgesF[y][x] / sf
|
||||
if t == 1:
|
||||
self.edgesY[y][x+1] = self.edgesY[y][x+1] / sf
|
||||
self.edgesF[y][x] = self.edgesF[y][x] / sf
|
||||
self.edgesH[y][x] = self.edgesH[y][x] / sf
|
||||
if t == 2:
|
||||
self.edgesX[y+1][x] = self.edgesX[y+1][x] / sf
|
||||
self.edgesH[y][x] = self.edgesH[y][x] / sf
|
||||
self.edgesG[y][x] = self.edgesG[y][x] / sf
|
||||
if t == 3:
|
||||
self.edgesY[y][x] = self.edgesY[y][x] / sf
|
||||
self.edgesG[y][x] = self.edgesG[y][x] / sf
|
||||
self.edgesE[y][x] = self.edgesE[y][x] / sf
|
||||
|
||||
|
||||
def calculatePixelCorners2(self):
|
||||
w = 20
|
||||
posMat = []
|
||||
CrVx = np.zeros((self.nCorners,self.nCorners))
|
||||
CrVy = np.zeros((self.nCorners,self.nCorners))
|
||||
ZrVx = np.zeros((self.nPixels,self.nPixels))
|
||||
ZrVy = np.zeros((self.nPixels,self.nPixels))
|
||||
|
||||
#boundary conditions matrix/vectors
|
||||
BCposMatX = []
|
||||
BCposMatY = []
|
||||
BCrVx = []
|
||||
BCrVy = []
|
||||
|
||||
for y in range(self.nCorners):
|
||||
for x in range(self.nCorners):
|
||||
BClineX = np.zeros((self.nCorners,self.nCorners))
|
||||
BClineY = np.zeros((self.nCorners,self.nCorners))
|
||||
if (x == 0 and y == 0) or \
|
||||
(x == 0 and y == self.nPixels) or \
|
||||
(x == self.nPixels and y == 0) or \
|
||||
(x == self.nPixels and y == self.nPixels):
|
||||
BClineX[y][x] = w
|
||||
BClineY[y][x] = w
|
||||
BCrVx.append(self.getCornerPos(y,x)[0] * w)
|
||||
BCrVy.append(self.getCornerPos(y,x)[1] * w)
|
||||
#print("bclinex shape {} zeros shape {}".format( BClineX.reshape((self.nCorners*self.nCorners,)).shape , np.zeros((self.nPixels*self.nPixels)).shape ))
|
||||
BCposMatX.append(np.hstack((BClineX.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
|
||||
BCposMatY.append(np.hstack((BClineY.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
|
||||
|
||||
elif x == 0 or x == self.nPixels:
|
||||
BClineX[y][x] = w
|
||||
#BClineY[y][x] = 1
|
||||
BCrVx.append(self.getCornerPos(y,x)[0] * w)
|
||||
#BCrVy.append(self.getCornerPos(y,x)[1])
|
||||
BCposMatX.append(np.hstack((BClineX.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
|
||||
#BCposMatY.append(np.hstack((BClineY.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
|
||||
elif y == 0 or y == self.nPixels:
|
||||
#BClineX[y][x] = 1
|
||||
BClineY[y][x] = w
|
||||
#BCrVx.append(self.getCornerPos(y,x)[0])
|
||||
BCrVy.append(self.getCornerPos(y,x)[1] * w)
|
||||
#BCposMatX.append(np.hstack((BClineX.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
|
||||
BCposMatY.append(np.hstack((BClineY.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
|
||||
|
||||
|
||||
eLength = 0
|
||||
|
||||
if x != 0:
|
||||
eLength += self.edgesX[y][x-1]
|
||||
if y != 0:
|
||||
eLength += self.edgesY[y-1][x]
|
||||
if x != self.nPixels:
|
||||
eLength += self.edgesX[y][x]
|
||||
if y != self.nPixels:
|
||||
eLength += self.edgesY[y][x]
|
||||
|
||||
if y != 0 and x != 0:
|
||||
eLength += self.edgesH[y-1][x-1]
|
||||
if y != self.nPixels and x != 0:
|
||||
eLength += self.edgesF[y][x-1]
|
||||
if y != 0 and x != self.nPixels:
|
||||
eLength += self.edgesG[y-1][x]
|
||||
if y != self.nPixels and x != self.nPixels:
|
||||
eLength += self.edgesE[y][x]
|
||||
|
||||
line = np.zeros((self.nCorners,self.nCorners))
|
||||
lineZ = np.zeros((self.nPixels,self.nPixels))
|
||||
|
||||
|
||||
if x != 0:
|
||||
line[y][x-1] = - self.edgesX[y][x-1]/eLength
|
||||
if y != 0:
|
||||
line[y-1][x] = - self.edgesY[y-1][x]/eLength
|
||||
if x != self.nPixels:
|
||||
line[y][x+1] = - self.edgesX[y][x]/eLength
|
||||
if y != self.nPixels:
|
||||
line[y+1][x] = - self.edgesY[y][x]/eLength
|
||||
|
||||
if y != 0 and x != 0:
|
||||
lineZ[y-1][x-1] = -self.edgesH[y-1][x-1]/eLength
|
||||
if y != self.nPixels and x != 0:
|
||||
lineZ[y][x-1] = -self.edgesF[y][x-1]/eLength
|
||||
if y != 0 and x != self.nPixels:
|
||||
lineZ[y-1][x] = -self.edgesG[y-1][x]/eLength
|
||||
if y != self.nPixels and x != self.nPixels:
|
||||
lineZ[y][x] = -self.edgesE[y][x]/eLength
|
||||
|
||||
|
||||
line[y][x] = 1
|
||||
CrVx[y][x] = 0
|
||||
CrVy[y][x] = 0
|
||||
posMat.append( \
|
||||
np.hstack(( \
|
||||
line.reshape((self.nCorners*self.nCorners,)), \
|
||||
lineZ.reshape((self.nPixels*self.nPixels,)) \
|
||||
)) \
|
||||
)
|
||||
|
||||
for y in range(self.nPixels):
|
||||
for x in range(self.nPixels):
|
||||
line = np.zeros((self.nCorners,self.nCorners))
|
||||
lineZ = np.zeros((self.nPixels,self.nPixels))
|
||||
|
||||
eLength = self.edgesE[y][x] + self.edgesF[y][x] +self.edgesG[y][x] +self.edgesH[y][x]
|
||||
line[y][x] = -self.edgesE[y][x] / eLength
|
||||
line[y][x+1] = -self.edgesF[y][x] / eLength
|
||||
line[y+1][x] = -self.edgesG[y][x] / eLength
|
||||
line[y+1][x+1] = -self.edgesH[y][x] / eLength
|
||||
|
||||
lineZ[y][x] = 1
|
||||
ZrVx[y][x] = 0
|
||||
ZrVy[y][x] = 0
|
||||
posMat.append( \
|
||||
np.hstack(( \
|
||||
line.reshape((self.nCorners*self.nCorners,)), \
|
||||
lineZ.reshape((self.nPixels*self.nPixels,)) \
|
||||
)) \
|
||||
)
|
||||
|
||||
CrVxFlat = CrVx.reshape((self.nCorners*self.nCorners,))
|
||||
CrVyFlat = CrVy.reshape((self.nCorners*self.nCorners,))
|
||||
ZrVxFlat = ZrVx.reshape((self.nPixels*self.nPixels,))
|
||||
ZrVyFlat = ZrVy.reshape((self.nPixels*self.nPixels,))
|
||||
posMat = np.asarray(posMat)
|
||||
|
||||
BCrVyFlat = np.asarray(BCrVy)
|
||||
BCrVxFlat = np.asarray(BCrVx)
|
||||
BCposMatX = np.asarray(BCposMatX)
|
||||
BCposMatY = np.asarray(BCposMatY)
|
||||
|
||||
print ("BCposMatY vy {} shape posMat {}".format(BCposMatY.shape,posMat.shape))
|
||||
|
||||
FinalrVy = np.hstack((CrVyFlat,ZrVyFlat,BCrVyFlat))
|
||||
FinalrVx = np.hstack((CrVxFlat,ZrVxFlat,BCrVxFlat))
|
||||
FinalposMatX = np.vstack((posMat,BCposMatX))
|
||||
FinalposMatY = np.vstack((posMat,BCposMatY))
|
||||
|
||||
print("posMat shape {}".format(posMat.shape))
|
||||
print("posMatX shape {}".format(FinalposMatX.shape))
|
||||
print("rVxFlat shape {}".format(FinalrVx.shape))
|
||||
|
||||
#print("posMat {}".format(FinalposMat))
|
||||
#print("rVxFlat {}".format(FinalrVx))
|
||||
|
||||
xPos = np.linalg.lstsq(FinalposMatX,FinalrVx)[0]
|
||||
yPos = np.linalg.lstsq(FinalposMatY,FinalrVy)[0]
|
||||
|
||||
print("xPosShape {} cutXPosShape {}".format(xPos.shape,xPos[:self.nCorners][:self.nCorners].shape))
|
||||
|
||||
self.cPosX = xPos[:self.nCorners*self.nCorners].reshape((self.nCorners,self.nCorners))
|
||||
self.cPosY = yPos[:self.nCorners*self.nCorners].reshape((self.nCorners,self.nCorners))
|
||||
|
||||
self.zPosX = xPos[self.nCorners*self.nCorners:].reshape((self.nPixels,self.nPixels))
|
||||
self.zPosY = yPos[self.nCorners*self.nCorners:].reshape((self.nPixels,self.nPixels))
|
||||
|
||||
|
||||
|
||||
def fill(self,yy,xx,count = 1):
|
||||
[y,x,t] = self.getPixel(yy,xx)
|
||||
self.counts[y][x][t] += count
|
||||
|
||||
def getCountDist(self):
|
||||
c = []
|
||||
for y in range(self.nPixels):
|
||||
for x in range(self.nPixels):
|
||||
c.append(self.counts[y][x])
|
||||
return c
|
||||
|
||||
def getPixel(self,yy,xx, debug = False):
|
||||
for y in range(self.nPixels):
|
||||
for x in range(self.nPixels):
|
||||
for t in range(4):
|
||||
[v1x,v1y,v2x,v2y,v3x,v3y] = self.getTriangleCorner(y,x,t)
|
||||
if self.pointInTriangle([xx,yy],[v1x,v1y],[v2x,v2y],[v3x,v3y]):
|
||||
return [y,x,t]
|
||||
|
||||
if not debug:
|
||||
raise Exception("not inside a pixel")
|
||||
else:
|
||||
print("no pixel found")
|
||||
return [0,0]
|
||||
|
||||
#http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-2d-triangle
|
||||
def trSign (self, p1, p2, p3):
|
||||
return (p1[0] - p3[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p3[1]);
|
||||
|
||||
def pointInTriangle (self,pt, v1, v2, v3):
|
||||
b1 = self.trSign(pt, v1, v2) < 0.0
|
||||
b2 = self.trSign(pt, v2, v3) < 0.0
|
||||
b3 = self.trSign(pt, v3, v1) < 0.0
|
||||
return ((b1 == b2) and (b2 == b3));
|
||||
|
||||
def getAvgCounts(self):
|
||||
return self.getTotalCounts() / self.nPixels/self.nPixels/4.
|
||||
|
||||
def getTotalCounts(self):
|
||||
tot = 0
|
||||
for y in range(self.nPixels):
|
||||
for x in range(self.nPixels):
|
||||
for t in range(4):
|
||||
tot += self.counts[y][x][t]
|
||||
return tot
|
||||
|
||||
#tl tr bl br
|
||||
def getPixelCorners(self,iy,ix):
|
||||
return self.getCornerPos(iy,ix) + self.getCornerPos(iy,ix+1) + self.getCornerPos(iy+1,ix) + self.getCornerPos(iy+1,ix+1)
|
||||
|
||||
def getTriangleCorner(self,iy,ix,tr):
|
||||
if tr == 0:
|
||||
return self.getCornerPos(iy,ix) + self.getCornerPos(iy,ix+1) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
if tr == 1:
|
||||
return self.getCornerPos(iy,ix+1) + self.getCornerPos(iy+1,ix+1) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
if tr == 2:
|
||||
return self.getCornerPos(iy+1,ix+1) + self.getCornerPos(iy+1,ix) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
if tr == 3:
|
||||
return self.getCornerPos(iy+1,ix) + self.getCornerPos(iy,ix) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
|
||||
def getCornerPos(self,iy,ix):
|
||||
return [self.cPosX[iy][ix],self.cPosY[iy][ix]]
|
||||
|
||||
def getXEdgePos(self,iy,ix):
|
||||
p1 = self.getCornerPos(iy,ix)
|
||||
p2 = self.getCornerPos(iy,ix+1)
|
||||
return [p1[0],p1[1],p2[0],p2[1]]
|
||||
|
||||
def getYEdgePos(self,iy,ix):
|
||||
p1 = self.getCornerPos(iy,ix)
|
||||
p2 = self.getCornerPos(iy+1,ix)
|
||||
return [p1[0],p1[1],p2[0],p2[1]]
|
||||
|
||||
def getEEdgePos(self,iy,ix):
|
||||
p1 = self.getCornerPos(iy,ix)
|
||||
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
|
||||
def getFEdgePos(self,iy,ix):
|
||||
p1 = self.getCornerPos(iy,ix+1)
|
||||
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
|
||||
def getGEdgePos(self,iy,ix):
|
||||
p1 = self.getCornerPos(iy+1,ix)
|
||||
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
|
||||
def getHEdgePos(self,iy,ix):
|
||||
p1 = self.getCornerPos(iy+1,ix+1)
|
||||
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
|
||||
|
||||
@ -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);
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
@ -0,0 +1,55 @@
|
||||
#ifndef INTERPOLATION_ETAVEL_H
|
||||
#define INTERPOLATION_ETAVEL_H
|
||||
|
||||
#include <slsInterpolation.h>
|
||||
#include "EtaVEL.h"
|
||||
//#include "TH2F.h"
|
||||
//#include "EtaVEL.cpp"
|
||||
//class EtaVEL;
|
||||
|
||||
class etaVELInterpolation: public etaInterpolationBase {
|
||||
|
||||
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
|
||||
Reference in New Issue
Block a user