diff --git a/slsDetectorCalibration/energyCalibration.cpp b/slsDetectorCalibration/energyCalibration.cpp index a4c78ffcf..596d9a0e4 100644 --- a/slsDetectorCalibration/energyCalibration.cpp +++ b/slsDetectorCalibration/energyCalibration.cpp @@ -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; diff --git a/slsDetectorCalibration/energyCalibration.h b/slsDetectorCalibration/energyCalibration.h index 7afda027a..f060b8d1a 100644 --- a/slsDetectorCalibration/energyCalibration.h +++ b/slsDetectorCalibration/energyCalibration.h @@ -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 diff --git a/slsDetectorCalibration/etaVEL/EVELAlg.C b/slsDetectorCalibration/etaVEL/EVELAlg.C new file mode 100644 index 000000000..5ae9a0c51 --- /dev/null +++ b/slsDetectorCalibration/etaVEL/EVELAlg.C @@ -0,0 +1,678 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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; ilGetEntryAndTree(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; ilGetEntryAndTree(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(); +} diff --git a/slsDetectorCalibration/etaVEL/EtaVEL.cpp b/slsDetectorCalibration/etaVEL/EtaVEL.cpp new file mode 100644 index 000000000..87f901cbd --- /dev/null +++ b/slsDetectorCalibration/etaVEL/EtaVEL.cpp @@ -0,0 +1,679 @@ +#include "EtaVEL.h" +#include + + +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 (xxmax) xx=max-1E-6; + // if (yymax) 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:" <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:" < 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]; + } + } + } +} + diff --git a/slsDetectorCalibration/etaVEL/EtaVEL.h b/slsDetectorCalibration/etaVEL/EtaVEL.h new file mode 100644 index 000000000..5c73d8a83 --- /dev/null +++ b/slsDetectorCalibration/etaVEL/EtaVEL.h @@ -0,0 +1,164 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +//#include + + +#include +#include +#include + +#include +#include + +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 diff --git a/slsDetectorCalibration/etaVEL/interpolationMacro.C b/slsDetectorCalibration/etaVEL/interpolationMacro.C new file mode 100644 index 000000000..aebbb04d7 --- /dev/null +++ b/slsDetectorCalibration/etaVEL/interpolationMacro.C @@ -0,0 +1,3109 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +//#include "interpolation_EtaVEL.h" +#include "interpolation_etaVEL.cpp" +#include "../energyCalibration.cpp" + +/* +Zum erstellen der correction map ist createGainAndEtaFile(...) in EVELAlg.C der entry point. +Zum erstellen des HR images ist createImage(...) der entry point. +*/ +#define YMIN 100 +#define XMIN 0 +#define NX 400 +#define NY 200 + +// sets nice plotting settings +int rootlogon () +{ + gStyle->SetDrawBorder(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetCanvasDefH(800); + gStyle->SetCanvasDefW(800); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetPaintTextFormat("5.2f"); + gStyle->SetLineWidth(2); + gStyle->SetTextSize(1.1); + gStyle->SetLabelSize(0.04,"xy"); + gStyle->SetTitleSize(0.05,"xy"); + gStyle->SetTitleOffset(1.0,"x"); + gStyle->SetTitleOffset(1.6,"y"); + gStyle->SetPadTopMargin(0.05); + gStyle->SetPadRightMargin(0.05); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetLegendBorderSize(1); + gStyle->SetFrameBorderMode(0); + gStyle->SetFrameFillColor(kWhite); + gStyle->SetTitleFillColor(kWhite); + gStyle->SetStatFontSize(0.03); + gStyle->SetStatBorderSize(1); + gStyle->SetStatFormat("6.4g"); + gStyle->SetStatX(0.95); + gStyle->SetStatY(0.95); + gStyle->SetStatW(0.2); + gStyle->SetStatH(0.2); + gStyle->SetStatColor(kWhite); + gStyle->SetTitleX(0.3); + gStyle->SetTitleY(0.98); + gStyle->SetTitleBorderSize(1); + gStyle->SetTitleFontSize(0.06); + gStyle->SetLegendBorderSize(1); + gROOT->SetStyle("Default"); + gROOT->ForceStyle(); + return 0; + +} + + +interpolation_EtaVEL *fillEta(TChain *ch, Long64_t nentries=-1, TH2F *gmap=NULL, THStack *hs=NULL){ + + // if (nentries<0) + // nentries=ch->GetEntries(); + + // cout << "Chain has " << nentries << " entries " << endl; + + + + Long64_t ie; + + Double_t dum[9], data[9], gain=1; + Int_t x, y, f, ix, iy; + interpolation_EtaVEL *inte= new interpolation_EtaVEL(400, 200, 25, -0.05,1.05); + int skip=0; + TEntryList *myelist=ch->GetEntryList(); + Long64_t listEntries = myelist->GetN(); + Double_t rr; + Long64_t treeEntry, chainEntry; + Int_t treenum; + TCanvas *c1=new TCanvas(); + int corn; + TH1D *h1=NULL, *h2=NULL, *h3=NULL; + TH2F *heta=NULL, *hetaL=NULL, *heta3=NULL, *heta3X=NULL, *hcorn=NULL ; + Double_t tot, totquad, sDum[2][2]; + Double_t etax, etay, etaxL, etayL, etax3, etay3,etax3x, etay3x,r ; + heta=inte->setEta((TH2F*)NULL); + + if (hs) { + hetaL=(TH2F*)(heta->Clone("hetaL")); + heta3=(TH2F*)(heta->Clone("heta3")); + heta3X=(TH2F*)(heta->Clone("heta3X")); + h1=new TH1D("h1","h1",1000,0,3); + h2=new TH1D("h2","h2",1000,0,3); + h3=new TH1D("h3","h3",1000,0,3); + hcorn=new TH2F("hcorn","hcorn",2,-1,1,2,-1,1); + hs->Add(heta); + hs->Add(hetaL); + hs->Add(heta3); + hs->Add(heta3X); + hs->Add(h1); + hs->Add(h2); + hs->Add(h3); + hs->Add(hcorn); + + } + + cout << "List has " << listEntries*1E-6 << "M entries " << endl; + + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + + for (Long64_t el = 0; el < listEntries; el++) { + treeEntry = myelist->GetEntryAndTree(el,treenum); + chainEntry = treeEntry+ch->GetTreeOffset()[treenum]; + + if (ch->GetEntry(chainEntry)) { + if (el%1000000==0) { + inte->DrawH(); + c1->Modified(); + c1->Update(); + cout << "+"; + } + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(x+ix, y+iy); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + if (skip==0) { + if (heta==NULL) + inte->addToFlatField(data,etax,etay); + else { + corn=slsInterpolation::calcQuad(data, tot, totquad, sDum); + if (tot>0) { + r=totquad/tot; + // cout << totquad<< " " << tot << endl; + if (r>0.8 && r<1.2) { + slsInterpolation::calcEta(totquad, sDum, etax, etay); + heta->Fill(etax,etay); + slsInterpolation::calcEtaL(totquad, corn, sDum, etaxL, etayL); + if (hetaL) hetaL->Fill(etaxL,etayL); + slsInterpolation::calcEta3(data, etax3, etay3, tot); + if (heta3) heta3->Fill(etax3,etay3); + slsInterpolation::calcEta3X(data, etax3x, etay3x, tot); + if (heta3X) heta3X->Fill(etax3x,etay3x); + + if (h1) h1->Fill(data[4]); + if (h2) h2->Fill(totquad); + if (h3) h3->Fill(tot); + if (hcorn) { + switch(corn) { + case TOP_LEFT: + hcorn->Fill(-0.5,0.5); + break; + case TOP_RIGHT: + hcorn->Fill(0.5,0.5); + break; + case BOTTOM_LEFT: + hcorn->Fill(-0.5,-0.5); + break; + case BOTTOM_RIGHT: + hcorn->Fill(0.5,-0.5); + break; + default: + break; + } + + } + } + } + } + } + } + } + + // inte->DrawH(); + // c1->Modified(); + // c1->Update(); + // cout << "+"; + return inte; +} + + + + +TH2F *makeSpectrum(TChain *ch, Long64_t nentries=-1) { + Long64_t ie; + Double_t dum[9], etax, etay; + Int_t x, y, f; + Double_t sum, totquad, sDum[2][2]; + // if (nentries<0) + // nentries=ch->GetEntries(); + // cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + TH2F *h1=new TH2F("h1","h1",1000,0,5000,NY*NX,0,NY*NX); + + TCanvas *c1=new TCanvas(); + h1->Draw("colz"); + TEntryList *myelist=ch->GetEntryList(); + Long64_t listEntries = myelist->GetN(); + + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + + Long64_t treeEntry, chainEntry; + Int_t treenum; + + cout << "List has " << listEntries*1E-6 << "M entries " << endl; + + for (Long64_t el = 0; el < listEntries; el++) { + treeEntry = myelist->GetEntryAndTree(el,treenum); + chainEntry = treeEntry+ch->GetTreeOffset()[treenum]; + + if (ch->GetEntry(chainEntry)) { + if (el%1000000==0) { + c1->Modified(); + c1->Update(); + cout << "+"; + } + h1->Fill(dum[4],(y-YMIN)*NX+x-XMIN); + } + } + return h1; +} + +TChain *discardDoubles(TChain *ch, int xmin=0, int xmax=400, int ymin=0, int ymax=400, Double_t smin=2000, Double_t smax=3000, Long64_t nentries=-1) { + + TEntryList *elist=new TEntryList(ch); + // TEntryList *elist1=new TEntryList(ch); + elist->SetName("nodoubles"); + if (nentries<0) + nentries=ch->GetEntries(); + Long64_t ie, ie0, iemax, i; + Double_t dum[9], tot; + Int_t x, y, f, ix, iy; + + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + + TH2I *hf=new TH2I("hf","frame",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax); + + Int_t ff=-1, bad=0; + Int_t ibad, igood, iph; + TCanvas *c1=new TCanvas(); + + for (ie=0; ieReset(); + ie0=ie; + iemax=ie0; + iph=0; + while (ff==-1 || f==ff) { + // cout << f << " " << ff << endl; + ch->GetEntry(iemax); + iemax++; + iph++; + if (ff<0) ff=f; + + for (ix=-1; ix<2; ix++) + for (iy=-1; iy<2; iy++) + hf->Fill(x+ix, y+iy); + // if (iph>160000) break; + } + // cout <<"*************************************"<< iemax<< " " << ie0 << " " << ff << " " << f << " " << iph<< endl; + hf->Draw("colz"); + c1->Modified(); + c1->Update(); + + ibad=0; + igood=0; + for (i=ie0; iGetEntry(i); + tot=0; + for (ix=-1; ix<2; ix++) { + for (iy=-1; iy<2; iy++) { + tot+=dum[ix+1+(iy+1)*3]; + if (hf->GetBinContent(x+ix-xmin+1, y+iy-ymin+1)>1.1) { + // cout << "Bad " << x+ix<< " " << y+iy << endl; + bad=1; + } + } + } + if (tot<=smin || tot>=smax) { + bad=1; + // // cout << "Bad Tot "<< tot << endl; + } + if (bad==0) { + elist->Enter(i,ch); + igood++; + } else { + ibad++; + } + + } + if (f%10000==0) + cout << "Frame " << f << " good " << igood << " bad " << ibad << endl; + + ff=-1; + ie=i; + // + } + + // delete hf; + ch->SetEntryList(elist); + return ch; +} + +TChain *discardDoublesQuad(TChain *ch, int xmin=0, int xmax=400, int ymin=0, int ymax=400, Double_t smin=2000, Double_t smax=3000, Long64_t nentries=-1) { + + TEntryList *elist=new TEntryList(ch); + // TEntryList *elist1=new TEntryList(ch); + elist->SetName("nodoubles"); + if (nentries<0) + nentries=ch->GetEntries(); + Long64_t ie, ie0, iemax, i; + Double_t dum[9], tot, totquad, sDum[2][2]; + Int_t x, y, f, ix, iy, corn; + + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + + TH2I *hf=new TH2I("hf","frame",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax); + + Int_t ff=-1, bad=0; + Int_t ibad, igood, iph; + // TCanvas *c1=new TCanvas(); + + for (ie=0; ieReset(); + ie0=ie; + iemax=ie0; + iph=0; + while (ff==-1 || f==ff) { + // cout << f << " " << ff << endl; + ch->GetEntry(iemax); + iemax++; + iph++; + if (ff<0) ff=f; + corn=slsInterpolation::calcQuad(dum, tot, totquad, sDum); + // for (ix=-1; ix<2; ix++) + // for (iy=-1; iy<2; iy++) + // hf->Fill(x+ix, y+iy); + // if (iph>160000) break; + hf->Fill(x,y); + switch(corn) { + case TOP_LEFT: + hf->Fill(x-1,y); + hf->Fill(x-1,y+1); + hf->Fill(x,y+1); + break; + case TOP_RIGHT: + hf->Fill(x+1,y); + hf->Fill(x+1,y+1); + hf->Fill(x,y+1); + break; + case BOTTOM_LEFT: + hf->Fill(x-1,y); + hf->Fill(x-1,y-1); + hf->Fill(x,y-1); + break; + case BOTTOM_RIGHT: + hf->Fill(x+1,y); + hf->Fill(x+1,y-1); + hf->Fill(x,y-1); + break; + default: + break; + } + } + // cout <<"*************************************"<< iemax<< " " << ie0 << " " << ff << " " << f << " " << iph<< endl; + // hf->Draw("colz"); + // c1->Modified(); + // c1->Update(); + + ibad=0; + igood=0; + for (i=ie0; iGetEntry(i); + tot=0; + for (ix=-1; ix<2; ix++) { + for (iy=-1; iy<2; iy++) { + tot+=dum[ix+1+(iy+1)*3]; + if (hf->GetBinContent(x+ix-xmin+1, y+iy-ymin+1)>1.1) { + // cout << "Bad " << x+ix<< " " << y+iy << endl; + bad=1; + } + } + } + if (tot<=smin || tot>=smax) { + bad=1; + // // cout << "Bad Tot "<< tot << endl; + } + if (bad==0) { + elist->Enter(i,ch); + igood++; + } else { + ibad++; + } + + } + if (f%10000==0) + cout << "Frame " << f << " good " << igood << " bad " << ibad << endl; + + ff=-1; + ie=i; + // + } + + delete hf; + ch->SetEntryList(elist); + return ch; +} + + + + + + + + + + + + + + + + + +TH2F *gainMap(TH2F *h1) { + TH1D *hh; + // fspixel->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","Corner"); + Double_t par[7], epar[7]; + TH2F *gmap=new TH2F("gmap","gmap",NX,XMIN,XMIN+NX,NY,YMIN,YMIN+NY); + int ix, iy; + energyCalibration *ecal=new energyCalibration(); + ecal->fixParameter(0,0); + ecal->fixParameter(1,0); + ecal->setScanSign(1); + ecal->setFitRange(2000,3000); + ecal->setPlotFlag(0); + for (int ich=0; ichGetNbinsY(); ich++) { + hh=h1->ProjectionX("hh",ich+1,ich+1); + if (hh->GetEntries()>10000) { + par[0]=0; + par[1]=0; + par[2]=2500; + par[3]=20; + par[4]=1000; + par[5]=0.1; + par[6]-0.1; + ecal->setStartParameters(par); + if (ecal->fitSpectrumPixel(hh,par, epar)) { + if (par[2]>2000 && par[2]<3000) { + ix=ich%NX+XMIN; + iy=ich/NX+YMIN; + gmap->SetBinContent(ix+1-XMIN,iy+1-YMIN,par[2]); + } + } + + } + delete hh; + } + + return gmap; + +} +TChain *m1() { + Long64_t nen=1000000000; //-1; + rootlogon(); + TChain *ch=new TChain("blank"); + + ch->Add("/local_zfs_raid/tomcat_20160528/trees/blank_*.root"); + discardDoublesQuad(ch, 0, 400, 100, 300, 2000, 3000, nen); + TFile *fout=new TFile("/local_zfs_raid/tomcat_20160528/trees/elist.root", "RECREATE"); + TEntryList *elist=ch->GetEntryList(); + elist->Write("nodoubles",TObject::kOverwrite); + fout->Close(); + + return ch; + +} + + +void mm() { + + int nen=1000000; //-1 + + rootlogon(); + // TFile *fout=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root", "RECREATE"); + + + TFile *fL=new TFile("/local_zfs_raid/tomcat_20160528/trees/elist.root"); + TEntryList *elist =(TEntryList*)fL->Get("nodoubles"); + + TFile *fg=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root"); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + + TChain *ch=new TChain("blank"); + + ch->Add("/local_zfs_raid/tomcat_20160528/trees/blank_*.root"); + + // cout << "Chain has " << ch->GetEntries()*1E-6 << "M entries " << endl; + cout << "creating entry list"<< endl; + // ch->Draw(">>elist", "Sum$(data)>2000 && Sum$(data)<3000", "entrylist"); + //TEntryList *elist = (TEntryList*)gDirectory->Get("elist"); + + // TEntryList *elist =(TEntryList*)fg->Get("elist"); + + ch->SetEntryList(elist); + + + + // fout->cd(); + // elist->Write("elist",TObject::kOverwrite); + + // cout << "creating spectrum"<< endl; + //TH2F *h1=makeSpectrum(ch); + //fout->cd(); + //h1->Write("hg",TObject::kOverwrite); + + + //cout << "creating gmap"<< endl; + //TH2F *gmap=gainMap(h1); + //new TCanvas(); + //gmap->Draw("colz"); + //fout->cd(); + //gmap->Write("gmap",TObject::kOverwrite); + + + cout << "creating eta & c."<< endl; + //TFile *fout=new TFile("/local_zfs_raid/tomcat_20160528/trees/eta_dists.root", "RECREATE"); THStack *hs=new THStack(); + // TFile *fg=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_blank.root"); + // TH2F *gmap=(TH2F*)fg->Get("gmap"); + // interpolation_EtaVEL *inte=fillEta(ch,nen, gmap, hs); + // cout << "Writing to file "<< endl; + // fout->cd(); + // inte->setEta((TH2F*)NULL)->Write(); + // inte->setEta((EtaVEL*)NULL)->Write(); + // int nh=-1; + // if (hs->GetHists())nh=hs->GetHists()->GetEntries(); + // for (int ih=0; ihGetHists()->At(ih))->Write(); + // // fout->Close(); + // cout << "Finished "<< endl; +} + +TH2F *fillFF(TChain *ch, interpolation_EtaVEL *eVel, Long64_t nentries=-1, TH2F *gmap=NULL, TH2F *hff=NULL) { + + + cout << "****************** Filling flat field " << nentries << endl; + Double_t dum[9], data[9], gain=1; + Int_t px, py; + Int_t x, y, f, ix, iy; + int skip=0; + Long64_t treeEntry, chainEntry; + Int_t treenum; + if (hff==NULL) + hff=new TH2F("hff","Flat field",25,-12.5,12.5,25,-12.5,12.5); + TAxis *ax=hff->GetXaxis(); + TAxis *ay=hff->GetYaxis(); + + + TEntryList *myelist=ch->GetEntryList(); + Long64_t listEntries = myelist->GetN(); + + cout << "List has " << listEntries*1E-6 << "M entries " << endl; + + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + + if (nentries<=0) nentries=listEntries; + + cout << "Nentries is " << nentries*1E-6 << "M entries " << endl; + + for (Long64_t el = 0; el < nentries; el++) { + // if (el%1000==0) cout << (el*100.)/nentries<< "% " << nentries<< endl; + skip=0; + treeEntry = myelist->GetEntryAndTree(el,treenum); + chainEntry = treeEntry+ch->GetTreeOffset()[treenum]; + + ch->GetEntry(chainEntry); + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(x+ix, y+iy); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + if (skip==0) { + eVel->getInterpolatedBin(data, px, py); + if (px>=0 && py>=0) hff->Fill(ax->GetBinCenter(px+1),ay->GetBinCenter(py+1)); + } + } + return hff; +} + + + + +void pos_dependent_etaX_etaY_normalization() { + + // int ok; + TFile *feta=new TFile("/local_zfs_raid/tomcat_20160528/trees/eta_dists.root"); + TH2F* heta=(TH2F*)feta->Get("heta"); + + + Double_t bsize=1./25.; + + + + cout << "Counts per bin required are "<< heta->Integral()/625. << endl; + cout << "Counts per projection required are "<< heta->Integral()/25. << endl; + + Double_t meanC=heta->Integral()/625.; + + + Double_t bx[26], by[26]; + + TH2I *hhx=new TH2I("hhx","hhx",heta->GetNbinsX(),-0.05, 1.05, heta->GetNbinsY(),-0.05, 1.05); + TH2I *hhy=new TH2I("hhy","hhy",heta->GetNbinsX(),-0.05, 1.05, heta->GetNbinsY(),-0.05, 1.05); + + TH2D *hff=new TH2D("hff","hff",25,0, 25, 25,0, 25); + + Double_t cc[50]; + for (int i=0; i<50; i++) cc[i]=i+0.5; + hhx->SetContour(50,cc); + hhy->SetContour(50,cc); + + + + TH1D *hx; + TH1D *hy; + + TH1D *hix; + TH1D *hiy; + int ii=0; + cout << "X"<GetYaxis()->FindBin(0.02); ibGetYaxis()->FindBin(1.02); ib++) { + for (int ib=0; ibGetNbinsY(); ib++) { + // if (ibGetYaxis()->FindBin(0.02) || ib>heta->GetYaxis()->FindBin(1.02)) + // else + //nb=0; + hx=heta->ProjectionX("px",ib+1,ib+1); + // while (hx->Integral()Integral() << endl; + // nb++; + // hx=heta->ProjectionX("px",ib+1,ib+1+nb); + // if ((ib+nb)>heta->GetNbinsY()) break; + // } + hix=(TH1D*)hx->Clone("hix"); + for (int ibx=0; ibxGetNbinsX(); ibx++) + hix->SetBinContent(ibx+1,hix->GetBinContent(ibx)+hix->GetBinContent(ibx+1)); + hix->Scale(1./hx->Integral()); + ii=1; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + // for (int ibb=0; ibb<=nb; ibb++) + hhx->SetBinContent(ibx+1, ib+1,ii-1); + if (hix->GetBinContent(ibx+1)>=((double)ii)*bsize) { + ii++; + } + } + delete hx; + delete hix; + + // ib=ib+nb; + + } + + cout << "Y"<GetXaxis()->FindBin(0.02); ibGetXaxis()->FindBin(1.02); ib++) { + for (int ib=0; ibGetNbinsX(); ib++) { + // if (ibGetXaxis()->FindBin(0.02) || ib>heta->GetXaxis()->FindBin(1.02)) + // else + // nb=0; + + hx=heta->ProjectionY("py",ib+1,ib+1); + // while (hx->Integral()ProjectionY("px",ib+1,ib+1+nb); + // if ((ib+nb)>heta->GetNbinsX()) break; + // } + hix=(TH1D*)hx->Clone("hix"); + for (int ibx=0; ibxGetNbinsX(); ibx++) + hix->SetBinContent(ibx+1,hix->GetBinContent(ibx)+hix->GetBinContent(ibx+1)); + hix->Scale(1./hx->Integral()); + ii=1; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + // for (int ibb=0; ibb<=nb; ibb++) + hhy->SetBinContent( ib+1, ibx+1,ii-1); + if (hix->GetBinContent(ibx+1)>=((double)ii)*bsize) { + ii++; + } + } + delete hx; + delete hix; + + // ib=ib+nb; + } + + new TCanvas(); + heta->Draw("colz"); + hhy->Draw("cont2 same"); + hhx->Draw("cont2 same"); + + Double_t chi_sq=0; + Double_t v; + cout << "FF"<GetXaxis()->FindBin(0.02); ibxGetXaxis()->FindBin(1.02); ibx++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + // for (int iby=heta->GetYaxis()->FindBin(0.02); ibyGetYaxis()->FindBin(1.02); iby++) { + hff->Fill(hhx->GetBinContent(ibx+1,iby+1)+0.5,hhy->GetBinContent(ibx+1,iby+1)+0.5,heta->GetBinContent(ibx+1,iby+1)); + } + } + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + v=hff->GetBinContent(ibx+1,iby+1); + if (v>0) + chi_sq+=(v-meanC)*(v-meanC)/v; + else + chi_sq+=(v-meanC)*(v-meanC); + } + } + cout <<"*******chi square is " << chi_sq/625. << endl; + + + new TCanvas(); + hff->Draw("colz"); + hff->Scale(1./meanC); + + + TH2F *he2=new TH2F("he2","he2",heta->GetNbinsX(),-0.05, 1.05, heta->GetNbinsY(),-0.05, 1.05); + + int ix, iy; + + // for (int ibx=heta->GetXaxis()->FindBin(0.02); ibxGetXaxis()->FindBin(1.02); ibx++) { + // for (int iby=heta->GetYaxis()->FindBin(0.02); ibyGetYaxis()->FindBin(1.02); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + ix=hhx->GetBinContent(ibx+1,iby+1); + iy=hhy->GetBinContent(ibx+1,iby+1); + he2->SetBinContent(ibx+1,iby+1,hff->GetBinContent(ix+1,iy+1)); + } + } + + new TCanvas(); + he2->Draw("colz"); + +} + + +void global_etaX_etaY_normalization() { + + // int ok; + TFile *feta=new TFile("/local_zfs_raid/tomcat_20160528/trees/eta_dists.root"); + TH2F* heta=(TH2F*)feta->Get("heta"); + + + Double_t bsize=1./25.; + + + + cout << "Counts per bin required are "<< heta->Integral()/625. << endl; + cout << "Counts per projection required are "<< heta->Integral()/25. << endl; + + Double_t meanC=heta->Integral()/625.; + + + Double_t bx[26], by[26]; + + TH2I *hhx=new TH2I("hhx","hhx",heta->GetNbinsX(),-0.05, 1.05, heta->GetNbinsY(),-0.05, 1.05); + TH2I *hhy=new TH2I("hhy","hhy",heta->GetNbinsX(),-0.05, 1.05, heta->GetNbinsY(),-0.05, 1.05); + + TH2D *hff=new TH2D("hff","hff",25,0, 25, 25,0, 25); + + Double_t cc[50]; + for (int i=0; i<50; i++) cc[i]=i+0.5; + hhx->SetContour(50,cc); + hhy->SetContour(50,cc); + + + + TH1D *hx; + TH1D *hy; + + TH1D *hix; + TH1D *hiy; + int ii=0; + cout << "X"<ProjectionX(); + hix=(TH1D*)hx->Clone("hix"); + for (int ibx=0; ibxGetNbinsX(); ibx++) + hix->SetBinContent(ibx+1,hix->GetBinContent(ibx)+hix->GetBinContent(ibx+1)); + hix->Scale(1./hx->Integral()); + // for (int ib=heta->GetYaxis()->FindBin(0.02); ibGetYaxis()->FindBin(1.02); ib++) { + for (int ib=0; ibGetNbinsY(); ib++) { + // if (ibGetYaxis()->FindBin(0.02) || ib>heta->GetYaxis()->FindBin(1.02)) + // else + //nb=0; + // hx=heta->ProjectionX("px",ib+1-0.5*nb,ib+1+0.5*nb); + // while (hx->Integral()Integral() << endl; + // nb++; + // hx=heta->ProjectionX("px",ib+1,ib+1+nb); + // if ((ib+nb)>heta->GetNbinsY()) break; + // } + ii=1; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + // for (int ibb=0; ibb<=nb; ibb++) + hhx->SetBinContent(ibx+1, ib+1,ii-1); + if (hix->GetBinContent(ibx+1)>=((double)ii)*bsize) { + ii++; + } + } + + // ib=ib+nb; + + } + + + delete hx; + delete hix; + + hx=heta->ProjectionY(); + + hix=(TH1D*)hx->Clone("hix"); + for (int ibx=0; ibxGetNbinsX(); ibx++) + hix->SetBinContent(ibx+1,hix->GetBinContent(ibx)+hix->GetBinContent(ibx+1)); + hix->Scale(1./hx->Integral()); + + + + cout << "Y"<GetXaxis()->FindBin(0.02); ibGetXaxis()->FindBin(1.02); ib++) { + for (int ib=0; ibGetNbinsX(); ib++) { + // if (ibGetXaxis()->FindBin(0.02) || ib>heta->GetXaxis()->FindBin(1.02)) + // else + // nb=0; + + ii=1; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + // for (int ibb=0; ibb<=nb; ibb++) + hhy->SetBinContent( ib+1, ibx+1,ii-1); + if (hix->GetBinContent(ibx+1)>=((double)ii)*bsize) { + ii++; + } + } + + // ib=ib+nb; + } + + delete hx; + delete hix; + + new TCanvas(); + heta->Draw("colz"); + hhy->Draw("cont2 same"); + hhx->Draw("cont2 same"); + + Double_t chi_sq=0; + Double_t v; + cout << "FF"<GetXaxis()->FindBin(0.02); ibxGetXaxis()->FindBin(1.02); ibx++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + // for (int iby=heta->GetYaxis()->FindBin(0.02); ibyGetYaxis()->FindBin(1.02); iby++) { + hff->Fill(hhx->GetBinContent(ibx+1,iby+1)+0.5,hhy->GetBinContent(ibx+1,iby+1)+0.5,heta->GetBinContent(ibx+1,iby+1)); + } + } + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + v=hff->GetBinContent(ibx+1,iby+1); + if (v>0) + chi_sq+=(v-meanC)*(v-meanC)/v; + else + chi_sq+=(v-meanC)*(v-meanC); + } + } + cout <<"*******chi square is " << chi_sq/625. << endl; + + + new TCanvas(); + hff->Draw("colz"); + hff->Scale(1./meanC); + + + TH2F *he2=new TH2F("he2","he2",heta->GetNbinsX(),-0.05, 1.05, heta->GetNbinsY(),-0.05, 1.05); + + int ix, iy; + + // for (int ibx=heta->GetXaxis()->FindBin(0.02); ibxGetXaxis()->FindBin(1.02); ibx++) { + // for (int iby=heta->GetYaxis()->FindBin(0.02); ibyGetYaxis()->FindBin(1.02); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + ix=hhx->GetBinContent(ibx+1,iby+1); + iy=hhy->GetBinContent(ibx+1,iby+1); + he2->SetBinContent(ibx+1,iby+1,hff->GetBinContent(ix+1,iy+1)); + } + } + + new TCanvas(); + he2->Draw("colz"); + +} + + + + + + + + + + + + + + + + +void hh() { + rootlogon(); + + + TFile *fL=new TFile("/local_zfs_raid/tomcat_20160528/trees/elist.root"); + TEntryList *elist =(TEntryList*)fL->Get("nodoubles"); + + TFile *fg=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root"); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + + TChain *ch=new TChain("blank"); + + ch->Add("/local_zfs_raid/tomcat_20160528/trees/blank_*.root"); + + + + ch->SetEntryList(elist); + + + + + int ok; + TFile *feta=new TFile("/local_zfs_raid/tomcat_20160528/trees/eta_dists.root"); + // EtaVEL* *eVel=(TH2F*)feta->Get("EtaVEL"); + TH2F* heta=(TH2F*)feta->Get("heta"); + heta->SetName("hheta"); + interpolation_EtaVEL *eVel=new interpolation_EtaVEL(400, 200, 25, -0.05,1.05); + int it=0; + char tit[10000]; + eVel->setEta(heta); + TH2F *hff=new TH2F("hff","Flat field",25,-12.5,12.5,25,-12.5,12.5); + hff->SetTitle("; #pi_{x} (#mum); #pi_{y} (#mum); Counts (A.U.)"); + hff->SetStats(kFALSE); + hff->SetMinimum(0); + // hff->SetMaximum(20E3); + hff->GetZaxis()->SetTitleOffset(1.3); + hff->GetZaxis()->SetLabelSize(0.02); + + // TH2F *hgrid=new TH2F("hgrid","",25,-12.5,12.5,25,-12.5,12.5); + // hgrid->SetTitle("; #pi_{x} (#mum); #pi_{y} (#mum); Counts (A.U.)"); + // hgrid->GetXaxis()->SetNdivisions(25); + // hgrid->GetYaxis()->SetNdivisions(25); + // hgrid->GetYaxis()->SetLabelOffset(999.); + // hgrid->GetXaxis()->SetLabelOffset(999.); + TGraph *g_chi=new TGraph(); + g_chi->SetMarkerStyle(20); + TCanvas *c1=new TCanvas("c1","c1",1000,500); + TCanvas *c2=new TCanvas("c2","c2",500,500); + c2->SetLogy(kTRUE); + c1->Divide(2,1); + + c1->cd(1); + TPad* pad = (TPad*)c1->GetPad(1); + //pad->SetLeftMargin(0.15); + // pad->SetBottomMargin(0.15); + // pad->SetTopMargin(0.2); + // pad->SetRightMargin(0.2); + // pad->SetFillColor(0); + pad->SetBorderMode(0); + pad->SetLogz(kTRUE); + + c1->cd(2); + pad = (TPad*)c1->GetPad(2); + pad->SetLeftMargin(0.1); + // pad->SetBottomMargin(0.15); + // pad->SetTopMargin(0.2); + pad->SetRightMargin(0.15); + // pad->SetFillColor(0); + pad->SetBorderMode(0); + // pad->SetLogz(kTRUE); + pad->SetGridx(kTRUE); + pad->SetGridy(kTRUE); + + + + + heta->SetStats(kFALSE); + // c1->cd(2); + // hff->Fill(0.,0.); + // hff->Draw("colz"); + // TPaletteAxis *palette1 =(TPaletteAxis*)hff->GetListOfFunctions()->FindObject("palette"); + // palette1->SetLabelSize(0.02); + + + while (ok==0 && it<10000) { + sprintf(tit,"Iteration %d; #eta_{x}; #eta_{y}; Counts (A.U.)",it); + hff->Reset(); + fillFF(ch, eVel, 100000, gmap, hff); + + heta->SetTitle(tit); + c1->cd(1); + eVel->DrawH(); + + c1->cd(2); + // hgrid->Draw(); + hff->Draw("colz"); + + c1->cd(); + c1->Modified(); + c1->Update(); + c1->SaveAs("test1.gif+"); + + + eVel->prepareInterpolation(ok,1); + cout << "**************************** Iteration number " << it << endl; + if (eVel->getChiSq()<1E10) + g_chi->SetPoint(it,it,eVel->getChiSq()); + else + g_chi->SetPoint(it,it,1E10); + it++; + if (g_chi->GetN()>0) { + c2->cd(); + g_chi->Draw("ALP"); + c2->Modified(); + c2->Update(); + } + } + + sprintf(tit,"Iteration %d",it); + hff->Reset(); + fillFF(ch, eVel, 100000, gmap, hff); + + heta->SetTitle(tit); + + c1->cd(1); + eVel->DrawH(); + + c1->cd(2); + // hgrid->Draw(); + hff->Draw("colz"); + + c1->Modified(); + c1->Update(); + c1->SaveAs("test1.gif+"); + + +} + +TH2D *interpGrating1D() { + + + + TFile *fg=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root"); + TH2F *gmap=NULL;//(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + TChain *ch=new TChain("grating_1d"); + ch->Add("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t*_v4.root"); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + if (iby==0) + hietaY->SetBinContent(ibx+1, iby+1,heta->GetBinContent(ibx+1, iby+1)/hetaX->GetBinContent(ibx+1)); + else + hietaY->SetBinContent(ibx+1, iby+1,hietaY->GetBinContent(ibx+1, iby)+heta->GetBinContent(ibx+1, iby+1)/hetaX->GetBinContent(ibx+1)); + + if (ibx==0) + hietaX->SetBinContent(ibx+1, iby+1,heta->GetBinContent(ibx+1, iby+1)/hetaY->GetBinContent(iby+1)); + else + hietaX->SetBinContent(ibx+1, iby+1,hietaX->GetBinContent(ibx, iby+1)+heta->GetBinContent(ibx+1, iby+1)/hetaY->GetBinContent(iby+1)); + + if (ibx==0) { + if (iby==0) + h1etaY->SetBinContent(iby+1,hetaY->GetBinContent(iby+1)/hetaY->Integral()); + else + h1etaY->SetBinContent(iby+1,h1etaY->GetBinContent(iby)+hetaY->GetBinContent(iby+1)/hetaY->Integral()); + + } + } + if (ibx==0) + h1etaX->SetBinContent(ibx+1,hetaX->GetBinContent(ibx+1)/hetaX->Integral()); + else + h1etaX->SetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral()); + + } + + + TFile *fout=new TFile("/local_zfs_raid/tomcat_20160528/trees/img_grating_1d_eta.root", "RECREATE"); + + + + new TCanvas(); + hietaY->Draw("colz"); + gPad->Modified(); + gPad->Update(); + + new TCanvas(); + hietaX->Draw("colz"); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaX->Draw(); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaY->Draw(); + gPad->Modified(); + gPad->Update(); + + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*25,0,400,200*25,100,300); + TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + TCanvas *c1=new TCanvas(); + imgLR->Draw("colz"); + TCanvas *c2=new TCanvas(); + imgHR->Draw("colz"); + TCanvas *c3=new TCanvas(); + imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + Double_t etamax=hietaX->GetXaxis()->GetBinCenter(hietaX->GetNbinsX()); + Double_t etamin=hietaX->GetXaxis()->GetBinCenter(0); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + //slsInterpolation::calcEta(dum, etax, etay, sum, totquad, sDum); + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (etaxetamax || etay>etamax || totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + imgHR->Fill(x-0.5+hietaX->GetBinContent(hietaX->GetXaxis()->FindBin(etax),hietaX->GetYaxis()->FindBin(etay)),y-0.5+hietaY->GetBinContent(hietaY->GetXaxis()->FindBin(etax),hietaY->GetYaxis()->FindBin(etay))); + imgHR1->Fill(x-0.5+h1etaX->GetBinContent(h1etaX->FindBin(etax)), y-0.5+h1etaY->GetBinContent(h1etaY->FindBin(etay))); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified();export PATH=/afs/psi.ch/project/sls_det_firmware/jungfrau_software/uClinux-2010_64bit/bfin-uclinux/bin:$PATH + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + } + + } + imgLR->Write("imgLR"); + imgHR->Write("imgHR"); + imgHR1->Write("imgHR1"); + + c1->cd(); + imgLR->DrawCopy("colz"); + c2->cd(); + imgHR->DrawCopy("colz"); + c3->cd(); + imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + + + + + + + + + +TH2D *interpBlank() { + + + const Double_t Delta_eta=0.025; + + + TFile *fg=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root"); + TH2F *gmap=NULL;//(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + TChain *ch=new TChain("blank"); + ch->Add("/local_zfs_raid/tomcat_20160528/trees/blank_t*.root"); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + + int binminX=heta->GetXaxis()->FindBin(-Delta_eta); + int binmaxX=heta->GetXaxis()->FindBin(1+Delta_eta); + int binminY=heta->GetYaxis()->FindBin(-Delta_eta); + int binmaxY=heta->GetYaxis()->FindBin(1+Delta_eta); + + + for (int ibx=binminX; ibxSetBinContent(ibx, iby,heta->GetBinContent(ibx, iby)/hetaX->GetBinContent(ibx)); + else + hietaY->SetBinContent(ibx, iby,hietaY->GetBinContent(ibx, iby-1)+heta->GetBinContent(ibx, iby)/hetaX->GetBinContent(ibx)); + + if (ibx==0) + hietaX->SetBinContent(ibx, iby,heta->GetBinContent(ibx, iby)/hetaY->GetBinContent(iby)); + else + hietaX->SetBinContent(ibx, iby,hietaX->GetBinContent(ibx-1, iby)+heta->GetBinContent(ibx, iby)/hetaY->GetBinContent(iby)); + + if (ibx==binminX) { + if (iby==0) + h1etaY->SetBinContent(iby+1,hetaY->GetBinContent(iby)/hetaY->Integral(binminY,binmaxY)); + else + h1etaY->SetBinContent(iby+1,h1etaY->GetBinContent(iby-1)+hetaY->GetBinContent(iby)/hetaY->Integral(binminY,binmaxY)); + + } + } + if (ibx==0) + h1etaX->SetBinContent(ibx+1,hetaX->GetBinContent(ibx+1)/hetaX->Integral(binminX,binmaxX)); + else + h1etaX->SetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral(binminX,binmaxX)); + + } + + + TFile *fout=new TFile("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta.root", "RECREATE"); + + + + new TCanvas(); + hietaY->Draw("colz"); + gPad->Modified(); + gPad->Update(); + + new TCanvas(); + hietaX->Draw("colz"); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaX->Draw(); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaY->Draw(); + gPad->Modified(); + gPad->Update(); + + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*25,0,400,200*25,100,300); + TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + TCanvas *c1=new TCanvas(); + imgLR->Draw("colz"); + TCanvas *c2=new TCanvas(); + imgHR->Draw("colz"); + TCanvas *c3=new TCanvas(); + imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + Double_t etamax=hietaX->GetXaxis()->GetBinCenter(hietaX->GetNbinsX()); + Double_t etamin=hietaX->GetXaxis()->GetBinCenter(0); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (etaxetamax || etay>etamax || totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + imgHR->Fill(x-0.5+hietaX->GetBinContent(hietaX->GetXaxis()->FindBin(etax),hietaX->GetYaxis()->FindBin(etay)),y-0.5+hietaY->GetBinContent(hietaY->GetXaxis()->FindBin(etax),hietaY->GetYaxis()->FindBin(etay))); + imgHR1->Fill(x-0.5+h1etaX->GetBinContent(h1etaX->FindBin(etax)), y-0.5+h1etaY->GetBinContent(h1etaY->FindBin(etay))); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified();export PATH=/afs/psi.ch/project/sls_det_firmware/jungfrau_software/uClinux-2010_64bit/bfin-uclinux/bin:$PATH + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + imgLR->Write("imgLR",TObject::kOverwrite); + imgHR->Write("imgHR",TObject::kOverwrite); + imgHR1->Write("imgHR1",TObject::kOverwrite); + } + + } + imgLR->Write("imgLR",TObject::kOverwrite); + imgHR->Write("imgHR",TObject::kOverwrite); + imgHR1->Write("imgHR1",TObject::kOverwrite); + + // c1->cd(); + // imgLR->DrawCopy("colz"); + // c2->cd(); + // imgHR->DrawCopy("colz"); + // c3->cd(); + // imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + + + + + + + +TH2D *interpSampleGrating() { + + + + TFile *fg=new TFile("/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root"); + TH2F *gmap=NULL;//(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + TChain *ch=new TChain("sample_grating_1d"); + ch->Add("/local_zfs_raid/tomcat_20160528/trees/sample_grating_1d_t*.root"); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + if (iby==0) + hietaY->SetBinContent(ibx+1, iby+1,heta->GetBinContent(ibx+1, iby+1)/hetaX->GetBinContent(ibx+1)); + else + hietaY->SetBinContent(ibx+1, iby+1,hietaY->GetBinContent(ibx+1, iby)+heta->GetBinContent(ibx+1, iby+1)/hetaX->GetBinContent(ibx+1)); + + if (ibx==0) + hietaX->SetBinContent(ibx+1, iby+1,heta->GetBinContent(ibx+1, iby+1)/hetaY->GetBinContent(iby+1)); + else + hietaX->SetBinContent(ibx+1, iby+1,hietaX->GetBinContent(ibx, iby+1)+heta->GetBinContent(ibx+1, iby+1)/hetaY->GetBinContent(iby+1)); + + if (ibx==0) { + if (iby==0) + h1etaY->SetBinContent(iby+1,hetaY->GetBinContent(iby+1)/hetaY->Integral()); + else + h1etaY->SetBinContent(iby+1,h1etaY->GetBinContent(iby)+hetaY->GetBinContent(iby+1)/hetaY->Integral()); + + } + } + if (ibx==0) + h1etaX->SetBinContent(ibx+1,hetaX->GetBinContent(ibx+1)/hetaX->Integral()); + else + h1etaX->SetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral()); + + } + + + TFile *fout=new TFile("/local_zfs_raid/tomcat_20160528/trees/img_sample_grating_1d_eta.root", "RECREATE"); + + + + new TCanvas(); + hietaY->Draw("colz"); + gPad->Modified(); + gPad->Update(); + + new TCanvas(); + hietaX->Draw("colz"); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaX->Draw(); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaY->Draw(); + gPad->Modified(); + gPad->Update(); + + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*25,0,400,200*25,100,300); + TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + TCanvas *c1=new TCanvas(); + imgLR->Draw("colz"); + TCanvas *c2=new TCanvas(); + imgHR->Draw("colz"); + TCanvas *c3=new TCanvas(); + imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + Double_t etamax=hietaX->GetXaxis()->GetBinCenter(hietaX->GetNbinsX()); + Double_t etamin=hietaX->GetXaxis()->GetBinCenter(0); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (etaxetamax || etay>etamax || totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + imgHR->Fill(x-0.5+hietaX->GetBinContent(hietaX->GetXaxis()->FindBin(etax),hietaX->GetYaxis()->FindBin(etay)),y-0.5+hietaY->GetBinContent(hietaY->GetXaxis()->FindBin(etax),hietaY->GetYaxis()->FindBin(etay))); + imgHR1->Fill(x-0.5+h1etaX->GetBinContent(h1etaX->FindBin(etax)), y-0.5+h1etaY->GetBinContent(h1etaY->FindBin(etay))); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified();export PATH=/afs/psi.ch/project/sls_det_firmware/jungfrau_software/uClinux-2010_64bit/bfin-uclinux/bin:$PATH + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + imgLR->Write("imgLR",TObject::kOverwrite); + imgHR->Write("imgHR",TObject::kOverwrite); + imgHR1->Write("imgHR1",TObject::kOverwrite); + } + + } + imgLR->Write("imgLR",TObject::kOverwrite); + imgHR->Write("imgHR",TObject::kOverwrite); + imgHR1->Write("imgHR1",TObject::kOverwrite); + + // c1->cd(); + // imgLR->DrawCopy("colz"); + // c2->cd(); + // imgHR->DrawCopy("colz"); + // c3->cd(); + // imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + +TH2D *interp(char *name, int ip, int nb=25) { + + char chainName[1000]; + sprintf(chainName,"/local_zfs_raid/tomcat_20160528/trees/%s_t*.root",name); + + const Double_t Delta_eta=0.05; + + + char gName[1000]; + //sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank_%d.root",ip); + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_gcorr_%d.root",ip); + TFile *fg=new TFile(gName); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + TChain *ch=new TChain(name); + ch->Add(chainName); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + hietaX->Reset(); + hietaY->Reset(); + + + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + h1etaX->Reset(); + h1etaY->Reset(); + + int binminX=heta->GetXaxis()->FindBin(-Delta_eta); + int binmaxX=heta->GetXaxis()->FindBin(1+Delta_eta); + int binminY=heta->GetYaxis()->FindBin(-Delta_eta); + int binmaxY=heta->GetYaxis()->FindBin(1+Delta_eta); + + + TH1D *hetaYp=heta->ProjectionX("py",binminY,binmaxY); + TH1D *hetaXp=heta->ProjectionY("px",binminX,binmaxX); + + for (int ibx=binminX; ibxSetBinContent(ibx, iby,hietaY->GetBinContent(ibx, iby-1)+heta->GetBinContent(ibx, iby)); + + hietaX->SetBinContent(ibx, iby,hietaX->GetBinContent(ibx-1, iby)+heta->GetBinContent(ibx, iby)); + + + + } + + } + + + for (int ibx=binminX; ibxSetBinContent(ibx,iby,hietaY->GetBinContent(ibx,iby)/hietaY->GetBinContent(ibx,binmaxY-1)); + hietaX->SetBinContent(ibx,iby,hietaX->GetBinContent(ibx,iby)/hietaX->GetBinContent(binmaxX-1,iby)); + + } + + } + + for (int iby=binminY; ibySetBinContent(iby+1,h1etaY->GetBinContent(iby)+hetaY->GetBinContent(iby+1)/hetaY->Integral(binminY,binmaxY)); + + for (int ibx=binminX; ibxSetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral(binminX,binmaxX)); + + for (int ibx=binmaxX; ibxGetNbinsX()+1; ibx++) + h1etaX->SetBinContent(ibx,1); + for (int iby=binmaxX; ibyGetNbinsY()+1; iby++) + h1etaY->SetBinContent(iby,1); + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaX->SetBinContent(ibx+1,iby+1,1); + } + } + + + for (int iby=0; ibyGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binminY+1)); + } + } + + + for (int iby=binmaxY-1; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binmaxY-1)); + } + } + + + + for (int iby=binmaxX-1; ibyGetNbinsY(); iby++) { + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaY->SetBinContent(ibx+1,iby+1,1); + } + } + + + + for (int ibx=0; ibxGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binminX+1, iby+1)); + } + } + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binmaxX-1, iby+1)); + } + } + + + + + // new TCanvas(); + // hietaY->Draw("colz"); + // gPad->Modified(); + // gPad->Update(); + + // new TCanvas(); + // hietaX->Draw("colz"); + // gPad->Modified(); + // gPad->Update(); + // new TCanvas(); + // h1etaX->Draw(); + // gPad->Modified(); + // gPad->Update(); + // new TCanvas(); + // h1etaY->Draw(); + // gPad->Modified(); + // gPad->Update(); + + + char foutName[1000]; + sprintf(foutName,"/local_zfs_raid/tomcat_20160528/trees/img_%s_eta_gcorr_nb%d.root",name,nb); + + + TFile *fout=new TFile(foutName, "RECREATE"); + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*nb,0,400,200*nb,100,300); + // TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + // TCanvas *c1=new TCanvas(); + // imgLR->Draw("colz"); + // TCanvas *c2=new TCanvas(); + // imgHR->Draw("colz"); + // TCanvas *c3=new TCanvas(); + // imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + Double_t etamax=hietaX->GetXaxis()->GetBinCenter(hietaX->GetNbinsX()); + Double_t etamin=hietaX->GetXaxis()->GetBinCenter(0); + + char name1[100],name2[100],name3[100]; + sprintf(name1,"%sLR",name); + sprintf(name2,"%sHR",name); + sprintf(name3,"%sHR1",name); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (etaxetamax || etay>etamax || totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + imgHR->Fill(x-0.5+hietaX->GetBinContent(hietaX->GetXaxis()->FindBin(etax),hietaX->GetYaxis()->FindBin(etay)),y-0.5+hietaY->GetBinContent(hietaY->GetXaxis()->FindBin(etax),hietaY->GetYaxis()->FindBin(etay))); + // imgHR1->Fill(x-0.5+h1etaX->GetBinContent(h1etaX->FindBin(etax)), y-0.5+h1etaY->GetBinContent(h1etaY->FindBin(etay))); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified();export PATH=/afs/psi.ch/project/sls_det_firmware/jungfrau_software/uClinux-2010_64bit/bfin-uclinux/bin:$PATH + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + // imgHR1->Write(name3,TObject::kOverwrite); + } + + } + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + // imgHR1->Write(name3,TObject::kOverwrite); + + // c1->cd(); + // imgLR->DrawCopy("colz"); + // c2->cd(); + // imgHR->DrawCopy("colz"); + // c3->cd(); + // imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + +TH2D *interpPos(char *name, int ip) { + + char chainName[1000]; + sprintf(chainName,"/local_zfs_raid/tomcat_20160528/trees/%s_t*.root",name); + + const Double_t Delta_eta=0.015; + + + char gName[1000]; + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank_%d.root",ip); + TFile *fg=new TFile(gName); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + // TH2F *heta=(TH2F*)fg->Get("heta"); + + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/hpos_anna_%d.root",ip); + TFile *fp=new TFile(gName); + TH2F *hpos=(TH2F*)fp->Get("hpos"); + + TChain *ch=new TChain(name); + ch->Add(chainName); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + char foutName[1000]; + sprintf(foutName,"/local_zfs_raid/tomcat_20160528/trees/img_%s_pos_gmap_v2.root",name); + + + TFile *fout=new TFile(foutName, "RECREATE"); + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*25,0,400,200*25,100,300); + TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + TCanvas *c1=new TCanvas(); + imgLR->Draw("colz"); + TCanvas *c2=new TCanvas(); + imgHR->Draw("colz"); + TCanvas *c3=new TCanvas(); + imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + int ppos, xpos, ypos; + char name1[100],name2[100],name3[100]; + sprintf(name1,"%sLR",name); + sprintf(name2,"%sHR",name); + sprintf(name3,"%sHR1",name); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + ppos=hpos->GetBinContent(hpos->GetXaxis()->FindBin(etax),hpos->GetYaxis()->FindBin(etay)); + xpos=ppos/25; + ypos=ppos%25; + imgHR->Fill(x-0.5+((Double_t)xpos+0.1)/25.,y-0.5+((Double_t)ypos+0.1)/25.); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified();export PATH=/afs/psi.ch/project/sls_det_firmware/jungfrau_software/uClinux-2010_64bit/bfin-uclinux/bin:$PATH + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + imgHR1->Write(name3,TObject::kOverwrite); + } + + } + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + imgHR1->Write(name3,TObject::kOverwrite); + + // c1->cd(); + // imgLR->DrawCopy("colz"); + // c2->cd(); + // imgHR->DrawCopy("colz"); + // c3->cd(); + // imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + + + + + +TH2F *iterate(TH2F *heta, TH2F *hpos, int i) { + int ipx,ipy; + + + TH2F *hpos2=(TH2F*)heta->Clone("hpos"); + hpos2->Reset(); + + if (i%2==0) { + + TH1F *hx=new TH1F("hx","hx",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax()); + TH1F *hix=(TH1F*)hx->Clone("hix"); + for ( ipy=0; ipy<25; ipy++) { + hx->Reset(); + hix->Reset(); + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + if ((((int)(hpos->GetBinContent(ibx+1,iby+1)))%25)==ipy) { + hx->SetBinContent(ibx+1,hx->GetBinContent(ibx+1)+heta->GetBinContent(ibx+1,iby+1)); + } + } + } + ipx=0; + for (int ibx=1; ibxGetNbinsX(); ibx++) { + hix->SetBinContent(ibx+1, hix->GetBinContent(ibx)+hx->GetBinContent(ibx+1)); + } + // new TCanvas(); + //hix->Draw(); + + hix->Scale(1./hix->GetBinContent(hix->GetNbinsX())); + for (int ibx=0; ibxGetNbinsX(); ibx++) { + if (hix->GetBinContent(ibx+1)>((float)(ipx+1))/25.) ipx++; + for (int iby=0; ibyGetNbinsY(); iby++) { + if ((((int)hpos->GetBinContent(ibx+1,iby+1))%25)==ipy) { + hpos2->SetBinContent(ibx+1,iby+1,25*ipx+ipy); + + } + } + } + } + + delete hx; + delete hix; + + + } else { + TH1F *hy=new TH1F("hy","hy",heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax()); + TH1F *hiy=(TH1F*)hy->Clone("hiy"); + for ( ipx=0; ipx<25; ipx++) { + hy->Reset(); + hiy->Reset(); + for (int iby=0; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + if ((((int)hpos->GetBinContent(ibx+1,iby+1))/25)==ipx) { + hy->SetBinContent(iby+1,hy->GetBinContent(iby+1)+heta->GetBinContent(ibx+1,iby+1)); + } + } + } + ipy=0; + for (int iby=1; ibyGetNbinsY(); iby++) { + hiy->SetBinContent(iby+1, hiy->GetBinContent(iby)+hy->GetBinContent(iby+1)); + } + hiy->Scale(1./hiy->GetBinContent(hiy->GetNbinsX())); + for (int iby=0; ibyGetNbinsY(); iby++) { + if (hiy->GetBinContent(iby+1)>((float)(ipy+1))/25.) ipy++; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + if ((((int)hpos->GetBinContent(ibx+1,iby+1))/25)==ipx) { + hpos2->SetBinContent(ibx+1,iby+1,ipx*25+ipy);//hpos2->GetBinContent(ibx+1,iby+1)+ipy); + } + } + } + } + delete hy; + delete hiy; + } + + + + // new TCanvas(); + // heta->Draw("colz"); + // hnet->Draw("same"); + + + // new TCanvas(); + // hpos2->Draw("colz"); + + + + + + + + + + + + return hpos2; + + } + + + +TH2F *etaTest() { + + TMultiGraph *mg=new TMultiGraph(); + TGraph *g=new TGraph(); + TGraph *g1=new TGraph(); + + mg->Add(g); + mg->Add(g1); + + const Double_t Delta_eta=0.1; + + + char gName[1000]; + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank.root"); + TFile *fg=new TFile(gName); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + + + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + hietaX->Reset(); + hietaY->Reset(); + + TH2F *hnet=(TH2F*)heta->Clone("hnet"); + hnet->Reset(); + TH2F *hpos=(TH2F*)heta->Clone("hpos"); + hpos->Reset(); + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + h1etaX->Reset(); + h1etaY->Reset(); + + int binminX=heta->GetXaxis()->FindBin(-Delta_eta); + int binmaxX=heta->GetXaxis()->FindBin(1+Delta_eta); + int binminY=heta->GetYaxis()->FindBin(-Delta_eta); + int binmaxY=heta->GetYaxis()->FindBin(1+Delta_eta); + + + TH1D *hetaYp=heta->ProjectionX("py",binminY,binmaxY); + TH1D *hetaXp=heta->ProjectionY("px",binminX,binmaxX); + + for (int ibx=binminX; ibxSetBinContent(ibx, iby,hietaY->GetBinContent(ibx, iby-1)+heta->GetBinContent(ibx, iby)); + + hietaX->SetBinContent(ibx, iby,hietaX->GetBinContent(ibx-1, iby)+heta->GetBinContent(ibx, iby)); + + + + } + + } + + + for (int ibx=binminX; ibxSetBinContent(ibx,iby,hietaY->GetBinContent(ibx,iby)/hietaY->GetBinContent(ibx,binmaxY-1)); + hietaX->SetBinContent(ibx,iby,hietaX->GetBinContent(ibx,iby)/hietaX->GetBinContent(binmaxX-1,iby)); + + } + + } + + for (int iby=binminY; ibySetBinContent(iby+1,h1etaY->GetBinContent(iby)+hetaY->GetBinContent(iby+1)/hetaY->Integral(binminY,binmaxY)); + + for (int ibx=binminX; ibxSetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral(binminX,binmaxX)); + + for (int ibx=binmaxX; ibxGetNbinsX()+1; ibx++) + h1etaX->SetBinContent(ibx,1); + for (int iby=binmaxX; ibyGetNbinsY()+1; iby++) + h1etaY->SetBinContent(iby,1); + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaX->SetBinContent(ibx+1,iby+1,1); + } + } + + + for (int iby=0; ibyGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binminY+1)); + } + } + + + for (int iby=binmaxY-1; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binmaxY-1)); + } + } + + + + for (int iby=binmaxX-1; ibyGetNbinsY(); iby++) { + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaY->SetBinContent(ibx+1,iby+1,1); + } + } + + + + for (int ibx=0; ibxGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binminX+1, iby+1)); + } + } + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binmaxX-1, iby+1)); + } + } + + + int ipx=0; + int ipy=0; + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + ipy=0; + for (int iby=0; ibyGetNbinsY(); iby++) { + if (hietaY->GetBinContent(ibx+1,iby+1)>((float)ipy+1)/25.) { + hnet->SetBinContent(ibx+1,iby+1,1); + ipy++; + } + hpos->SetBinContent(ibx+1,iby+1,ipy); + } + } + + for (int iby=0; ibyGetNbinsY(); iby++) { + ipx=0; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + if (hietaX->GetBinContent(ibx+1,iby+1)>((float)ipx+1)/25.) ipx++; + hpos->SetBinContent(ibx+1,iby+1,hpos->GetBinContent(ibx+1,iby+1)+25*ipx); + } + } + new TCanvas(); + heta->Draw("colz"); + + + for (int iby=1; ibyGetNbinsY(); iby++) { + for (int ibx=1; ibxGetNbinsX(); ibx++) { + if (hpos->GetBinContent(ibx+1,iby+1)!=hpos->GetBinContent(ibx,iby+1)) hnet->SetBinContent(ibx+1,iby+1,1); + if (hpos->GetBinContent(ibx+1,iby+1)!=hpos->GetBinContent(ibx+1,iby)) hnet->SetBinContent(ibx+1,iby+1,1); + else hnet->SetBinContent(ibx+1,iby+1,0); + } + } + + + hnet->DrawCopy("same"); + + + TH2F *hff=new TH2F("hff","hff",25,0,25,25,0,25); + Double_t mm=heta->Integral()/(hff->GetNbinsX()*hff->GetNbinsY()); + + int ibp; + for (int iby=0; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + ibp=hpos->GetBinContent(ibx+1,iby+1); + hff->Fill(ibp/25,ibp%25,heta->GetBinContent(ibx+1,iby+1)); + } + } + + Double_t chi2=0; + Double_t maxS=0; + for (int ibx=0; ibx<25; ibx++) + for (int iby=0; iby<25; iby++) { + chi2+=(hff->GetBinContent(ibx+1,iby+1)-mm)*(hff->GetBinContent(ibx+1,iby+1)-mm)/mm; + if ((hff->GetBinContent(ibx+1,iby+1)-mm)*(hff->GetBinContent(ibx+1,iby+1))>maxS) maxS=(hff->GetBinContent(ibx+1,iby+1)-mm)*(hff->GetBinContent(ibx+1,iby+1)-mm); + } + chi2/=625.; + maxS=TMath::Sqrt(maxS/mm); + + cout << "Start Chi2 "<< chi2 << " maxS "<< maxS << endl; + + g->SetPoint(0,0,chi2); + g1->SetPoint(0,0,maxS); + + new TCanvas(); + char tit[1000]; + sprintf(tit,"Start Chi2=%d",(int)chi2); + hff->SetTitle(tit); + hff->DrawCopy("colz"); + + + TH2F *hpos2; + int i; + new TCanvas(); + for (i=0; i<500; i++) { + hpos2=iterate(heta,hpos,i); + delete hpos; + hpos=hpos2; + + + hff->Reset(); + + for (int iby=0; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + ibp=hpos->GetBinContent(ibx+1,iby+1); + hff->Fill(ibp/25,ibp%25,heta->GetBinContent(ibx+1,iby+1)); + } + } + + chi2=0; + maxS=0; + for (int ibx=0; ibx<25; ibx++) + for (int iby=0; iby<25; iby++) { + chi2+=(hff->GetBinContent(ibx+1,iby+1)-mm)*(hff->GetBinContent(ibx+1,iby+1)-mm)/mm; + if ((hff->GetBinContent(ibx+1,iby+1)-mm)*(hff->GetBinContent(ibx+1,iby+1))>maxS) maxS=(hff->GetBinContent(ibx+1,iby+1)-mm)*(hff->GetBinContent(ibx+1,iby+1)-mm); + + } + chi2/=625.; + maxS=TMath::Sqrt(maxS/mm); + + cout << "Iteration " << i << " Chi2 "<< chi2 << " maxS "<< maxS<< endl; + + g->SetPoint(i+1,i+1,chi2); + g1->SetPoint(i+1,i+1,maxS); + if (i%100==0) { + + g->Draw("ALP"); + g1->Draw("LP"); + gPad->Modified(); + gPad->Update(); + } + + } + + new TCanvas(); + heta->Draw("colz"); + + + for (int iby=1; ibyGetNbinsY(); iby++) { + for (int ibx=1; ibxGetNbinsX(); ibx++) { + if (hpos->GetBinContent(ibx+1,iby+1)!=hpos->GetBinContent(ibx,iby+1)) hnet->SetBinContent(ibx+1,iby+1,1); + else if (hpos->GetBinContent(ibx+1,iby+1)!=hpos->GetBinContent(ibx+1,iby)) hnet->SetBinContent(ibx+1,iby+1,1); + else hnet->SetBinContent(ibx+1,iby+1,0); + } + } + + + hnet->DrawCopy("same"); + + new TCanvas(); + sprintf(tit,"Iteration %d Chi2=%d",i,(int)chi2); + hff->SetTitle(tit); + hff->DrawCopy("colz"); + + TFile *fffff=new TFile("/local_zfs_raid/tomcat_20160528/trees/hpos_anna.root","RECREATE"); + hpos->Write("hpos"); + fffff->Close(); + + + + return hpos; + + +} + + +TH2D *interpV2(char *name, int ip, double sc=30./25.) { + + char chainName[1000]; + sprintf(chainName,"/local_zfs_raid/tomcat_20160528/trees/%s_t*.root",name); + + const Double_t Delta_eta=0.05; + + + char gName[1000]; + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank_%d.root",ip); + TFile *fg=new TFile(gName); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + TChain *ch=new TChain(name); + ch->Add(chainName); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + hietaX->Reset(); + hietaY->Reset(); + + + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + h1etaX->Reset(); + h1etaY->Reset(); + + int binminX=heta->GetXaxis()->FindBin(-Delta_eta); + int binmaxX=heta->GetXaxis()->FindBin(1+Delta_eta); + int binminY=heta->GetYaxis()->FindBin(-Delta_eta); + int binmaxY=heta->GetYaxis()->FindBin(1+Delta_eta); + + + TH1D *hetaYp=heta->ProjectionX("py",binminY,binmaxY); + TH1D *hetaXp=heta->ProjectionY("px",binminX,binmaxX); + + for (int ibx=binminX; ibxSetBinContent(ibx, iby,hietaY->GetBinContent(ibx, iby-1)+heta->GetBinContent(ibx, iby)); + + hietaX->SetBinContent(ibx, iby,hietaX->GetBinContent(ibx-1, iby)+heta->GetBinContent(ibx, iby)); + + + + } + + } + + + for (int ibx=binminX; ibxSetBinContent(ibx,iby,hietaY->GetBinContent(ibx,iby)/hietaY->GetBinContent(ibx,binmaxY-1)); + hietaX->SetBinContent(ibx,iby,hietaX->GetBinContent(ibx,iby)/hietaX->GetBinContent(binmaxX-1,iby)); + + } + + } + + for (int iby=binminY; ibySetBinContent(iby+1,h1etaY->GetBinContent(iby)+hetaY->GetBinContent(iby+1)/hetaY->Integral(binminY,binmaxY)); + + for (int ibx=binminX; ibxSetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral(binminX,binmaxX)); + + for (int ibx=binmaxX; ibxGetNbinsX()+1; ibx++) + h1etaX->SetBinContent(ibx,1); + for (int iby=binmaxX; ibyGetNbinsY()+1; iby++) + h1etaY->SetBinContent(iby,1); + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaX->SetBinContent(ibx+1,iby+1,1); + } + } + + + for (int iby=0; ibyGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binminY+1)); + } + } + + + for (int iby=binmaxY-1; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binmaxY-1)); + } + } + + + + for (int iby=binmaxX-1; ibyGetNbinsY(); iby++) { + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaY->SetBinContent(ibx+1,iby+1,1); + } + } + + + + for (int ibx=0; ibxGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binminX+1, iby+1)); + } + } + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binmaxX-1, iby+1)); + } + } + + + + + new TCanvas(); + hietaY->Draw("colz"); + gPad->Modified(); + gPad->Update(); + + new TCanvas(); + hietaX->Draw("colz"); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaX->Draw(); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaY->Draw(); + gPad->Modified(); + gPad->Update(); + + + char foutName[1000]; + sprintf(foutName,"/local_zfs_raid/tomcat_20160528/trees/img_%s_eta_gmap_v2.root",name); + + + TFile *fout=new TFile(foutName, "RECREATE"); + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*25,0,400,200*25,100,300); + TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + TCanvas *c1=new TCanvas(); + imgLR->Draw("colz"); + TCanvas *c2=new TCanvas(); + imgHR->Draw("colz"); + TCanvas *c3=new TCanvas(); + imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + Double_t etamax=hietaX->GetXaxis()->GetBinCenter(hietaX->GetNbinsX()); + Double_t etamin=hietaX->GetXaxis()->GetBinCenter(0); + + char name1[100],name2[100],name3[100]; + sprintf(name1,"%sLR",name); + sprintf(name2,"%sHR",name); + sprintf(name3,"%sHR1",name); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (etaxetamax || etay>etamax || totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + imgHR->Fill(x-0.5+hietaX->GetBinContent(hietaX->GetXaxis()->FindBin(etax),hietaX->GetYaxis()->FindBin(etay))*sc-0.5*sc,y-0.5+hietaY->GetBinContent(hietaY->GetXaxis()->FindBin(etax),hietaY->GetYaxis()->FindBin(etay))*sc-0.5*sc); + imgHR1->Fill(x-0.5+h1etaX->GetBinContent(h1etaX->FindBin(etax))*sc-0.5*sc, y-0.5+h1etaY->GetBinContent(h1etaY->FindBin(etay)*sc-0.5*sc)); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + imgHR1->Write(name3,TObject::kOverwrite); + } + + } + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + imgHR1->Write(name3,TObject::kOverwrite); + + // c1->cd(); + // imgLR->DrawCopy("colz"); + // c2->cd(); + // imgHR->DrawCopy("colz"); + // c3->cd(); + // imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + + +TH2D *interpV3(char *name, int ip, double sc=26./25.) { + + char chainName[1000]; + sprintf(chainName,"/local_zfs_raid/tomcat_20160528/trees/%s_t*.root",name); + + const Double_t Delta_eta=0.05; + + + char gName[1000]; + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank_%d.root",ip); + TFile *fg=new TFile(gName); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + TH2F *heta=(TH2F*)fg->Get("heta"); + + TChain *ch=new TChain(name); + ch->Add(chainName); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + + + + TH2F *hietaX=(TH2F*)heta->Clone("hietaX"); + TH2F *hietaY=(TH2F*)heta->Clone("hietaY"); + hietaX->SetTitle("hietaX"); + hietaY->SetTitle("hietaY"); + hietaX->Reset(); + hietaY->Reset(); + + + + TH1D *hetaY=heta->ProjectionX(); + TH1D *hetaX=heta->ProjectionY(); + + TH1D *h1etaX=(TH1D*)hetaX->Clone("h1etaX"); + TH1D *h1etaY=(TH1D*)hetaY->Clone("h1etaY"); + h1etaX->SetTitle("h1etaX"); + h1etaY->SetTitle("h1etaY"); + h1etaX->Reset(); + h1etaY->Reset(); + + int binminX=heta->GetXaxis()->FindBin(-Delta_eta); + int binmaxX=heta->GetXaxis()->FindBin(1+Delta_eta); + int binminY=heta->GetYaxis()->FindBin(-Delta_eta); + int binmaxY=heta->GetYaxis()->FindBin(1+Delta_eta); + + + TH1D *hetaYp=heta->ProjectionX("py",binminY,binmaxY); + TH1D *hetaXp=heta->ProjectionY("px",binminX,binmaxX); + + for (int ibx=binminX; ibxSetBinContent(ibx, iby,hietaY->GetBinContent(ibx, iby-1)+heta->GetBinContent(ibx, iby)); + + hietaX->SetBinContent(ibx, iby,hietaX->GetBinContent(ibx-1, iby)+heta->GetBinContent(ibx, iby)); + + + + } + + } + + + for (int ibx=binminX; ibxSetBinContent(ibx,iby,hietaY->GetBinContent(ibx,iby)/hietaY->GetBinContent(ibx,binmaxY-1)); + hietaX->SetBinContent(ibx,iby,hietaX->GetBinContent(ibx,iby)/hietaX->GetBinContent(binmaxX-1,iby)); + + } + + } + + for (int iby=binminY; ibySetBinContent(iby+1,h1etaY->GetBinContent(iby)+hetaY->GetBinContent(iby+1)/hetaY->Integral(binminY,binmaxY)); + + for (int ibx=binminX; ibxSetBinContent(ibx+1,h1etaX->GetBinContent(ibx)+hetaX->GetBinContent(ibx+1)/hetaX->Integral(binminX,binmaxX)); + + for (int ibx=binmaxX; ibxGetNbinsX()+1; ibx++) + h1etaX->SetBinContent(ibx,1); + for (int iby=binmaxX; ibyGetNbinsY()+1; iby++) + h1etaY->SetBinContent(iby,1); + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaX->SetBinContent(ibx+1,iby+1,1); + } + } + + + for (int iby=0; ibyGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binminY+1)); + } + } + + + for (int iby=binmaxY-1; ibyGetNbinsY(); iby++) { + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaX->SetBinContent(ibx+1,iby+1,hietaX->GetBinContent(ibx+1,binmaxY-1)); + } + } + + + + for (int iby=binmaxX-1; ibyGetNbinsY(); iby++) { + + for (int ibx=0; ibxGetNbinsX(); ibx++) { + hietaY->SetBinContent(ibx+1,iby+1,1); + } + } + + + + for (int ibx=0; ibxGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binminX+1, iby+1)); + } + } + + + for (int ibx=binmaxX-1; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + hietaY->SetBinContent(ibx+1,iby+1,hietaY->GetBinContent(binmaxX-1, iby+1)); + } + } + + + + + new TCanvas(); + hietaY->Draw("colz"); + gPad->Modified(); + gPad->Update(); + + new TCanvas(); + hietaX->Draw("colz"); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaX->Draw(); + gPad->Modified(); + gPad->Update(); + new TCanvas(); + h1etaY->Draw(); + gPad->Modified(); + gPad->Update(); + + + char foutName[1000]; + sprintf(foutName,"/local_zfs_raid/tomcat_20160528/trees/img_%s_eta_gmap_v2.root",name); + + + TFile *fout=new TFile(foutName, "RECREATE"); + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + TH2D *imgLR=new TH2D("imgLR","imgLR",400,0,400,200,100,300); + TH2D *imgHR=new TH2D("imgHR","imgHR",400*25,0,400,200*25,100,300); + TH2D *imgHR1=new TH2D("imgHR1","imgHR1",400*25,0,400,200*25,100,300); + + + TCanvas *c1=new TCanvas(); + imgLR->Draw("colz"); + TCanvas *c2=new TCanvas(); + imgHR->Draw("colz"); + TCanvas *c3=new TCanvas(); + imgHR1->Draw("colz"); + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + Double_t etamax=hietaX->GetXaxis()->GetBinCenter(hietaX->GetNbinsX()); + Double_t etamin=hietaX->GetXaxis()->GetBinCenter(0); + + char name1[100],name2[100],name3[100]; + sprintf(name1,"%sLR",name); + sprintf(name2,"%sHR",name); + sprintf(name3,"%sHR1",name); + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + } + if (etaxetamax || etay>etamax || totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + imgLR->Fill(x,y); + imgHR->Fill(x-0.5+hietaX->GetBinContent(hietaX->GetXaxis()->FindBin(etax),hietaX->GetYaxis()->FindBin(etay))*sc,y-0.5+hietaY->GetBinContent(hietaY->GetXaxis()->FindBin(etax),hietaY->GetYaxis()->FindBin(etay))*sc); + imgHR1->Fill(x-0.5+h1etaX->GetBinContent(h1etaX->FindBin(etax))*sc, y-0.5+h1etaY->GetBinContent(h1etaY->FindBin(etay)*sc)); + } + if (ie%1000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + } + if (ie%10000000==0) { + // c1->Modified(); + // c1->Update(); + // c2->Modified(); + // c2->Update(); + // c3->Modified(); + // c3->Update(); + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + imgHR1->Write(name3,TObject::kOverwrite); + } + + } + imgLR->Write(name1,TObject::kOverwrite); + imgHR->Write(name2,TObject::kOverwrite); + imgHR1->Write(name3,TObject::kOverwrite); + + // c1->cd(); + // imgLR->DrawCopy("colz"); + // c2->cd(); + // imgHR->DrawCopy("colz"); + // c3->cd(); + // imgHR1->DrawCopy("colz"); + + fout->Close(); + +} + +TH2D *calcEtaGcorr(int ip) { + char name[1000]; + strcpy(name,"blank"); + char chainName[1000]; + sprintf(chainName,"/local_zfs_raid/tomcat_20160528/trees/%s_t*.root",name); + + const Double_t Delta_eta=0.05; + + + char gName[1000]; + sprintf(gName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_blank_%d.root",ip); + TFile *fg=new TFile(gName); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + + TH2F *heta=(TH2F*)fg->Get("heta"); + + TH2F *heta_new=(TH2F*)heta->Clone("heta_gmap"); + + heta_new->Reset(); + + TChain *ch=new TChain(name); + ch->Add(chainName); + + // TFile *ft=new TFile("/local_zfs_raid/tomcat_20160528/trees/grating_1d_t0_v4.root"); + // TTree *ch=(TTree*)ft->Get("grating_1d"); + + // ch->Draw("y:x>>hh(400,0,400,400,0,400)","","colz"); + // gPad->Modified(); + // gPad->Update(); + + + + char foutName[1000]; + sprintf(foutName,"/local_zfs_raid/tomcat_20160528/trees/gmap_eta_gcorr.root"); + + + TFile *fout=new TFile(foutName, "RECREATE"); + Double_t data[9], gain=1; + Double_t dum[9], sDum[2][2], etax, etay, sum, totquad; + Int_t x,y,f; + int ix, iy, skip; + + Long64_t nentries=ch->GetEntries(); + ch->SetBranchAddress("iFrame",&f); + ch->SetBranchAddress("x",&x); + ch->SetBranchAddress("y",&y); + ch->SetBranchAddress("data",dum); + Long64_t ie; + + + cout << "Chain has " << nentries*1E-6 << "M entries " << endl; + + + + for (ie=0; ieGetEntry(ie)) { + skip=0; + for (ix=0; ix<3; ix++) { + if (skip==0) { + for (iy=0; iy<3; iy++) { + if (gmap==NULL) { + data[ix+iy*3]=dum[ix+iy*3]; + } else { + gain=gmap->GetBinContent(gmap->GetXaxis()->FindBin(x), gmap->GetYaxis()->FindBin(y)); + if (gain>2000 && gain<3000) { + data[ix+iy*3]=dum[ix+iy*3]/gain; + // cout << dum[ix+iy*3] << " " << gain << " " << data[ix+iy*3] << endl; + } else + skip=1; + } + } + } + } + + } + if (skip==0) { + slsInterpolation::calcEta(data, etax, etay, sum, totquad, sDum); + if (totquad/sum<0.8 || totquad/sum>1.2) skip=1; + + if (skip==0) { + heta_new->Fill(etax,etay); + } + } + if (ie%10000000==0) { + cout << " " << ((float)ie)*100./((float)nentries)<< "%" << endl; + gmap->Write("gmap",TObject::kOverwrite); + heta_new->Write("heta",TObject::kOverwrite); + // imgHR1->Write(name3,TObject::kOverwrite); + } + + } + gmap->Write("gmap",TObject::kOverwrite); + heta_new->Write("heta",TObject::kOverwrite); + heta_new->DrawCopy("colz"); + fout->Close(); + +} diff --git a/slsDetectorCalibration/etaVEL/interpolation_EtaVEL.h b/slsDetectorCalibration/etaVEL/interpolation_EtaVEL.h new file mode 100644 index 000000000..93e2eaddd --- /dev/null +++ b/slsDetectorCalibration/etaVEL/interpolation_EtaVEL.h @@ -0,0 +1,52 @@ +#ifndef INTERPOLATION_ETAVEL_H +#define INTERPOLATION_ETAVEL_H + +#include +#include "EtaVEL.h" +#include "TH2F.h" +//#include "EtaVEL.cpp" +//class EtaVEL; + +class interpolation_EtaVEL: public slsInterpolation { + + public: + interpolation_EtaVEL(int nx=40, int ny=160, int ns=25, double etamin=-0.02, double etamax=1.02, int p=0); + ~interpolation_EtaVEL(); + + + //create eta distribution, eta rebinnining etc. + //returns flat field image + void prepareInterpolation(int &ok){prepareInterpolation(ok,10000);}; + void prepareInterpolation(int &ok, int maxit); + + //create interpolated image + //returns interpolated image + + //return position inside the pixel for the given photon + void getInterpolatedPosition(Int_t x, Int_t y, Double_t *data, Double_t &int_x, Double_t &int_y); + void getInterpolatedBin(Double_t *cluster, Int_t &int_x, Int_t &int_y); + + + + int addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay); + int addToFlatField(Double_t etax, Double_t etay); + int setPlot(int p=-1) {if (p>=0) plot=p; return plot;}; + int WriteH(){newEta->Write("newEta"); heta->Write("heta");}; + EtaVEL *setEta(EtaVEL *ev){if (ev) {delete newEta; newEta=ev;} return newEta;}; + TH2F *setEta(TH2F *ev){if (ev) {delete heta; heta=ev;} return heta;}; + void iterate(); + void DrawH(); + double getChiSq(){return newEta->getChiSq();}; + + + + protected: + EtaVEL *newEta; + TH2F *heta; + int plot; + + // ClassDefNV(interpolation_EtaVEL,1); + // #pragma link C++ class interpolation_EtaVEL-; +}; + +#endif diff --git a/slsDetectorCalibration/etaVEL/interpolation_etaVEL.cpp b/slsDetectorCalibration/etaVEL/interpolation_etaVEL.cpp new file mode 100644 index 000000000..1cdd3c718 --- /dev/null +++ b/slsDetectorCalibration/etaVEL/interpolation_etaVEL.cpp @@ -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++Modified(); + gPad->Update(); + } + if (newEta->converged==1) ok=1; else ok=0; +} + +int interpolation_EtaVEL::addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay) { + Double_t sum, totquad, sDum[2][2]; + int corner =calcEta(cluster, etax, etay, sum, totquad, sDum); + //check if it's OK...should redo it every time? + //or should we fill a finer histogram and afterwards re-fill the newEta? + addToFlatField(etax, etay); + return corner; +} + +int interpolation_EtaVEL::addToFlatField(Double_t etax, Double_t etay) { + // newEta->fill(etaX,etaY); + heta->Fill(etax,etay); + return 0; +} + +void interpolation_EtaVEL::iterate() { + cout << " -------------- newEta refilled"<< endl; + for (int ibx=0; ibxGetNbinsX(); ibx++) { + for (int iby=0; ibyGetNbinsY(); iby++) { + newEta->fill(heta->GetXaxis()->GetBinCenter(ibx+1),heta->GetYaxis()->GetBinCenter(iby+1),heta->GetBinContent(ibx+1,iby+1)); + } + } + newEta->updatePixelPos(); + cout << " -------------- pixelPosition updated"<< endl; +} + +void interpolation_EtaVEL::DrawH() { + heta->Draw("col"); + (newEta->plotPixelBorder())->Draw(); +} + + +void interpolation_EtaVEL::getInterpolatedPosition(Int_t x, Int_t y, Double_t *cluster, Double_t &int_x, Double_t &int_y) { + + Double_t etax, etay, sum, totquad, sDum[2][2]; + + int corner =calcEta(cluster, etax, etay, sum, totquad, sDum); + + int bin = newEta->findBin(etax,etay); + if (bin<=0) { + int_x=-1; + int_y=-1; + return; + } + double subX = ((double)(newEta->getXBin(bin))+.5)/((double)newEta->getNPixels()); + double subY = ((double)(newEta->getYBin(bin))+.5)/((double)newEta->getNPixels()); + + double dX, dY; + switch (corner) { + case TOP_LEFT: + dX=-1.; + dY=+1.; + break; + case TOP_RIGHT: + dX=+1.; + dY=+1.; + break; + case BOTTOM_LEFT: + dX=-1.; + dY=-1.; + break; + case BOTTOM_RIGHT: + dX=+1.; + dY=-1.; + break; + default: + dX=0; + dY=0; + } + + int_x=((double)x)+ subX+0.5*dX; + int_y=((double)y)+ subY+0.5*dY; + + // cout << corner << " " << subX<< " " << subY << " " << dX << " " << dY << " " << int_x << " " << int_y << endl; + +}; + + +// void interpolation_EtaVEL::Streamer(TBuffer &b){newEta->Streamer(b);}; +void interpolation_EtaVEL::getInterpolatedBin(Double_t *cluster, Int_t &int_x, Int_t &int_y) { + + Double_t etax, etay, sum, totquad, sDum[2][2]; + + int corner =calcEta(cluster, etax, etay, sum, totquad, sDum); + + int bin = newEta->findBin(etax,etay); + if (bin<0) { + int_x=-1; + int_y=-1; + return; + } + int_x=newEta->getXBin(bin); + int_y=newEta->getYBin(bin); + + + +}; + diff --git a/slsDetectorCalibration/etaVEL/slsInterpolation.h b/slsDetectorCalibration/etaVEL/slsInterpolation.h new file mode 100644 index 000000000..4a23ff8a1 --- /dev/null +++ b/slsDetectorCalibration/etaVEL/slsInterpolation.h @@ -0,0 +1,224 @@ +#ifndef SLS_INTERPOLATION_H +#define SLS_INTERPOLATION_H + +#include +#include +#include + +#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 diff --git a/slsDetectorCalibration/moench03TCtbData.h b/slsDetectorCalibration/moench03TCtbData.h new file mode 100644 index 000000000..6895078e1 --- /dev/null +++ b/slsDetectorCalibration/moench03TCtbData.h @@ -0,0 +1,154 @@ +#ifndef MOENCH03TCTBDATA_H +#define MOENCH03TCTBDATA_H +#include "slsDetectorData.h" + + + +class moench03CtbData : public slsDetectorData { + + 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(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=2*400*400) + cout << "Error: pointer " << dataMap[row][col] << " out of range "<< endl; + + } + } + for (int i=0; i0) { + iframe++; + // cout << ib << "-" << endl; + return (char*)afifo_cont; + } else { + delete [] afifo_cont; + return NULL; + } + } + return NULL; + }; + + + + +}; + + + +#endif diff --git a/slsDetectorCalibration/slsDetectorData.h b/slsDetectorCalibration/slsDetectorData.h index 3a6fae403..36f7a4218 100644 --- a/slsDetectorCalibration/slsDetectorData.h +++ b/slsDetectorCalibration/slsDetectorData.h @@ -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);}; /** diff --git a/slsDetectorCalibration/slsReceiverData.h b/slsDetectorCalibration/slsReceiverData.h index d40cbff4b..c66294cd3 100644 --- a/slsDetectorCalibration/slsReceiverData.h +++ b/slsDetectorCalibration/slsReceiverData.h @@ -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 (np0){ //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 <