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