Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "TGraphErrors.h"
#include "TSystem.h"
#include "TMath.h"

#define COROUT
#define SHOW_SLICES
#define NCH 20
#define USED_CH 2
#define TDC_BIN (25./1000.) //1 TDC bin in ns
//~ #define FIRST_CH 0

double g_corpars[NCH][3];

int mapii[] = {  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19 };
int mapij[] = {  1, 0,16,16,16,16,16,16,16,16,  6, 6, 6, 6, 6, 6, 6, 6, 6, 6 };

// root -l "cor.cpp(\"run045\",64,1,20,\"s\")"
void cor(char *fname="cscan3", int range=1, const int nslicex=64, int slmin=1, const int projmin=5, char *fite="s")
{
        double fitw=1.0;
        char skip='n';
        char fullname[256];
       
        //get ROOT file with histograms
        char fnameroot[1024];
        TFile * rootfile;
        TDirectory *dir;
       
        sprintf(fnameroot, "root/%s.root", fname);
        rootfile = (TFile *) gROOT->FindObject(fname);
        if(rootfile==NULL) rootfile = new TFile(fnameroot);
        if(rootfile==NULL) {
          printf("Cannot open root file %s!!!\n",fnameroot);
          return;
        }
        dir = (TDirectory*) rootfile;
       
       
        // set draw style
        gStyle->SetOptStat("ne");
        gStyle->SetPalette(1, 0);
       
        gStyle->SetPaperSize(TStyle::kA4);
        gStyle->SetStatBorderSize(1);
        gStyle->SetFrameBorderMode(0);
        gStyle->SetFrameFillColor(0);
        gStyle->SetCanvasBorderMode(0);
        gStyle->SetPadBorderMode(0);
        gStyle->SetPadColor(0);
        gStyle->SetCanvasColor(0);
        gStyle->SetStatColor(0);
        gStyle->SetOptFit(11);
       
        gStyle->SetPadRightMargin(0.1);
       
    char hname[256];
    TH1F *hp1d; TH2F *hp2d;

    int ch=0;
    int nbinsx, deltax, bin0, bin1;
    double binc;
    TH1F *hprojx[nslicex];  
   
    TF1 *f_fit = 0;
    TF1 *fg = new TF1("fg", "gaus(0)");
    fg->SetParNames("Constant","Mean","Sigma");
    TF1 *fgg = new TF1("fgg", "gaus(0)+gaus(3)");
    fgg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2");
    TF1 *fggg = new TF1("fgg", "gaus(0)+gaus(3)+gaus(6)");
    fggg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2","Constant3","Mean3","Sigma3");
    TF1 *fitsqrt = new TF1("fitsqrt","[0]+[1]/TMath::Sqrt(x-[2])");
   
    TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)");
    fcg->SetParNames("C", "Constant","Mean","Sigma");
   
    TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)");
        fcgg->SetParName(0,"Const");
        fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
        fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
   
    TGraphErrors *gcor[NCH];
    double fitcenter, fitampl;

    double fitsigmas[2];
    int fitminsigmai;
    double fiterror, fitmean;

#ifdef COROUT
  FILE *fp;
  //open file to save parameters
  sprintf(fullname, "data/%s_cor.txt", fname);
  if( (fp=fopen(fullname, "wt")) == NULL )
    printf("Cannot create parameter file %s !!!\n", fullname);
  else
    printf("Created parameter file %s\n", fullname);
#endif  
 

         TCanvas *c1 = new TCanvas("Cor.Fits", "Cor.Fits", 0, 0, 1300, 450);
   c1->Divide(2,1);  

 
#ifdef SHOW_SLICES
        TCanvas *c0 = new TCanvas("Slices", "Slices", 0, 520, 1300, 470);
        c0->Divide(5,5);
