diff --git a/slsDetectorCalibration/demoCreateTree.C b/slsDetectorCalibration/demoCreateTree.C new file mode 100644 index 000000000..d14aa8b02 --- /dev/null +++ b/slsDetectorCalibration/demoCreateTree.C @@ -0,0 +1,31 @@ +{ + + //.L moenchReadData.C + + + + + + + TFile *fout; + THStack *hs2N; + + fout=new TFile("/scratch/outfile.root","RECREATE"); + + hs2N=moenchReadData("/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_cds_g4_f00000%%04d000_0.raw",0,20,1500,-500,2500,1,0.,1,159,1,159, 0,1); + hs2N->SetName("cds_g4"); + hs2N->SetTitle("cds_g4"); + (TH2F*)(hs2N->GetHists()->At(0))->Write(); + + (TH2F*)(hs2N->GetHists()->At(1))->Write(); + (TH2F*)(hs2N->GetHists()->At(2))->Write(); + (TH2F*)(hs2N->GetHists()->At(3))->Write(); + (TH2F*)(hs2N->GetHists()->At(4))->Write(); + + + fout->Close(); + + + +} + diff --git a/slsDetectorCalibration/energyCalibration.cpp b/slsDetectorCalibration/energyCalibration.cpp index df548f1cf..cca067d14 100644 --- a/slsDetectorCalibration/energyCalibration.cpp +++ b/slsDetectorCalibration/energyCalibration.cpp @@ -289,7 +289,9 @@ energyCalibration::~energyCalibration(){ -TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch) { + + +TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch, int direction) { if (h2==NULL || nch==0) return NULL; @@ -299,19 +301,28 @@ TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch) { double val=-1; - h1=new TH1F("median","Median",h2->GetYaxis()->GetNbins(),h2->GetYaxis()->GetXmin(),h2->GetYaxis()->GetXmax()); - - - for (int ib=0; ibGetXaxis()->GetNbins(); ib++) { - for (int ich=0; ichGetBinContent(ch0+ich+1,ib+1); + if (direction==0) { + h1=new TH1F("median","Median",h2->GetYaxis()->GetNbins(),h2->GetYaxis()->GetXmin(),h2->GetYaxis()->GetXmax()); + for (int ib=0; ibGetXaxis()->GetNbins(); ib++) { + for (int ich=0; ichGetBinContent(ch0+ich+1,ib+1); + } + val=energyCalibrationFunctions::median(x, nch); + h1->SetBinContent(ib+1,val); + } + } else if (direction==1) { + h1=new TH1F("median","Median",h2->GetXaxis()->GetNbins(),h2->GetXaxis()->GetXmin(),h2->GetXaxis()->GetXmax()); + for (int ib=0; ibGetYaxis()->GetNbins(); ib++) { + for (int ich=0; ichGetBinContent(ib+1,ch0+ich+1); + } + val=energyCalibrationFunctions::median(x, nch); + h1->SetBinContent(ib+1,val); } - val=energyCalibrationFunctions::median(x, nch); - h1->SetBinContent(ib+1,val); } - return h1; + delete [] x; - + return h1; } @@ -480,6 +491,7 @@ TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Doub eoff=fitfun->GetParError(0); gain=funcs->setScanSign()*mypar[1]; + off=mypar[0]; return gr; diff --git a/slsDetectorCalibration/energyCalibration.h b/slsDetectorCalibration/energyCalibration.h index 904705aa3..7afda027a 100644 --- a/slsDetectorCalibration/energyCalibration.h +++ b/slsDetectorCalibration/energyCalibration.h @@ -251,7 +251,15 @@ class energyCalibration { #ifdef MYROOT - static TH1F* createMedianHistogram(TH2F* h2, int ch0, int nch); + /** + Creates an histogram with the median of nchannels starting from a specified one. the direction on which it is mediated can be selected (defaults to x=0) + \param h2 2D histogram on which the median will be calculated + \param ch0 starting channel + \param nch number of channels to be mediated + \param direction can be either 0 (x, default) or 1 (y) + \returns a TH1F histogram with the X-axis as a clone of the h2 Y (if direction=0) or X (if direction=0) axis, and on the Y axis the median of the counts of the mediated channels f h2 + */ + static TH1F* createMedianHistogram(TH2F* h2, int ch0, int nch, int direction=0); /** sets the s-curve fit range diff --git a/slsDetectorCalibration/gMapDemo.C b/slsDetectorCalibration/gMapDemo.C new file mode 100644 index 000000000..51addf69a --- /dev/null +++ b/slsDetectorCalibration/gMapDemo.C @@ -0,0 +1,11 @@ +{ + //.L energyCalibration.cpp+ + //.L gainMap.C+ + TFile fin("/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_cds_g4.root"); + TH2F *h2=fin.Get("h2"); + TH2F *gMap=gainMap(h2,4); + gMap->Draw("colz"); + + + +} diff --git a/slsDetectorCalibration/gainMap.C b/slsDetectorCalibration/gainMap.C new file mode 100644 index 000000000..7923fe687 --- /dev/null +++ b/slsDetectorCalibration/gainMap.C @@ -0,0 +1,228 @@ +#define MYROOT +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + + +#define FOPT "0" + +TH2F *gainMap(TH2F *h2, float g) { + // int npx, npy; + int npx=160, npy=160; + // det->getDetectorSize(npx, npy); + + TH2F *gMap=new TH2F("gmap",h2->GetTitle(), npx, -0.5, npx-0.5, npy, -0.5, npy-0.5); + + Double_t ens[3]={0,8,17.5}, eens[3]={0.,0.1,0.1}; + Double_t peaks[3], epeaks[3]; + + + + int ibin; + TH1D *px; + + + energyCalibration *enCal=new energyCalibration(); + enCal->setPlotFlag(0); + // enCal->setChargeSharing(0); + enCal->setScanSign(1); + + Double_t gain, off, egain, eoff; + + + TList *functions; + TPolyMarker *pm ; + Double_t *peakX, *peakY; + TSpectrum *s=new TSpectrum(); + Double_t mypar[10], emypar[10]; + Double_t prms, np; + int iit=0; + TGraph *glin; + Double_t peakdum, hpeakdum; + + for (int ix=1; ixProjectionX("px",ibin+1,ibin+1); + prms=10*g; + iit=0; + np=s->Search(px,prms,"",0.2); + while (np !=2) { + if (np>2) + prms+=0.5*prms; + else + prms-=0.5*prms; + iit++; + if (iit>=10) + break; + np=s->Search(px,prms,"",0.2); + } + if (np!=2) + cout << "peak search could not converge " << ibin << endl; + if (np==2) { + pm=NULL; + functions=px->GetListOfFunctions(); + if (functions) + pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); + if (pm) { + peakX=pm->GetX(); + peakY=pm->GetY(); + + + if (peakX[0]>peakX[1]) { + peakdum=peakX[0]; + hpeakdum=peakY[0]; + peakX[0]= peakX[1]; + peakY[0]= peakY[1]; + peakX[1]= peakdum; + peakY[1]= hpeakdum; + + } + + cout << "("<< ix << "," << iy << ") " << endl; + for (int ip=0; ipsetFitRange(peakX[ip]-10*g,peakX[ip]+10*g); + mypar[0]=0; + mypar[1]=0; + mypar[2]=peakX[ip]; + mypar[3]=g*10; + mypar[4]=peakY[ip]; + mypar[5]=0; + + + + enCal->setStartParameters(mypar); + enCal->fitSpectrum(px,mypar,emypar); + + + peaks[ip+1]=mypar[2]; + epeaks[ip+1]=emypar[2]; + } + + peaks[0]=0; + epeaks[0]=1; + + // for (int i=0; i<3; i++) cout << i << " " << ens[i] << " " << eens[i]<< " "<< peaks[i]<< " " << epeaks[i] << endl; + + glin= enCal->linearCalibration(3,ens,eens,peaks,epeaks,gain,off,egain,eoff); + + // cout << "Gain " << gain << " off " << off << endl; + if (off>-10 && off<10) { + gMap->SetBinContent(ix+1,iy+1,gain); + gMap->SetBinError(ix+1,iy+1,egain); + } + if (glin) + delete glin; + } + } + + + } + } + + return gMap; +} + + +TH2F *noiseMap(TH2F *h2) { + // int npx, npy; + int npx=160, npy=160; + // det->getDetectorSize(npx, npy); + + TH2F *nMap=new TH2F("nmap",h2->GetTitle(), npx, -0.5, npx-0.5, npy, -0.5, npy-0.5); + + int ibin; + TH1D *px; + + for (int ix=0; ixProjectionX("px",ibin+1,ibin+1); + px->Fit("gaus", FOPT); + if (px->GetFunction("gaus")) { + nMap->SetBinContent(ix+1,iy+1,px->GetFunction("gaus")->GetParameter(2)); + } + // delete px; + } + } + + return nMap; +} + + +THStack *noiseHistos(char *tit) { + char fname[10000]; + + sprintf(fname,"/data/moench_xbox_20140116/noise_map_%s.root",tit); + TFile *fn=new TFile(fname); + TH2F *nmap=(TH2F*)fn->Get("nmap"); + + if (nmap==NULL) { + cout << "No noise map in file " << fname << endl; + + return NULL; + } + + sprintf(fname,"/data/moench_xbox_20140113/gain_map_%s.root",tit); + TFile *fg=new TFile(fname); + TH2F *gmap=(TH2F*)fg->Get("gmap"); + + if (gmap==NULL) { + cout << "No gain map in file " << fname << endl; + + return NULL; + } + + nmap->Divide(gmap); + nmap->Scale(1000./3.6); + + THStack *hs=new THStack(tit,tit); + hs->SetTitle(tit); + + TH1F *h; + char hname[100]; + + cout << tit << endl; + for (int is=0; is<4; is++) { + sprintf(hname,"h%ds",is+1); + + h=new TH1F(hname,tit,500,0,500); + hs->Add(h); + // cout << hs->GetHists()->GetEntries() << endl; + for (int ix=40*is+2; ix<40*(is+1)-2; ix++) { + + for (int iy=2; iy<158; iy++) { + if (ix<100 || ix>120) + h->Fill(nmap->GetBinContent(ix+1,iy+1)); + } + } + cout << is+1 << "SC: " << "" << h->GetMean() << "+-" << h->GetRMS(); + h->Fit("gaus","0Q"); + h->SetLineColor(is+1); + if (h->GetFunction("gaus")) { + h->GetFunction("gaus")->SetLineColor(is+1); + cout << " or " << h->GetFunction("gaus")->GetParameter(1) << "+-" << h->GetFunction("gaus")->GetParError(1); + } + cout << endl; + } + + // cout << hs->GetHists()->GetEntries() << endl; + + return hs; + + +} diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index e98fea494..bf7173f0b 100644 --- a/slsDetectorCalibration/moenchReadData.C +++ b/slsDetectorCalibration/moenchReadData.C @@ -26,7 +26,7 @@ using namespace std; #define NR 160 -//#define MY_DEBUG 1 +#define MY_DEBUG 1 #ifdef MY_DEBUG #include #endif diff --git a/slsDetectorCalibration/raedNoiseData.C b/slsDetectorCalibration/raedNoiseData.C index 0e7f24f21..b2021b29f 100644 --- a/slsDetectorCalibration/raedNoiseData.C +++ b/slsDetectorCalibration/raedNoiseData.C @@ -1,21 +1,24 @@ #include "moenchReadData.C" + void raedNoiseData(char *tit, int sign=1){ + + char fname[1000]; char f[1000]; TFile *fout; THStack *hs2N; - sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s.root",tit); + sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s_0.root",tit); fout=new TFile(fname,"RECREATE"); sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s_f00000%%04d000_0.raw",tit); - hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0); + hs2N=moenchReadData(fname,0,3000,1500,-500,2500,1,0.,1,159,1,159, 0,1); hs2N->SetName(tit); hs2N->SetTitle(tit); (TH2F*)(hs2N->GetHists()->At(0))->Write();