Subversion Repositories f9daq

#include "TROOT.h"

#include "TFile.h"

#include "TBenchmark.h"

#include "TH1F.h"

#include "TH2F.h"

#include "TH3F.h"

#include "TCanvas.h"

#include "TStyle.h"

#include "TPad.h"

#include "TF1.h"

#include "TGraph.h"

#include "TSpectrum.h"

#include "stdio.h"

#include "THStack.h"

#include "TPaveText.h"

#include "include/RTUtil.h"

double getNoise(TH2F*, int, int);

int biasScan(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, double par1=0.0, double par2=0.0, bool debug = false)

{

  const int c_nChannels = 64;

  const double c_xOffset = 0; // mm

  const double c_yOffset = 0;

  int map[8][8]={{32,34,53,55,40,42,61,63},

                 {48,50,37,39,56,58,45,47},

                 {33,35,52,54,41,43,60,62},

                 {49,51,36,38,57,59,44,46},

                 {17,19,4,6,25,27,12,14},

                 {1,3,20,22,9,11,28,30},

                 {16,18,5,7,24,26,13,15},

                 {0,2,21,23,8,10,29,31}

                };

  // Set draw style              

  RTSetStyle(gStyle);

  // open the file with histograms

  char fnameroot[256];

  TFile* rootfile;

        sprintf(fnameroot, "root/%s.root", filename);

        rootfile = (TFile *) gROOT->FindObject(filename);

        if(rootfile==NULL) rootfile = new TFile(fnameroot);

        if(rootfile==NULL) {

          printf("Cannot open root file %s!!!\n",fnameroot);

          return(0);

        }

  if(strstr(plopt, "tdc") != NULL) {

    TCanvas *canvas21 = new TCanvas("canvas21","canvas21",1600,800);

    TH3F* h0 = (TH3F*) rootfile->Get("h_correctedTDC");

          Int_t binsX = h0->GetXaxis()->GetNbins();

    Int_t minX  = h0->GetXaxis()->GetFirst();

    Int_t maxX  = h0->GetXaxis()->GetLast()+1;

    Int_t binsY = h0->GetYaxis()->GetNbins();

    Int_t minY  = h0->GetYaxis()->GetFirst();

    Int_t maxY  = h0->GetYaxis()->GetLast()+1;

    Int_t binsZ = h0->GetZaxis()->GetNbins();

    Int_t minZ  = h0->GetZaxis()->GetFirst();

    Int_t maxZ  = h0->GetZaxis()->GetLast()+1;

    Double_t xLowUser  = h0->GetXaxis()->GetBinLowEdge(minX);

    Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);

    Double_t yLowUser  = h0->GetYaxis()->GetBinLowEdge(minY);

    Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);

    canvas21->Divide(2);

    //h0->Draw();

    canvas21->cd(1);

    //h0->Project3D("xy")->Draw("colz");

    //canvas21->cd(2);

    //h0->Project3D("yz")->Draw("colz");

    //canvas21->cd(3);

    TH2D* h_correctedTDC = (TH2D*) h0->Project3D("xz");

    //h0->GetZaxis()->SetRangeUser(-5,5);

    //h0->Project3D("xzo")->Draw("colz");

    //h_correctedTDC->SetTitle("; t [ns]; Channel");

    h_correctedTDC->Draw("colz");

    /*

    TH1F* tdc1 = new TH1F("tdc1",";TDC [ns];Events",binsZ, minZ, maxZ);

    //for(int j=minY; j<maxY; j++) {

      for(int k=minZ; k<maxZ; k++) {

        double signal = h0->GetBinContent(1,1,k);

        tdc1->Fill(k,signal);

        //}

        }

    tdc1->Draw();

    */

    TH2D* h_timeWalk = (TH2D*) h0->Project3D("zy");

    canvas21->cd(2);

    //h_timeWalk->SetTitle(";Threshold [V]; t [ns]");

    h_timeWalk->Draw("colz");

  }

  if (strstr(plopt, "bias") != NULL) {

    TCanvas* canvas01 = new TCanvas("canvas01","",800, 800);

    canvas01->cd();

    TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");

    h_bias->Draw("colz");

    TCanvas* canvas02 = new TCanvas("canvas02","",800, 800);

    //canvas02->Divide(2,2);

    canvas02->cd();

    TH3F* h_biasThreshold = (TH3F*) rootfile->Get("h_biasThreshold");

    //h_biasThreshold->Draw();

    //canvas02->cd(2);

    //h_biasThreshold->GetXaxis()->SetRange(1,5);

    h_biasThreshold->Project3D("zy");

    h_biasThreshold_zy->DrawCopy("colz");

    TCanvas* canvas03 = new TCanvas("canvas03","",1200, 1200);

    canvas03->Divide(4,4);

    int x1=0;

    int x2=3;

    int y1=0;

    int y2=3;

    for(int j=x1; j<=x2; j++){

        for(int k=y1; k<=y2; k++){

                int channel = map[j][k];

                int canvasPosition = 4*(4-k-1) + 4 - j;

                canvas03->cd(canvasPosition);

                h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);

                char histoName[128];

                sprintf(histoName, "Ch %d", channel);

                h_biasThreshold->Project3D("zy");

                //sprintf(option, "h_biasThresholsd_zy%f", channel);

                h_biasThreshold_zy->SetTitle(histoName);

                h_biasThreshold_zy->DrawCopy("colz");

        }

    }

    TCanvas* canvas04 = new TCanvas("canvas04","",1200, 1200);

    canvas04->Divide(4,4);

    int x1=0;

        int x2=3;

        int y1=4;

        int y2=7;

        for(int j=x1; j<=x2; j++){

                for(int k=y1; k<=y2; k++){

                        int channel = map[j][k];

                        int canvasPosition = 4*(8-k-1) + 4 - j;

                        canvas04->cd(canvasPosition);

                        h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);

                        char histoName[128];

                        sprintf(histoName, "Ch %d", channel);

                        h_biasThreshold->Project3D("zy");

                        //sprintf(option, "h_biasThresholsd_zy%f", channel);

                        h_biasThreshold_zy->SetTitle(histoName);

                        h_biasThreshold_zy->DrawCopy("colz");

                }

        }

        TCanvas* canvas05 = new TCanvas("canvas05","",1200, 1200);

        canvas05->Divide(4,4);

        int x1=4;

        int x2=7;

        int y1=0;

        int y2=3;

        for(int j=x1; j<=x2; j++){

                for(int k=y1; k<=y2; k++){

                        int channel = map[j][k];

                        int canvasPosition = 4*(4-k-1) + 8 - j;

                        canvas05->cd(canvasPosition);

                        h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);

                        char histoName[128];

                        sprintf(histoName, "Ch %d", channel);

                        h_biasThreshold->Project3D("zy");

                        //sprintf(option, "h_biasThresholsd_zy%f", channel);

                        h_biasThreshold_zy->SetTitle(histoName);

                        h_biasThreshold_zy->DrawCopy("colz");

                }

        }

        TCanvas* canvas06 = new TCanvas("canvas06","",1200, 1200);

        canvas06->Divide(4,4);

        int x1=4;

        int x2=7;

        int y1=4;

        int y2=7;

        for(int j=x1; j<=x2; j++){

                for(int k=y1; k<=y2; k++){

                        int channel = map[j][k];

                        int canvasPosition = 4*(8-k-1) + 8 - j;

                        canvas06->cd(canvasPosition);

                        h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);

                        char histoName[128];

                        sprintf(histoName, "Ch %d", channel);

                        h_biasThreshold->Project3D("zy");

                        //sprintf(option, "h_biasThresholsd_zy%f", channel);

                        h_biasThreshold_zy->SetTitle(histoName);

                        h_biasThreshold_zy->DrawCopy("colz");

                }

        }

  }

  if (strstr(plopt, "th") != NULL) {

   TCanvas* canvas1 = new TCanvas("canvas1","",800,800);

   canvas1->cd();

   //gPad->SetLogz();

   //canvas1->SetLogz();

   TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");

   h_threshold->SetTitleOffset(1.5,"y");

   //h_threshold->GetZaxis()->SetRangeUser(0,1000);

   h_threshold->SetContour(20);

   h_threshold->Draw("colz");

   TCanvas* canvas2 = new TCanvas("canvas2","",1600,800);

   canvas2->Divide(2);

   TH1D* h_projection1 = h_threshold->ProjectionY("Ch 37",38,38);

   canvas2->cd(1);

   //h_projection1->GetYaxis()->SetRangeUser(0,1000);

   //gPad->SetLogy();

   h_projection1->Draw();

   TH1D* h_projection2 = h_threshold->ProjectionY("Ch 38",39,39);

   canvas2->cd(2);

   //gPad->SetLogy();

   h_projection2->Draw();

  }

  if(strstr(plopt,"16") != NULL)

  {

  TH2F* h_threshold = (TH2F*) rootfile->Get("h_bias");

  TCanvas* canvas = new TCanvas("canvas","",1600,1600);

  TCanvas* canvas11 = new TCanvas("canvas11","",800,800);

  canvas->cd();

  TVirtualPad *main = new TPad("main","main",0,0,1,1,10,1);

  main->Draw();

  main->cd();

  main->Divide(4,4);

  for(int i=chXstart; i<chXend; i++) {

    for(int j=(int)par1; j<(int)par2; j++) {

      int channel = map[i][j];

      TH1D* h_projection = h_threshold->ProjectionY("",channel+1,channel+1);

      int canvasPosition = i-chXend+4*((int)par2-j)+1;

      printf(" %d ", canvasPosition);

      main->cd(canvasPosition);

      gPad->SetMargin(0.08, 0.08, 0.08, 0.08);

      char name[128];

      sprintf(name,"Channel %d",channel);

      h_projection->SetTitle(name);

      h_projection->GetYaxis()->SetRangeUser(0,10000);

      h_projection->DrawCopy("pe1x0");

      canvas11->cd();

      if (canvasPosition == 1) {

        h_projection->SetMarkerStyle(kFullDotMedium);

        h_projection->GetYaxis()->SetRangeUser(0,10000);

        h_projection->DrawCopy("Pl");

        }

      else {

        h_projection->SetMarkerStyle(kFullDotMedium);

        h_projection->SetLineColor(canvasPosition);

        h_projection->SetMarkerColor(canvasPosition);

        h_projection->SetTitle("");

        h_projection->DrawCopy("plsame");

        }

    }

  }

  printf("\n");

  }

  if(strstr(plopt, "1ch") != NULL)

  {

    TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");

    TCanvas* canvas = new TCanvas("canvas","",800,800);

    canvas->cd();

    int channel = map[chXstart][chXend];

    TH1D* h_projection = h_bias->ProjectionY("",channel+1,channel+1);

    char name[128];

    sprintf(name,"Channel %d",channel);

    h_projection->SetTitle(name);

    h_projection->SetMarkerStyle(kFullDotMedium);

    //h_projection->SetMarkerSize(8);

    gPad->SetMargin(0.1,0.1,0.1,0.1);

    TF1* f_linear = new TF1("f_linear", "[0] + [1]*x", par1, par2);

    f_linear->SetParameter(0, -1);

    f_linear->SetParameter(1, 0.1);

    f_linear->SetParNames("p0", "p1");

    h_projection->Fit(f_linear,"q");

    h_projection->Fit(f_linear,"r");

    h_projection->Draw("pe1x0");

    double p0 = f_linear->GetParameter(0);

    double p1 = f_linear->GetParameter(1);

    printf("p0 %f p1 %f V_b = %f [V]\n", p0, p1, -p0/p1);

    TPaveText *text = new TPaveText(0.15,0.5,0.35,0.56);

    text->SetFillColor(0);

    text->SetBorderSize(0.1);

    text->SetTextSize(0.030);

    char string[128];

    sprintf(string, "V_b = %2.2f", -p0/p1);

    text->AddText(string);

    text->Draw("same");

    text->Paint();

  }

  if( strstr(plopt, "all") != NULL ) {

    TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);

    TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);

    canvas2->Divide(8,8);

    canvas3->Divide(8,8);

    TH1F* h_hitsx;

    TH1F* h_hitsy;

    for(int i=0; i<c_nChannels; i++) {

      canvas2->cd(i+1);

      char hname[128];

      sprintf(hname, "hnhitsx%d", i);

      h_hitsx = (TH1F*)rootfile->Get(hname);

      h_hitsx->Draw();

      canvas3->cd(i+1);

      sprintf(hname, "hnhitsy%d", i);

      h_hitsy = (TH1F*)rootfile->Get(hname);

      h_hitsy->Draw();  

    }

  }

  if( strstr(plopt, "x") != NULL ) {

    TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);

    TH1F* h_hitsx;

    canvas10->cd();

    char hname[128];

    sprintf(hname, "hnhitsx%d", chXstart);

    h_hitsx = (TH1F*)rootfile->Get(hname);

    h_hitsx->Draw("colz");

  }

  if( strstr(plopt, "y") != NULL ) {

    TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);

    TH1F* h_hitsy;

    canvas11->cd();

    char hname[128];

    sprintf(hname, "hnhitsy%d", chXstart);

    h_hitsy = (TH1F*)rootfile->Get(hname);

    h_hitsy->Draw();

  }

  /** Draws the sum of the channels

    * Each channel is a 2d plot

    * Intended for the study of 1 channel

    */

  if (strstr(plopt, "2d") != NULL) {

    int nX = chXend - chXstart + 1;

    int nY = (int)par2 - (int)par1 + 1;

    TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX*400,nY*400);

    printf("nx %d ny %d\n",nX,nY);

    canvas7->Divide(nX,nY);

    for(int i=(int)par1; i<=(int)par2; i++) {

      for(int j=chXstart; j<=chXend; j++) {

      //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;

      int canvasPosition = nX*((int)par2-i)+(j-chXstart)+1;

      if (debug) printf("canvas %d\n",canvasPosition);

        canvas7->cd(canvasPosition);

        char hname[128];

        int chPosition = map[j][i];

        sprintf(hname, "h2d%d", chPosition);

        TH2F* h_2d = (TH2F*)rootfile->Get(hname);

        h_2d->Draw("colz");

      } //x

    }

    // Number of photoelectrons - Poissonian correction

    TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);

    canvas8->cd();

    gStyle->SetOptStat(0);

      char hname[128];

      int chPosition = map[1][2];

      sprintf(hname, "h2d%d", chPosition);

      TH2F* h_2d = (TH2F*)rootfile->Get(hname);

      Int_t binsX = h_2d->GetXaxis()->GetNbins();

      Int_t minX  = h_2d->GetXaxis()->GetFirst();

      Int_t maxX  = h_2d->GetXaxis()->GetLast()+1;

      Int_t binsY = h_2d->GetYaxis()->GetNbins();

      Int_t minY  = h_2d->GetYaxis()->GetFirst();

      Int_t maxY  = h_2d->GetYaxis()->GetLast()+1;

      Double_t xLowUser  = h_2d->GetXaxis()->GetBinLowEdge(minX);

      Double_t xUpUser   = h_2d->GetXaxis()->GetBinUpEdge(maxX);

      Double_t yLowUser  = h_2d->GetYaxis()->GetBinLowEdge(minY);

      Double_t yUpUser =   h_2d->GetYaxis()->GetBinUpEdge(maxY);

      if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);

      TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);

      double noise = getNoise(h_2d, 1, 89);

      if(debug) printf("Noise = %f\n", noise);

      for(int k=minX; k<=maxX; k++) {

        for(int j=minY; j<=maxY; j++) {

          double signal = h_2d->GetBinContent(k,j); // detected

          //double p = ((signal - noise) > 1) ? (signal-noise) : 1;

          double p = signal - noise;

          p /= 10000.;

          double p0 = 1.0 - p; // events with zero photons

          //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);

          double eta = -log(p0); // constant of the poissonian statistics

          if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);

          double x = k*(xUpUser-xLowUser)/double(binsX);

          double y = j*(yUpUser-yLowUser)/double(binsY);

          h_corrected->Fill(x+3,y+8, eta);

        }

       }

      h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");

      h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);

      h_corrected->Draw("colz");

      // collection efficiency

      int nPoints =0;

      double efficiency=0;

      for (int i=18; i<=58; i++) {

        for (int j=19; j<=59; j++) {

          double signal = h_corrected->GetBinContent(i,j);

          if(debug) printf("signal %f\n",signal);

          efficiency += signal;

          nPoints++;

        }

      }  

    printf("Signal sum = %f\n   # of points = %d\n",efficiency,nPoints);

  }

  /** Draws the sum of channel signals

    * Each channel is a 2d ('h2d') histogram

    * Suitable for 8x8 chs scan

    */

  if( strstr(plopt, "sum") != NULL ) {

    int nX = 7 - 0 + 1;

    int nY = 7 - 0 + 1;

          TCanvas* canvas12 = new TCanvas("canvas12","c2",8*200, 8*200);

          canvas12->cd();

          gStyle->SetOptStat(0);

    // final histogram parameters

    TH2F* h0 = (TH2F*) rootfile->Get("h2d0");

          Int_t binsX = h0->GetXaxis()->GetNbins();

    Int_t minX  = h0->GetXaxis()->GetFirst();

    Int_t maxX  = h0->GetXaxis()->GetLast()+1;

    Int_t binsY = h0->GetYaxis()->GetNbins();

    Int_t minY  = h0->GetYaxis()->GetFirst();

    Int_t maxY  = h0->GetYaxis()->GetLast()+1;

    Double_t xLowUser  = h0->GetXaxis()->GetBinLowEdge(minX);

    Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);

    Double_t yLowUser  = h0->GetYaxis()->GetBinLowEdge(minY);

    Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);

    if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);

    TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);

    TH2F* h[9];

    // 2d histogram noise subtraction and poisson scaling

    for(int i=0; i<=7; i++) {

      for(int j=chXstart; j<=chXend; j++) {

        int chPosition = map[j][i];

        char hname[128];

              sprintf(hname, "h2d%d", chPosition);  

              int histogram = nX*(i-0)+(j-chXstart);

              h[histogram] = (TH2F *) rootfile->Get(hname);

              int noise = getNoise(h[histogram], 1, 170);

              if (debug) printf("noise: %d\n",noise);

              for(int k=minX; k<=maxX; k++) {

          for(int l=minY; l<=maxY; l++) {

            int signal = h[histogram]->GetBinContent(k,l); // detected

            //double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;

            double p = signal - noise;

            p /= 10000.;

            double p0 = 1.0 - p; // events with zero photons

            //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);

            double eta = -log(p0);

            //printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);

            double x = k*(xUpUser-xLowUser)/double(binsX);

            double y = l*(yUpUser-yLowUser)/double(binsY);

            h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);

          }

        }    

      }

    }

    h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");

    //h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);

    h_corrected->Draw("colz");

    TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);

    canvas13->Divide(2);

    canvas13->cd(1);

    //h[16]->Draw("colz");

    canvas13->cd(2);

    //h[8]->Draw("colz");

  }

  return(0);

}

/** Function calculates the noise from one channel

  * it runs through the bins along x and returns the average value

  */

double getNoise(TH2F* histogram, int yStart, int yEnd)

{

  double noise=0;

  int count=0;

  for(int j=yStart; j<yEnd; j++) {

    double value = histogram->GetBinContent(j,2);

    //if (noise < value) noise = value;

    noise += value;

    count++;

    }

  return (noise/double(count));

}