#endif
  const int c0pads = 25;    
   
             
 if(nslicex==0) {
#ifdef COROUT
  for(int j=0;j<NCH; j++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);  
  gSystem->Exit(1);
  return;
#endif
  } else {
#ifdef FIRST_CH
    ch=FIRST_CH;
#else
    for(ch=0;ch<USED_CH; ch++)    
#endif    
    {      
     
      gcor[ch] = new TGraphErrors(nslicex);
      sprintf(hname,"TDC Corection ch%d;QDC;TDC [ps]",ch);
      gcor[ch]->SetTitle(hname);
      if(range == 0) {
                sprintf(hname, "hdiffcor_low_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);     
          } else if(range == 2) {
                sprintf(hname, "hdiffcor_hi_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);
          } else {
                sprintf(hname, "hdiffcor_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);
          }
      //~ sprintf(hname, "hcor%d",ch); hp2d = (TH2F *) dir->Get(hname);
      nbinsx=hp2d->GetNbinsX();
      deltax=nbinsx/nslicex;
     
      c0->Clear();c0->Divide(8,4);
        for(int i=slmin; i<nslicex;     i++) {
          bin0=i*deltax+1;
          bin1=(i+1)*deltax+1;
          binc=(hp2d->GetXaxis()->GetBinCenter(bin0)+hp2d->GetXaxis()->GetBinCenter(bin1))/2.0;
         
          sprintf(hname, "hprojx%d",i);
          hprojx[i]=(TH1F*)hp2d->ProjectionY(hname, bin0, bin1);
         
          //~ (hprojx[i]->GetXaxis())->SetRangeUser(-1, 2);
        if( hprojx[i]->GetMaximum() > projmin ) {
         
          fitcenter = hprojx[i]->GetBinCenter(hprojx[i]->GetMaximumBin());
          fitampl = hprojx[i]->GetMaximum();
         
//------/
        if(i<1024) {
           
            //~ if( (ch<8) && (i<3) && (fitcenter < 0.4) ) fitcenter = 0.7;
           
            f_fit = fcg;
            f_fit->SetParameters(fitampl*0.1, fitampl*0.9, fitcenter, 0.1);
            //~
            //~ f_fit = fcgg;
            //~ f_fit->SetParameters(fitampl*0.01,
                                                                        //~ fitampl*0.4, 0, 0.4,
                                                                        //~ fitampl*0.4, 0, 0.2);
           
            //~ f_fit = fgg;
            //~ f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1, fitampl*0.1, fitcenter, 0.5);
            //~
            //~ f_fit = fg;
            //~ f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1);

            //~ f_fit->SetParLimits(0, 0, fitampl*3);
                 
      #ifdef SHOW_SLICES                         
            c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
            //~ hprojx[i]->Fit(f_fit,"QL", "", fitcenter-1.0*fitw, fitcenter+0.5*fitw);
            hprojx[i]->Fit(f_fit,"QL", "");
                  //c0->Update();
                  //if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
            //if( (skip=='q') || (skip=='Q') ) return;
            //hprojx[i]->DrawClone();
            //fgg->SetRange(-1,2);
            //fgg->DrawClone("LSAME");
      #else
            //~ hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
            hprojx[i]->Fit(f_fit,"0QL", "");
      #endif   
     
            //~ fitsigmas[0] = TMath::Abs(f_fit->GetParameter(2));
            //~ fitsigmas[1] = TMath::Abs(f_fit->GetParameter(5));
                  //~ fitminsigmai = TMath::LocMin(2,fitsigmas);
                  //~ fitmean = f_fit->GetParameter(1 + fitminsigmai*3);
                  //~ fiterror = TMath::Abs(f_fit->GetParameter(2 + fitminsigmai*3));
                  //~ printf("fitminsigmai = %d {%lf, %lf}\n", fitminsigmai, fitsigmas[0], fitsigmas[1]);
              //~
                  //~ fitmean = f_fit->GetParameter(1);
                  //~ fiterror = TMath::Abs(f_fit->GetParameter(2));
              //~
             
                  // fcg
                  fitmean = f_fit->GetMaximumX(fitcenter-1.0*fitw, fitcenter+0.5*fitw);
                  fiterror = TMath::Abs(f_fit->GetParameter(3));
                 
                  //~ // fcgg
                  //~ fitmean = f_fit->GetMaximumX(fitcenter-1.0*fitw, fitcenter+1.0*fitw);
                  //~ double HalfMax = (f_fit->GetMaximum(fitcenter-1.0*fitw, fitcenter+1.0*fitw) + f_fit->GetParameter(0))/2.0;
                          //~ double MinHalfMax = f_fit->GetX(HalfMax, fitcenter-1.0*fitw, fitcenter);
                          //~ double MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitcenter+1.0*fitw);
                          //~ double FWHMc = MaxHalfMax - MinHalfMax;
                  //~ fiterror = TMath::Abs(FWHMc);
//~ //------/
        } else {          
                  fgg->SetParameters(fitampl, fitcenter, 0.05,
                                   fitampl/15, fitcenter, 0.5);
           
            fgg->SetParLimits(0, fitampl*0.3, fitampl*1.1);
            fgg->SetParLimits(3, fitampl*0.01, fitampl*0.7);
           
            fgg->SetParLimits(1, -1, 1);
            fgg->SetParLimits(4, -1, 1);
           
            fgg->SetParLimits(2, 0, 0.300);
            fgg->SetParLimits(5, 0.3, 4.0);
                 
      #ifdef SHOW_SLICES                         
            c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
            hprojx[i]->Fit(fgg,"QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
                  //c0->Update();
                  //if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
            //if( (skip=='q') || (skip=='Q') ) return;
            //hprojx[i]->DrawClone();
            //fgg->SetRange(-1,2);
            //fgg->DrawClone("LSAME");
      #else
            hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
      #endif   
     
            fitsigmas[0] = TMath::Abs(fgg->GetParameter(2));
            fitsigmas[1] = TMath::Abs(fgg->GetParameter(5));
                  fitminsigmai = TMath::LocMin(2,fitsigmas);
                  //printf("fitminsigmai = %d {%.1lf, %.1lf}\n", fitminsigmai, fitsigmas[0]*1000, fitsigmas[1]*1000);
                  fitmean = fgg->GetParameter(1 + fitminsigmai*3);
                  fiterror = TMath::Abs(fgg->GetParameter(2 + fitminsigmai*3));
//------/
        } //if(i<4)
                   
                  printf(">>> ADC %3d - %3d -> %5.1lf ps\n", bin0, bin1, fiterror*1000.);
                 
                  gcor[ch]->SetPoint(i, binc, fitmean);  
                  // weight = slice fit sigma (default)
                  if( strcmp(fite, "s")==0 )              
                    gcor[ch]->SetPointError(i, 0.0, fiterror);
                  // weight = slice maximum bin content ~ n of events
                  else if( strcmp(fite, "h")==0 )
                    gcor[ch]->SetPointError(i, 0.0, 2./sqrt(hprojx[i]->GetMaximum()));
                  // weight = both
                  else if( strcmp(fite, "b")==0 )
                    gcor[ch]->SetPointError(i, 0, 100.*fiterror/hprojx[i]->GetMaximum());
                  else if( strcmp(fite, "k")==0 )
                    gcor[ch]->SetPointError(i, 0, 1./sqrt(hprojx[i]->GetMaximum()));
               
                } //if( hprojx[i]->GetMaximum() > projmin )            
        }// for(int i=slmin; i<nslicex; i++)
#ifdef SHOW_SLICES     
      c0->Update();              
        sprintf(fullname, "gif/%s_ch_%d.gif", fname, ch); c0->SaveAs(fullname);
#endif 


        c1->cd(ch+1); gPad->SetLogz();
     
      fitsqrt->SetLineWidth(2);
        gStyle->SetOptStat("e");
      gStyle->SetOptFit(11);
      (gcor[ch]->GetYaxis())->SetRangeUser((hp2d->GetYaxis())->GetXmin(), (hp2d->GetYaxis())->GetXmax());
      //~ (gcor[ch]->GetYaxis())->SetRangeUser(-2,3);
      //gcor[ch]->Draw("AP");  
   
      fitsqrt->SetParameters(0, 10, 50);
     
      fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),(gcor[ch]->GetXaxis())->GetXmax());
      //~ fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),3*256.);
      //fitsqrt->SetParameters(-0.25, 3, 55);
      //fitsqrt->SetParLimits(0,-4,8);
      //~ fitsqrt->SetParLimits(1,0,20);
      //~ fitsqrt->SetParLimits(2,-100,2000);
      fitsqrt->SetParLimits(2,0,2000);
//~ #ifdef FIRST_CH      
      //~ if(ch==FIRST_CH) gcor[ch]->Fit(fitsqrt, "0Q","");
//~ #else      
      //~ if(ch==0) gcor[ch]->Fit(fitsqrt, "0Q","");
//~ #endif
      gcor[ch]->Fit(fitsqrt, "0QR","");
      //(gcor[ch]->GetYaxis())->UnZoom();  
      sprintf(fullname, "SiPM_%d; Charge [200 fC/bin]; TDC [ns]", ch+1);  
          hp2d->SetTitle(fullname);
          //~ (hp2d->GetXaxis())->SetRangeUser(0.,4*256.);
          //~ (hp2d->GetXaxis())->SetRangeUser(0.,3*256.);
          //(hp2d->GetYaxis())->SetRangeUser(-1.,3.);
          //~ (hp2d->GetYaxis())->SetRangeUser(-2.5,2.5);
                hp2d->DrawCopy("COLZ");
                fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),8*256.);
      //~ fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),(gcor[ch]->GetXaxis())->GetXmax());
          fitsqrt->DrawCopy("LSAME");
      gcor[ch]->Draw("PSAME"); 
     
      sprintf(fullname, "gif/%s_cor.gif", fname); c1->SaveAs(fullname);
     
          printf("[ch.%d - ch.%d] %lf %lf\n", mapii[ch], mapij[ch], fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
        #ifdef COROUT
          fprintf(fp,"%lf %lf %lf\n", fitsqrt->GetParameter(0), fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
        #endif
          g_corpars[ch][0]=fitsqrt->GetParameter(0);
            g_corpars[ch][1]=fitsqrt->GetParameter(1);
            g_corpars[ch][2]=fitsqrt->GetParameter(2);
           
     }         
  }
 
#ifdef COROUT  
        for(int i=USED_CH;i<NCH; i++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);
        if(fp) fclose(fp);
#endif

        //rootfile->Close();
}