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 "stdio.h"

#include "include/RTUtil.h"

double getNoise(TH2F*, int, int);

int sipm(char filename[256] = "test", char plopt[256]="all", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, 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}

                };

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

        }

        // set draw style

        RTSetStyle(gStyle);

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

  }

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

  }

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

   /*TCanvas *canvas4 = new TCanvas("canvas1","canvas1",1000,1000);

   int nChannels = chYend-chYstart+1;

   int ncols = nChannels/2;

   printf("nch %d nch\\2 %d\n", nChannels, ncols);

   canvas4->Divide(2,ncols);

   TH1F* h_hitsy;

   for(int i=chYstart; i<=chYend; i++){

     canvas4->cd(i-chYstart+1);

     char hname[128];

     int chPosition = map[0][i];

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

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

     h_hitsy->Draw();

   }*/

   TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);

   canvas4->cd();

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

     TH1F* h_hitsx;

     char hname[128];

     int chPosition = map[i][chYstart];

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

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

     h_hitsx->SetTitle("Scan X;x [mm]; Entries");

     h_hitsx->GetYaxis()->SetTitleOffset(1.3);

     h_hitsx->SetStats(0);

     if (i == chXstart)

      h_hitsx->Draw();

     else {

      h_hitsx->SetLineColor(i+1);

      h_hitsx->Draw("same");

     }

   }

   //sprintf(fullname, "ps/%s_Yshare.eps", filename);

         //canvas4->SaveAs(fullname);

   TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);

   canvas5->cd();

   for(int i=chYstart; i<=chYend; i++) {

     TH1F* h_hitsy;

     char hname[128];

     int chPosition = map[chXstart][i];

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

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

     h_hitsy->SetTitle("Scan Y;y [mm]; Entries");

     h_hitsy->GetYaxis()->SetTitleOffset(1.3);

     h_hitsy->SetStats(0);

     if (i == chYstart)

      h_hitsy->Draw();

     else {

      h_hitsy->SetLineColor(i+1);

      h_hitsy->Draw("same");

     }

   }

   //sprintf(fullname, "ps/%s_Yshare.eps", filename);

         //canvas5->SaveAs(fullname);

  }

  /** Draw one channel in x and y

    * Without and with the noise subtraction

    */

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

    TCanvas *canvas51 = new TCanvas("canvas51","canvas51",500,500);

    canvas51->cd();

    // Get the filled histogram from scan

    TH1F* h_x;

    char name[32];

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

    h_x = (TH1F*)rootfile->Get(name);

    h_x->DrawCopy();

    TH1F* h_y;

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

    h_y = (TH1F*)rootfile->Get(name);

    // Create and fill corrected histogram

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

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

    //Double_t xLowUser  = h_x->GetXaxis()->GetBinLowEdge(1);

    //Double_t xUpUser = h_x->GetXaxis()->GetBinUpEdge(binsX);

    Double_t xLowUser = h_x->GetXaxis()->GetXmin();

    Double_t xUpUser = h_x->GetXaxis()->GetXmax();

    TH1F* h_xCorrected = new TH1F("h_xCorrected",";;",binsX, xLowUser, xUpUser);

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

    double noise = h_x->GetBinContent(2);

    for (int jk = 0; jk < binsX; jk++) {

      double signal = h_x->GetBinContent(jk+1);

      signal -= noise;

      double x = h_x->GetXaxis()->GetBinCenter(jk+1);

      h_xCorrected->Fill(x, signal);

    }

    h_xCorrected->SetLineColor(2);

    h_xCorrected->DrawCopy("same");

    //Create and fill corrected histogram

    Int_t binsY = h_y->GetXaxis()->GetNbins();

    Double_t yLowUser = h_y->GetXaxis()->GetXmin();

    Double_t yUpUser = h_y->GetXaxis()->GetXmax();

    TH1F* h_yCorrected = new TH1F("h_yCorrected",";;",binsY, yLowUser, yUpUser);

    noise = h_y->GetBinContent(2);

    for (int jk = 0; jk < binsY; jk++) {

      double signal = h_y->GetBinContent(jk+1);

      signal -= noise;

      double x = h_y->GetXaxis()->GetBinCenter(jk+1);

      h_yCorrected->Fill(x, signal);

    }

    TCanvas *canvas52 = new TCanvas("canvas52","canvas52",500,500);

    canvas52->cd();

    h_y->DrawCopy();

    h_yCorrected->SetLineColor(2);

    h_yCorrected->DrawCopy("same");

  }

  /** Draws the signal from 8 channels in x-row

    * for one specific y bin, so the background and cross-talk

    * can be estimated.

    * Draws also a 2d scan of these channels.

    */

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

    TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);

    canvas6->cd(0);

    gStyle->SetOptStat(0);

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

    //! 1-dimension position in x vs. hits

    TH2F* h[8];

    TH1F* h_line[8];

    for(int j=0; j<8; j++) {

      h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);

    }

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

      int chPosition = map[j][chYstart];

      char hname[128];

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

      int histogram = j;

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

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

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

        int l=chYstart*20+12;

        //for(int l=12; l<=16; l++) {

          double signal = h[histogram]->GetBinContent(k,l);

          //signal -= noise;

          //signal /= 5*10000.0;

          double eta = -log(1 - signal);

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

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

          h_line[j]->Fill(x-c_xOffset, signal);

        //}

      }

      if (j == chXstart) {

        h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");

        //h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);

        //h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);

        h_line[j]->Draw("");

      }

      else {

        h_line[j]->SetLineColor(j+1);

        h_line[j]->Draw("same");

      }

    }

    //! 2d scan

    TCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);

    canvas61->cd();

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

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

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

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

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

        int chPosition = map[i][chYstart];

        char hname[128];

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

              int histogram = i;

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

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

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

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

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

            //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 signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;

            double p = signal - noise;

            p /= 10000.0;

            double eta = -log(1 - p);

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

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

            double x = h[histogram]->GetXaxis()->GetBinCenter(k);

            double y = h[histogram]->GetYaxis()->GetBinCenter(l);

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

          }

        }

    }

    h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");

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

    h_corrected->Draw("colz");

  }

  /** Draws the sum of the channels

    * Each channel is a 2d plot

    * Intended for the study of 1 channel

    *

    * Input: scanned channel coordinates (x,y)

    */

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

    int nX = 3;

    int nY = 3;

    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=chXend-1; i<=chXend+1; i++) {

      for(int j=chXstart-1; j<=chXstart+1; j++) {

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

      int canvasPosition = nX*(chXend+1-i) + (j-chXstart+1) +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[chXstart][chXend];

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

      Double_t xUpUser   = h_2d->GetXaxis()->GetXmax();

      Double_t yLowUser  = h_2d->GetYaxis()->GetXmin();

      Double_t yUpUser =   h_2d->GetYaxis()->GetXmax();

      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_corrected = new TH2F("h_corrected","h_corrected",binsX, minX, maxX, binsY, minY, maxY);

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

      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 = xLowUser + k*(xUpUser - xLowUser) / double(binsX);

          double x = h_2d->GetXaxis()->GetBinCenter(k);

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

          double y = h_2d->GetYaxis()->GetBinCenter(j);

          if (debug) printf("x=%f y=%f\n", x, y);

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

        }

       }

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

      h_corrected->SetTitle(";x[mm];y[mm]");

      h_corrected->GetZaxis()->SetRangeUser(-0.01,0.5);

      h_corrected->SetContour(50);

      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 = chXend - chXstart + 1;

    int nY = chYend - chYstart + 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[nX*nY];

    // 2d histogram noise subtraction and poisson scaling

    for(int i=chYstart; i<=chYend; 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-chYstart)+(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, eta);

          }

        }    

      }

    }

    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");

  }

  /** Draws the beam profile and fits it with error function

    * on some background function

    */

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

  TCanvas* canvas9 = new TCanvas("canvas9","canvas9", 500,500);

  canvas9->cd();

  char hname[128];

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

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

  h_laser->Draw();

  h_laser->SetStats(1);

  TF1* err = new TF1("err","[0]+[1]*TMath::Erf((x-[2])/[3])",17.02,17.30);

  err->SetParameter(0,2500);

  err->SetParameter(1, h_laser->GetMaximum());

  err->SetParameter(2, h_laser->GetBinCenter(h_laser->GetMaximumBin()));

  err->SetParameter(3, 0.001);

  h_laser->Fit(err,"qr");

  h_laser->Fit(err,"lr");

  double sigma = err->GetParameter(3);

  printf("sigma = %2.0f um,     FWHM = %2.0f um\n", sigma*1000, 2.35*sigma*1000);

  }

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

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

    canvas1->cd();

    TH2I* h_map = new TH2I("h_map","h_map",8,-0.5,7.5,8,-0.5,7.5);

    for (int i=7; i>=0; i--) {

      //for (int j=7; j>=0; j--) printf("(%d, %d) ", j,i);

      for (int j=7; j>=0; j--) {

        printf("%2d (%d %d)\n", map[j][i], j*10080+5040, i*10080+5040);

        h_map->Fill(j,i, map[j][i]);

      }

      printf("\n");

    }

    gStyle->SetOptStat(0);

    h_map->Draw("text");  

  }

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

}