Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  1. #include "TROOT.h"
  2. #include "TFile.h"
  3. #include "TBenchmark.h"
  4. #include "TH1F.h"
  5. #include "TH2F.h"
  6. #include "TCanvas.h"
  7. #include "TStyle.h"
  8. #include "TPad.h"
  9. #include "TF1.h"
  10. #include "TGraph.h"
  11. #include "TSpectrum.h"
  12. #include "TGraphErrors.h"
  13.  
  14. #define COROUT
  15. #define SHOW_SLICES
  16. #define NCH 20
  17. #define USED_CH 1
  18. #define TDC_BIN (25./1000.) //1 TDC bin in ns
  19. //~ #define FIRST_CH 1
  20.  
  21. double g_corpars[NCH][3];
  22.  
  23. int mapii[] = {  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19 };
  24. int mapij[] = {  1, 0,16,16,16,16,16,16,16,16,  6, 6, 6, 6, 6, 6, 6, 6, 6, 6 };
  25.  
  26. // root -l "cor.cpp(\"run045\",64,1,20,\"s\")"
  27. void cor_laser(char *fname="cscan3", const int nslicex=64, int slmin=1, const int projmin=5, char *fite="s")
  28. {
  29.         double fitw=1.0;
  30.         char skip='n';
  31.        
  32.         char fullname[256];
  33.        
  34.         //get ROOT file with histograms
  35.         char fnameroot[1024];
  36.         TFile * rootfile;
  37.         TDirectory *dir;
  38.        
  39.         sprintf(fnameroot, "root/%s.root", fname);
  40.         rootfile = (TFile *) gROOT->FindObject(fname);
  41.         if(rootfile==NULL) rootfile = new TFile(fnameroot);
  42.         if(rootfile==NULL) {
  43.           printf("Cannot open root file %s!!!\n",fnameroot);
  44.           return;
  45.         }
  46.         dir = (TDirectory*) rootfile;
  47.        
  48.        
  49.         // set draw style
  50.         gStyle->SetOptStat("ne");
  51.         gStyle->SetPalette(1, 0);
  52.        
  53.         gStyle->SetPaperSize(TStyle::kA4);
  54.         gStyle->SetStatBorderSize(1);
  55.         gStyle->SetFrameBorderMode(0);
  56.         gStyle->SetFrameFillColor(0);
  57.         gStyle->SetCanvasBorderMode(0);
  58.         gStyle->SetPadBorderMode(0);
  59.         gStyle->SetPadColor(0);
  60.         gStyle->SetCanvasColor(0);
  61.         gStyle->SetStatColor(0);
  62.         gStyle->SetOptFit(11);
  63.        
  64.         gStyle->SetPadRightMargin(0.1);
  65.        
  66.     char hname[256];
  67.     TH1F *hp1d; TH2F *hp2d;
  68.  
  69.     int ch=0;
  70.     int nbinsx, deltax, bin0, bin1;
  71.     double binc;
  72.     TH1F *hprojx[nslicex];  
  73.    
  74.     TF1 *f_fit = 0;
  75.     TF1 *fg = new TF1("fg", "gaus(0)");
  76.     fg->SetParNames("Constant","Mean","Sigma");
  77.     TF1 *fgg = new TF1("fgg", "gaus(0)+gaus(3)");
  78.     fgg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2");
  79.     TF1 *fggg = new TF1("fgg", "gaus(0)+gaus(3)+gaus(6)");
  80.     fggg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2","Constant3","Mean3","Sigma3");
  81.     TF1 *fitsqrt = new TF1("fitsqrt","[0]+[1]/TMath::Sqrt(x-[2])");
  82.    
  83.     TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)");
  84.     fcg->SetParNames("C", "Constant","Mean","Sigma");
  85.    
  86.     TGraphErrors *gcor[NCH];
  87.     double fitcenter, fitampl;
  88.  
  89.     double fitsigmas[2];
  90.     int fitminsigmai;
  91.     double fiterror, fitmean;
  92.  
  93. #ifdef COROUT
  94.   FILE *fp;
  95.   //open file to save parameters
  96.   sprintf(fullname, "data/%s_cor.txt", fname);
  97.   if( (fp=fopen(fullname, "wt")) == NULL )
  98.     printf("Cannot create parameter file %s !!!\n", fullname);
  99.   else
  100.     printf("Created parameter file %s\n", fullname);
  101. #endif  
  102.  
  103.  
  104.          TCanvas *c1 = new TCanvas("Cor.Fits", "Cor.Fits", 0, 0, 1300, 700);
  105.    c1->Divide(2,1);  
  106.  
  107.  
  108. #ifdef SHOW_SLICES
  109.         TCanvas *c0 = new TCanvas("Slices", "Slices", 0, 0, 1300, 700);
  110.         c0->Divide(5,5);
  111. #endif
  112.   const int c0pads = 25;    
  113.    
  114.              
  115.  if(nslicex==0) {
  116. #ifdef COROUT
  117.   for(int j=0;j<NCH; j++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);  
  118.   gSystem->Exit(1);
  119. #endif
  120.   } else {
  121. #ifdef FIRST_CH
  122.     ch=FIRST_CH;
  123. #else
  124.     for(ch=0;ch<USED_CH; ch++)    
  125. #endif    
  126.     {      
  127.      
  128.       gcor[ch] = new TGraphErrors(nslicex);
  129.       sprintf(hname,"TDC Corection ch%d;QDC;TDC [ps]",ch);
  130.       gcor[ch]->SetTitle(hname);
  131.       //~ sprintf(hname, "hdiffcor_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);       
  132.       sprintf(hname, "hcor%d",ch); hp2d = (TH2F *) dir->Get(hname);    
  133.       nbinsx=hp2d->GetNbinsX();
  134.       deltax=nbinsx/nslicex;
  135.      
  136.       c0->Clear();c0->Divide(8,4);
  137.         for(int i=slmin; i<nslicex-1;   i++) {
  138.           bin0=i*deltax+1;
  139.           bin1=(i+1)*deltax+1;
  140.           binc=(hp2d->GetXaxis()->GetBinCenter(bin0)+hp2d->GetXaxis()->GetBinCenter(bin1))/2.0;
  141.          
  142.           sprintf(hname, "hprojx%d",i);
  143.           hprojx[i]=(TH1F*)hp2d->ProjectionY(hname, bin0, bin1);
  144.          
  145.           //~ (hprojx[i]->GetXaxis())->SetRangeUser(-1, 2);
  146.         if( hprojx[i]->GetMaximum() > projmin ) {
  147.          
  148.           fitcenter = hprojx[i]->GetBinCenter(hprojx[i]->GetMaximumBin());
  149.           fitampl = hprojx[i]->GetMaximum();
  150.          
  151. //------/
  152.         if(i<1024) {
  153.            
  154.             //~ if( (ch<8) && (i<3) && (fitcenter < 0.4) ) fitcenter = 0.7;
  155.            
  156.             //~ f_fit = fcg;
  157.             //~ f_fit->SetParameters(fitampl*0.1, fitampl*0.9, fitcenter, 0.1);
  158.            
  159.             f_fit = fgg;
  160.             f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1, fitampl*0.1, fitcenter, 0.5);
  161.             //~
  162.             //~ f_fit = fg;
  163.             //~ f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1);
  164.  
  165.             //~ f_fit->SetParLimits(0, 0, fitampl*3);
  166.                  
  167.       #ifdef SHOW_SLICES                         
  168.             c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
  169.             hprojx[i]->Fit(f_fit,"QL", "", fitcenter-1.0*fitw, fitcenter+0.5*fitw);
  170.             //~ hprojx[i]->Fit(f_fit,"QL", "");
  171.                   //c0->Update();
  172.                   //if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
  173.             //if( (skip=='q') || (skip=='Q') ) return;
  174.             //hprojx[i]->DrawClone();
  175.             //fgg->SetRange(-1,2);
  176.             //fgg->DrawClone("LSAME");
  177.       #else
  178.             //~ hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
  179.             hprojx[i]->Fit(f_fit,"0QL", "");
  180.       #endif   
  181.      
  182.             //~ fitsigmas[0] = TMath::Abs(f_fit->GetParameter(2));
  183.             //~ fitsigmas[1] = TMath::Abs(f_fit->GetParameter(5));
  184.                   //~ fitminsigmai = TMath::LocMin(2,fitsigmas);
  185.                   //~ fitmean = f_fit->GetParameter(1 + fitminsigmai*3);
  186.                   //~ fiterror = TMath::Abs(f_fit->GetParameter(2 + fitminsigmai*3));
  187.                   //~ printf("fitminsigmai = %d {%lf, %lf}\n", fitminsigmai, fitsigmas[0], fitsigmas[1]);
  188.               //~
  189.                   //~ fitmean = f_fit->GetParameter(1);
  190.                   //~ fiterror = TMath::Abs(f_fit->GetParameter(2));
  191.               //~
  192.                   fitmean = f_fit->GetMaximumX(fitcenter-1.0*fitw, fitcenter+0.5*fitw);
  193.                   fiterror = TMath::Abs(f_fit->GetParameter(2));
  194. //------/
  195.         } else {          
  196.                   fgg->SetParameters(fitampl, fitcenter, 0.05,
  197.                                    fitampl/15, fitcenter, 0.5);
  198.            
  199.             fgg->SetParLimits(0, fitampl*0.3, fitampl*1.1);
  200.             fgg->SetParLimits(3, fitampl*0.01, fitampl*0.7);
  201.            
  202.             fgg->SetParLimits(1, -1, 1);
  203.             fgg->SetParLimits(4, -1, 1);
  204.            
  205.             fgg->SetParLimits(2, 0, 0.300);
  206.             fgg->SetParLimits(5, 0.3, 4.0);
  207.                  
  208.       #ifdef SHOW_SLICES                         
  209.             c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
  210.             hprojx[i]->Fit(fgg,"QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
  211.                   //c0->Update();
  212.                   //if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
  213.             //if( (skip=='q') || (skip=='Q') ) return;
  214.             //hprojx[i]->DrawClone();
  215.             //fgg->SetRange(-1,2);
  216.             //fgg->DrawClone("LSAME");
  217.       #else
  218.             hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
  219.       #endif   
  220.      
  221.             fitsigmas[0] = TMath::Abs(fgg->GetParameter(2));
  222.             fitsigmas[1] = TMath::Abs(fgg->GetParameter(5));
  223.                   fitminsigmai = TMath::LocMin(2,fitsigmas);
  224.                   //printf("fitminsigmai = %d {%.1lf, %.1lf}\n", fitminsigmai, fitsigmas[0]*1000, fitsigmas[1]*1000);
  225.                   fitmean = fgg->GetParameter(1 + fitminsigmai*3);
  226.                   fiterror = TMath::Abs(fgg->GetParameter(2 + fitminsigmai*3));
  227. //------/
  228.         } //if(i<4)
  229.                    
  230.                   printf(">>> ADC %3d - %3d -> %5.1lf ps\n", bin0, bin1, fiterror*1000.);
  231.                  
  232.                   gcor[ch]->SetPoint(i, binc, fitmean);  
  233.                   // weight = slice fit sigma (default)
  234.                   if( strcmp(fite, "s")==0 )              
  235.                     gcor[ch]->SetPointError(i, 0.0, fiterror);
  236.                   // weight = slice maximum bin content ~ n of events
  237.                   else if( strcmp(fite, "h")==0 )
  238.                     gcor[ch]->SetPointError(i, 0.0, 10./sqrt(hprojx[i]->GetMaximum()));
  239.                   // weight = both
  240.                   else if( strcmp(fite, "b")==0 )
  241.                     gcor[ch]->SetPointError(i, 0, 100.*fiterror/hprojx[i]->GetMaximum());
  242.                 } //if( hprojx[i]->GetMaximum() > projmin )            
  243.         }// for(int i=slmin; i<nslicex; i++)
  244. #ifdef SHOW_SLICES     
  245.       c0->Update();              
  246.         sprintf(fullname, "gif/%s_ch_%d.gif", fname, ch); c0->SaveAs(fullname);
  247. #endif 
  248.  
  249.  
  250.         c1->cd(ch+1); gPad->SetLogz();
  251.      
  252.       fitsqrt->SetLineWidth(2);
  253.       gStyle->SetOptFit();
  254.       (gcor[ch]->GetYaxis())->SetRangeUser((hp2d->GetYaxis())->GetXmin(), (hp2d->GetYaxis())->GetXmax());
  255.       //gcor[ch]->Draw("AP");  
  256.    
  257.       fitsqrt->SetParameters(0, 10, 50);
  258.      
  259.       //fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),(gcor[ch]->GetXaxis())->GetXmax());
  260.       fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),3*256.);
  261.       //fitsqrt->SetParameters(-0.25, 3, 55);
  262.       //fitsqrt->SetParLimits(0,-4,8);
  263.       //~ fitsqrt->SetParLimits(1,0,20);
  264.       //~ fitsqrt->SetParLimits(2,0,200);
  265. //~ #ifdef FIRST_CH      
  266.       //~ if(ch==FIRST_CH) gcor[ch]->Fit(fitsqrt, "0Q","");
  267. //~ #else      
  268.       //~ if(ch==0) gcor[ch]->Fit(fitsqrt, "0Q","");
  269. //~ #endif
  270.       gcor[ch]->Fit(fitsqrt, "0QR","");
  271.       //(gcor[ch]->GetYaxis())->UnZoom();      
  272.           hp2d->SetTitle("TDCvsQDC; QDC; TDC [ns]");
  273.           //~ (hp2d->GetXaxis())->SetRangeUser(0.,4*256.);
  274.           //(hp2d->GetYaxis())->SetRangeUser(-1.,3.);
  275.           //~ (hp2d->GetYaxis())->SetRangeUser(-3.,3.);
  276.                 hp2d->DrawCopy("COLZ");
  277.                 fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),4*256.);
  278.           fitsqrt->DrawCopy("LSAME");
  279.       gcor[ch]->Draw("PSAME"); 
  280.      
  281.       sprintf(fullname, "gif/%s_cor.gif", fname); c1->SaveAs(fullname);
  282.      
  283.           printf("[ch.%d - ch.%d] %lf %lf\n", mapii[ch], mapij[ch], fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
  284.         #ifdef COROUT
  285.           fprintf(fp,"%lf %lf %lf\n", fitsqrt->GetParameter(0), fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
  286.         #endif
  287.           g_corpars[ch][0]=fitsqrt->GetParameter(0);
  288.             g_corpars[ch][1]=fitsqrt->GetParameter(1);
  289.             g_corpars[ch][2]=fitsqrt->GetParameter(2);
  290.            
  291.      }         
  292.   }
  293.  
  294. #ifdef COROUT  
  295.         for(int i=USED_CH;i<NCH; i++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);
  296.         if(fp) fclose(fp);
  297. #endif
  298.  
  299.         //rootfile->Close();
  300. }
  301.