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 "include/RTUtil.h"

double getNoise(TH2F*, int, int);

int biasScan(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)
{
  const int c_nChannels = 64;
  const double c_xOffset = 1; // mm
  const double c_yOffset = 0.7;
  
  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}
                };
                
                
  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");
  }
  
  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);
  canvas->cd();
  TVirtualPad *main = new TPad("main","main",0,0,1,1,10,1);
  main->Draw();
  main->Divide(4,4);
  for(int i=chXstart; i<chXend; i++) {
    for(int j=chYstart; j<chYend; j++) {
      int channel = map[i][j];
      TH1D* h_projection = h_threshold->ProjectionY("",channel+1,channel+1);
      int canvasPosition = i-chXend+4*(chYend-j)+1;
      printf(" %d ", canvasPosition);
      main->cd(canvasPosition);
      char name[128];
      sprintf(name,"Channel %d",channel);
      h_projection->SetTitle(name);
      h_projection->DrawCopy("ex0");
    }
  }
  printf("\n");
  }
  
  if(strstr(plopt, "1ch") != NULL)
  {
    TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");
    TCanvas* canvas = new TCanvas("canvas","",1600,1600);
    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->DrawCopy("ex0");
  }
  
  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();
  }
  
  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);
  }
  
  /** 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 = 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);
            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
    */
  if (strstr(plopt, "2d") != NULL) {
    
    int nX = chXend - chXstart + 1;
    int nY = chYend - chYstart + 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=chYstart; i<=chYend; i++) {
      for(int j=chXstart; j<=chXend; j++) {
      //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
      int canvasPosition = nX*(chYend-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 = 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[9];
    
    // 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, 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));
}