diff --git a/slsDetectorCalibration/energyCalibration.cpp b/slsDetectorCalibration/energyCalibration.cpp index cca067d14..1e5df8c0c 100644 --- a/slsDetectorCalibration/energyCalibration.cpp +++ b/slsDetectorCalibration/energyCalibration.cpp @@ -447,9 +447,11 @@ TF1* energyCalibration::fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t fitfun= h1->GetFunction(fname); - fitfun->GetParameters(mypar); - for (int ip=0; ip<6; ip++) { - emypar[ip]=fitfun->GetParError(ip); + if (fitfun) { + fitfun->GetParameters(mypar); + for (int ip=0; ip<6; ip++) { + emypar[ip]=fitfun->GetParError(ip); + } } return fitfun; } diff --git a/slsDetectorCalibration/gainMap.C b/slsDetectorCalibration/gainMap.C index 3c7abef3c..2ef890454 100644 --- a/slsDetectorCalibration/gainMap.C +++ b/slsDetectorCalibration/gainMap.C @@ -16,22 +16,25 @@ using namespace std; #define FOPT "0" -TH2F *gainMap(TH2F *h2, float g) { +THStack *gainMap(TH2F *h2, float g) { // int npx, npy; int npx=160, npy=160; // det->getDetectorSize(npx, npy); + THStack *hs=new THStack(); + TH2F *gMap=new TH2F("gmap",h2->GetTitle(), npx, -0.5, npx-0.5, npy, -0.5, npy-0.5); + TH2F *nMap=new TH2F("nmap",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]; - + hs->Add(gMap); + hs->Add(nMap); + Double_t ens[1]={20.}, eens[1]={20.}; + Double_t peaks[1], epeaks[1]; int ibin; TH1D *px; - energyCalibration *enCal=new energyCalibration(); enCal->setPlotFlag(0); // enCal->setChargeSharing(0); @@ -50,94 +53,48 @@ TH2F *gainMap(TH2F *h2, float g) { 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(); + enCal->setFitRange(50,3000); + //enCal->setChargeSharing(0); + if (px) { + enCal->fixParameter(0,0); + enCal->fixParameter(1,0); + enCal->fixParameter(5,0); + mypar[0]=0; + mypar[1]=0; + mypar[2]=px->GetBinCenter(px->GetMaximumBin()); + mypar[3]=10; + mypar[4]=px->GetMaximum(); + mypar[5]=0; + enCal->setStartParameters(mypar); + enCal->fitSpectrum(px,mypar,emypar); - 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; - } + cout << ix << " " << iy << " " << mypar[2] << endl; } - - } + if (mypar[2]>0) { + gMap->SetBinContent(ix+1,iy+1,mypar[2]/ens[0]); + gMap->SetBinError(ix+1,iy+1,emypar[2]/ens[0]); + nMap->SetBinContent(ix+1,iy+1,mypar[3]); + nMap->SetBinError(ix+1,iy+1,emypar[3]); + } + } } - return gMap; + return hs; } -TH2F *noiseMap(TH2F *h2, TH2F *gMap=NULL) { +TH2F *noiseMap(TH2F *h2) { // int npx, npy; int npx=160, npy=160; // det->getDetectorSize(npx, npy); @@ -152,90 +109,37 @@ TH2F *noiseMap(TH2F *h2, TH2F *gMap=NULL) { cout << ix << " " << iy << " " << ibin << endl; ibin=h2->GetYaxis()->FindBin(ix*npy+iy); px=h2->ProjectionX("px",ibin,ibin); - px->Fit("gaus", FOPT); + px->Fit("gaus", FOPT,"",-100,100); if (px->GetFunction("gaus")) { - nMap->SetBinContent(ix+1,iy+1,px->GetFunction("gaus")->GetParameter(2)); + if (px->GetFunction("gaus")->GetParameter(1)>-5 && px->GetFunction("gaus")->GetParameter(1)<5) + nMap->SetBinContent(ix+1,iy+1,px->GetFunction("gaus")->GetParameter(2)); } // delete px; } } - if (gMap) { - nMap->Divide(gMap); - nMap->Scale(1000./3.6); - } + return nMap; } - -THStack *noiseHistos(char *tit) { - char fname[10000]; - - sprintf(fname,"/data/moench_xbox_201401_trees/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_201401_trees/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; - - -} - -THStack *noiseHistos(TH2F *nmap) { +THStack *noiseHistos(TH2F *nmap, TH2F *gmap=NULL) { char tit[1000]; + + if (nmap==NULL) { + cout << "No noise map" << endl; + + return NULL; + } + + if (gmap) { + nmap->Divide(gmap); + nmap->Scale(1000./3.6); + } + + strcpy(tit,nmap->GetTitle()); THStack *hs=new THStack(tit,tit); hs->SetTitle(tit); @@ -260,7 +164,7 @@ THStack *noiseHistos(TH2F *nmap) { 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 << " or " << h->GetFunction("gaus")->GetParameter(1) << "+-" << h->GetFunction("gaus")->GetParameter(2); } cout << endl; } diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index f829ed791..619126ca0 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 @@ -79,6 +79,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb TH2F *h3; TH2F *hetaX; TH2F *hetaY; + TH2D *clustHist; THStack *hs=new THStack("hs",fformat); @@ -117,13 +118,19 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb #ifdef MY_DEBUG - + quadrant quad; + tall->Branch("q",&quad,"q/I"); + + TCanvas *myC; TH2F *he; TCanvas *cH1; TCanvas *cH2; TCanvas *cH3; + int quadrants[5]; + for(int i = 0; i < 5; i++){ quadrants[i] = 0; } + if (hitfinder) { myC=new TCanvas(); he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax); @@ -138,6 +145,8 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb cH3=new TCanvas(); cH3->SetLogz(); h3->Draw("colz"); + + clustHist= new TH2D("clustHist","clustHist",3,-1.5,1.5,3,-1.5,1.5); } #endif filter->newDataSet(); @@ -145,7 +154,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb for (int irun=runmin; irunreadNextFrame(filebin))) { @@ -189,9 +198,17 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb h3->Fill(filter->getClusterTotal(3), iy+NR*ix); iFrame=decoder->getFrameNumber(buff); + +#ifdef MY_DEBUG + for(int cx=-1; cx <2; cx++) + for(int cy=-1; cy < 2; cy++) + clustHist->Fill(cx,cy,filter->getClusterElement(cx,cy)); + + quad = filter->getQuadrant(); + quadrants[quad]++; +#endif tall->Fill(); - - + } @@ -217,10 +234,10 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb #endif nf++; - cout << "=" ; + //cout << "=" ; delete [] buff; } - cout << nph << endl; + //cout << nph << endl; if (filebin.is_open()) filebin.close(); else @@ -233,6 +250,14 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb delete decoder; cout << "Read " << nf << " frames" << endl; + +#ifdef MY_DEBUG + cout << "quadrants: " ; + for(int i = 0; i<5; i++) cout << i << ": " << quadrants[i] << " || " ; + cout << endl; + cout << "Read Events " << nph << endl; +#endif + return hs; } diff --git a/slsDetectorCalibration/moenchReadDataMT.C b/slsDetectorCalibration/moenchReadDataMT.C index 2411e73ae..54773f6b8 100644 --- a/slsDetectorCalibration/moenchReadDataMT.C +++ b/slsDetectorCalibration/moenchReadDataMT.C @@ -9,23 +9,30 @@ typedef struct task_s{ int runmin; int runmax; int treeIndex; + int tNumber; + double xTalkFactor; } Task; +double hc = 0.1; // read - out crosstalk + void *moenchMakeTreeTask(void *p){ TThread::Lock(); char fname[1000]; Task *t = (Task *)p; sprintf(fname,"%s%s_%i.root",t->tdir,t->tname,t->treeIndex); TFile *f = new TFile(fname,"RECREATE"); - cout << "Call moenchReadData(" << t->fformat << "," << t->tname << "," << t->runmin<< "," << t->runmax <<")" << endl; + double xTalkC = t->xTalkFactor; //((double)t->tNumber) * 0.01; Calibrated with data in /data/moench_trieste_20140313_trees/xtalkScan + cout << "Call moenchReadData(" << t->fformat << "," << t->tname << "," << t->runmin<< "," << t->runmax <<") with xTalkC: " << xTalkC << endl; TThread::UnLock(); - moenchReadData(t->fformat,t->tname,t->runmin,t->runmax); - f->Close(); + moenchReadData(t->fformat,t->tname,t->runmin,t->runmax, 1500, -500, 1000, 1, xTalkC); + if(f && f->IsOpen()){ + f->Close(); + } return 0; } -void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runoffset, int nThreads, int treeIndexStart=0){ +void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runoffset, int nThreads, int treeIndexStart=0, double xTalkFactor=0.044){ char threadName[1000]; TThread *threads[nThreads]; for(int i = 0; i < nThreads; i++){ @@ -37,6 +44,8 @@ void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runo t->runmin = runmin + i*runoffset; t->runmax = runmin + (i+1)*runoffset - 1; t->treeIndex = treeIndexStart + i; + t->tNumber = i; + t->xTalkFactor = xTalkFactor; cout << "start thread " << i << " start: " << t->runmin << " end " << t->runmax << endl; threads[i] = new TThread(threadName, moenchMakeTreeTask, t); threads[i]->Run(); @@ -45,15 +54,17 @@ void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runo for(int i = 0; i < nThreads; i++){ threads[i]->Join(); } - + + cout << " ( DONE ) " << endl; + } //to compile: g++ -DMYROOT -DMYROOT1 -g `root-config --cflags --glibs` -o moenchReadDataMT moenchReadDataMT.C int main(int argc, char **argv){ - if(argc != 8){ - cout << "Usage: " << argv[0] << " fformat tit tdir runmin runoffset nThreads treeIndexStart" << endl; + if(argc < 8){ + cout << "Usage: " << argv[0] << " fformat tit tdir runmin runoffset nThreads treeIndexStart [xTalkFactor]" << endl; exit(-1); } @@ -64,8 +75,12 @@ int main(int argc, char **argv){ int runoffset = atoi(argv[5]); int nThreads = atoi(argv[6]); int treeIndexStart = atoi(argv[7]); + double xTalkFactor = 0.044; - moenchReadDataMT(fformat, tit, tdir,runmin,runoffset,nThreads,treeIndexStart); + if(argc == 9) + xTalkFactor = atof(argv[8]); + + moenchReadDataMT(fformat, tit, tdir,runmin,runoffset,nThreads,treeIndexStart, xTalkFactor); } diff --git a/slsDetectorCalibration/raedNoiseData.C b/slsDetectorCalibration/raedNoiseData.C index c4b505df3..141c0f6ef 100644 --- a/slsDetectorCalibration/raedNoiseData.C +++ b/slsDetectorCalibration/raedNoiseData.C @@ -14,12 +14,12 @@ void raedNoiseData(char *tit, int sign=1){ TFile *fout; THStack *hs2N; - sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s.root",tit); + sprintf(fname,"/data/moench_trieste_calibration_trees/flat_20keV_%s.root",tit); fout=new TFile(fname,"RECREATE"); - sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s_f00000%%04d000_0.raw",tit); + sprintf(fname,"/data/moench_trieste_calibration/flat_20keV_%s_f00000%%04d000_0.raw",tit); - hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0,1); + hs2N=moenchReadData(fname,tit,0,2000,1500,-500,2500,sign,0.,1,159,1,159, 0,1); hs2N->SetName(tit); hs2N->SetTitle(tit); (TH2F*)(hs2N->GetHists()->At(0))->Write(); @@ -100,10 +100,10 @@ void raedNoiseDataN(char *tit, int sign=1){ TFile *fout; THStack *hs2N; - sprintf(fname,"/data/moench_xbox_20140116/noise_%s.root",tit); + sprintf(fname,"/data/moench_noise_20140327_trees/noise_%s.root",tit); fout=new TFile(fname,"RECREATE"); - sprintf(fname,"/data/moench_xbox_20140116/noise_%s_f00000%%04d000_0.raw",tit); + sprintf(fname,"/data/moench_noise_20140327/noise_%s_f00000%%04d000_0.raw",tit); hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0,0); hs2N->SetName(tit); diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h index 24d58f929..0f382ec9c 100644 --- a/slsDetectorCalibration/singlePhotonDetector.h +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -9,7 +9,7 @@ #include "commonModeSubtraction.h" -//#define MYROOT1 +#define MYROOT1 #ifdef MYROOT1 #include