Subversion Repositories f9daq

#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();

}