Subversion Repositories f9daq

Rev

Rev 30 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #include "TROOT.h"
  2. #include "TFile.h"
  3. #include "TBenchmark.h"
  4. #include "TH1F.h"
  5. #include "TH2F.h"
  6. #include "TCanvas.h"
  7. #include "TStyle.h"
  8. #include "TPad.h"
  9. #include "TF1.h"
  10. #include "TGraph.h"
  11. #include "TSpectrum.h"
  12. #include "stdio.h"
  13.  
  14. #include "include/RTUtil.h"
  15.  
  16. double getNoise(TH2F*, int, int);
  17.  
  18. int plotThreshold(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)
  19. {
  20.   const int c_nChannels = 64;
  21.   const double c_xOffset = 1; // mm
  22.   const double c_yOffset = 0.7;
  23.  
  24.   int map[8][8]={{32,34,53,55,40,42,61,63},
  25.                  {48,50,37,39,56,58,45,47},
  26.                  {33,35,52,54,41,43,60,62},
  27.                  {49,51,36,38,57,59,44,46},
  28.                  {17,19,4,6,25,27,12,14},
  29.                  {1,3,20,22,9,11,28,30},
  30.                  {16,18,5,7,24,26,13,15},
  31.                  {0,2,21,23,8,10,29,31}
  32.                 };
  33.  
  34.   char fnameroot[256];
  35.   TFile* rootfile;
  36.         sprintf(fnameroot, "root/%s.root", filename);
  37.         rootfile = (TFile *) gROOT->FindObject(filename);
  38.         if(rootfile==NULL) rootfile = new TFile(fnameroot);
  39.         if(rootfile==NULL) {
  40.           printf("Cannot open root file %s!!!\n",fnameroot);
  41.           return(0);
  42.         }
  43.  
  44.         // set draw style
  45.         gStyle->SetOptStat("ne");
  46.         gStyle->SetPalette(1, 0);
  47.        
  48.         gStyle->SetPaperSize(TStyle::kA4);
  49.         gStyle->SetStatBorderSize(1);
  50.         gStyle->SetFrameBorderMode(0);
  51.         gStyle->SetFrameFillColor(0);
  52.         gStyle->SetCanvasBorderMode(0);
  53.         gStyle->SetPadBorderMode(0);
  54.         gStyle->SetPadColor(0);
  55.         gStyle->SetCanvasColor(0);
  56.         gStyle->SetStatColor(0);
  57.         gStyle->SetOptFit(11);
  58.         gStyle->SetOptStat();
  59.         gStyle->SetPadRightMargin(0.15);
  60.         gStyle->SetPadLeftMargin(0.12);
  61.         //gStyle->SetTitleYOffset(1.4);
  62.  
  63.   if (strstr(plopt, "th") != NULL) {
  64.    TCanvas* canvas1 = new TCanvas("canvas1","",500,500);
  65.    canvas1->cd();
  66.    TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
  67.    h_threshold->Draw("colz");
  68.   }
  69.  
  70.   if( strstr(plopt, "all") != NULL ) {
  71.     TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);
  72.     TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);
  73.     canvas2->Divide(8,8);
  74.     canvas3->Divide(8,8);
  75.     TH1F* h_hitsx;
  76.     TH1F* h_hitsy;
  77.     for(int i=0; i<c_nChannels; i++) {
  78.       canvas2->cd(i+1);
  79.       char hname[128];
  80.       sprintf(hname, "hnhitsx%d", i);
  81.       h_hitsx = (TH1F*)rootfile->Get(hname);
  82.       h_hitsx->Draw();
  83.       canvas3->cd(i+1);
  84.       sprintf(hname, "hnhitsy%d", i);
  85.       h_hitsy = (TH1F*)rootfile->Get(hname);
  86.       h_hitsy->Draw();  
  87.     }
  88.   }
  89.  
  90.   if( strstr(plopt, "x") != NULL ) {
  91.     TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);
  92.     TH1F* h_hitsx;
  93.     canvas10->cd();
  94.     char hname[128];
  95.     sprintf(hname, "hnhitsx%d", chXstart);
  96.     h_hitsx = (TH1F*)rootfile->Get(hname);
  97.     h_hitsx->Draw();
  98.   }
  99.  
  100.   if( strstr(plopt, "y") != NULL ) {
  101.     TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);
  102.     TH1F* h_hitsy;
  103.     canvas11->cd();
  104.     char hname[128];
  105.     sprintf(hname, "hnhitsy%d", chXstart);
  106.     h_hitsy = (TH1F*)rootfile->Get(hname);
  107.     h_hitsy->Draw();
  108.   }
  109.  
  110.   if( strstr(plopt, "share") != NULL ) {
  111.    /*TCanvas *canvas4 = new TCanvas("canvas1","canvas1",1000,1000);
  112.    int nChannels = chYend-chYstart+1;
  113.    int ncols = nChannels/2;
  114.    printf("nch %d nch\\2 %d\n", nChannels, ncols);
  115.    canvas4->Divide(2,ncols);
  116.    TH1F* h_hitsy;
  117.    for(int i=chYstart; i<=chYend; i++){
  118.      canvas4->cd(i-chYstart+1);
  119.      char hname[128];
  120.      int chPosition = map[0][i];
  121.      sprintf(hname, "hnhitsy%d", chPosition);
  122.      h_hitsy = (TH1F*)rootfile->Get(hname);
  123.      h_hitsy->Draw();
  124.    }*/
  125.    
  126.    TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);
  127.    canvas4->cd();
  128.    for(int i=chXstart; i<=chXend; i++) {
  129.      TH1F* h_hitsx;
  130.      char hname[128];
  131.      int chPosition = map[i][chYstart];
  132.      sprintf(hname, "hnhitsx%d", chPosition);
  133.      h_hitsx = (TH1F*)rootfile->Get(hname);
  134.      h_hitsx->SetTitle("Scan X;x [mm]; Entries");
  135.      h_hitsx->GetYaxis()->SetTitleOffset(1.3);
  136.      h_hitsx->SetStats(0);
  137.      if (i == chXstart)
  138.       h_hitsx->Draw();
  139.      else {
  140.       h_hitsx->SetLineColor(i+1);
  141.       h_hitsx->Draw("same");
  142.      }
  143.    }
  144.    //sprintf(fullname, "ps/%s_Yshare.eps", filename);
  145.          //canvas4->SaveAs(fullname);
  146.    
  147.    TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);
  148.    canvas5->cd();
  149.    for(int i=chYstart; i<=chYend; i++) {
  150.      TH1F* h_hitsy;
  151.      char hname[128];
  152.      int chPosition = map[chXstart][i];
  153.      sprintf(hname, "hnhitsy%d", chPosition);
  154.      h_hitsy = (TH1F*)rootfile->Get(hname);
  155.      h_hitsy->SetTitle("Scan Y;y [mm]; Entries");
  156.      h_hitsy->GetYaxis()->SetTitleOffset(1.3);
  157.      h_hitsy->SetStats(0);
  158.      if (i == chYstart)
  159.       h_hitsy->Draw();
  160.      else {
  161.       h_hitsy->SetLineColor(i+1);
  162.       h_hitsy->Draw("same");
  163.      }
  164.    }
  165.    //sprintf(fullname, "ps/%s_Yshare.eps", filename);
  166.          //canvas5->SaveAs(fullname);
  167.   }
  168.  
  169.   /** Draws the signal from 8 channels in x-row
  170.     * for one specific y bin, so the background and cross-talk
  171.     * can be estimated.
  172.     * Draws also a 2d scan of these channels.
  173.     */
  174.   if (strstr(plopt, "line") != NULL) {
  175.     TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);
  176.     canvas6->cd(0);
  177.     gStyle->SetOptStat(0);
  178.    
  179.     TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
  180.           Int_t binsX = h0->GetXaxis()->GetNbins();
  181.     Int_t minX  = h0->GetXaxis()->GetFirst();
  182.     Int_t maxX  = h0->GetXaxis()->GetLast()+1;
  183.     Int_t binsY = h0->GetYaxis()->GetNbins();
  184.     Int_t minY  = h0->GetYaxis()->GetFirst();
  185.     Int_t maxY  = h0->GetYaxis()->GetLast()+1;
  186.     Double_t xLowUser  = h0->GetXaxis()->GetBinLowEdge(minX);
  187.     Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
  188.     Double_t yLowUser  = h0->GetYaxis()->GetBinLowEdge(minY);
  189.     Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
  190.     if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
  191.    
  192.     //! 1-dimension position in x vs. hits
  193.     TH2F* h[8];
  194.     TH1F* h_line[8];
  195.     for(int j=0; j<8; j++) {
  196.       h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);
  197.     }
  198.    
  199.     for(int j=chXstart; j<=chXend; j++) {
  200.       int chPosition = map[j][chYstart];
  201.       char hname[128];
  202.       sprintf(hname, "h2d%d", chPosition);  
  203.       int histogram = j;
  204.       h[histogram] = (TH2F *) rootfile->Get(hname);
  205.       int noise = getNoise(h[histogram], 1, 160);
  206.       for(int k=minX; k<=maxX; k++) {
  207.         int l=chYstart*20+12;
  208.         //for(int l=12; l<=16; l++) {
  209.           double signal = h[histogram]->GetBinContent(k,l);
  210.           //signal -= noise;
  211.           //signal /= 5*10000.0;
  212.           double eta = -log(1 - signal);
  213.           double x = k*(xUpUser-xLowUser)/double(binsX);
  214.           //double y = l*(yUpUser-yLowUser)/double(binsY);
  215.           h_line[j]->Fill(x-c_xOffset, signal);
  216.         //}
  217.       }
  218.       if (j == chXstart) {
  219.         h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");
  220.         //h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);
  221.         //h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);
  222.         h_line[j]->Draw("");
  223.       }
  224.       else {
  225.         h_line[j]->SetLineColor(j+1);
  226.         h_line[j]->Draw("same");
  227.       }
  228.     }
  229.    
  230.    
  231.     //! 2d scan
  232.     TCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);
  233.     canvas61->cd();
  234.     TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
  235.     for(int i=chXstart; i<=chXend; i++) {
  236.       //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
  237.       //int canvasPosition = nX*(chYend-i)+chXstart+1;
  238.       //if (debug) printf("canvas %d\n",canvasPosition);
  239.         int chPosition = map[i][chYstart];
  240.         char hname[128];
  241.               sprintf(hname, "h2d%d", chPosition);  
  242.               int histogram = i;
  243.               h[histogram] = (TH2F *) rootfile->Get(hname);
  244.               int noise = getNoise(h[histogram], 1, 100);
  245.              
  246.               for(int k=minX; k<=maxX; k++) {
  247.           for(int l=minY; l<=maxY; l++) {
  248.             int signal = h[histogram]->GetBinContent(k,l); // detected
  249.             //p /= 10000.;
  250.             //double p0 = 1.0 - p; // events with zero photons
  251.             //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
  252.             //double eta = -log(p0);
  253.             //printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);      
  254.             //double signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;
  255.             double p = signal - noise;
  256.             p /= 10000.0;
  257.             double eta = -log(1 - p);
  258.             double x = k*(xUpUser-xLowUser)/double(binsX);
  259.             double y = l*(yUpUser-yLowUser)/double(binsY);
  260.             h_corrected->Fill(x-c_xOffset, y-c_yOffset, eta);
  261.           }
  262.         }
  263.     }
  264.     h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");
  265.     h_corrected->GetZaxis()->SetRangeUser(-0.05,0.30);
  266.     h_corrected->Draw("colz");
  267.    
  268.   }
  269.  
  270.   /** Draws the sum of the channels
  271.     * Each channel is a 2d plot
  272.     * Intended for the study of 1 channel
  273.     */
  274.   if (strstr(plopt, "2d") != NULL) {
  275.    
  276.     int nX = chXend - chXstart + 1;
  277.     int nY = chYend - chYstart + 1;
  278.     TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX*400,nY*400);
  279.     printf("nx %d ny %d\n",nX,nY);
  280.     canvas7->Divide(nX,nY);
  281.     for(int i=chYstart; i<=chYend; i++) {
  282.       for(int j=chXstart; j<=chXend; j++) {
  283.       //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
  284.       int canvasPosition = nX*(chYend-i)+(j-chXstart)+1;
  285.       if (debug) printf("canvas %d\n",canvasPosition);
  286.         canvas7->cd(canvasPosition);
  287.         char hname[128];
  288.         int chPosition = map[j][i];
  289.         sprintf(hname, "h2d%d", chPosition);
  290.         TH2F* h_2d = (TH2F*)rootfile->Get(hname);
  291.         h_2d->Draw("colz");
  292.       } //x
  293.     }
  294.    
  295.     // Number of photoelectrons - Poissonian correction
  296.     TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);
  297.     canvas8->cd();
  298.     gStyle->SetOptStat(0);
  299.       char hname[128];
  300.       int chPosition = map[1][2];
  301.       sprintf(hname, "h2d%d", chPosition);
  302.       TH2F* h_2d = (TH2F*)rootfile->Get(hname);
  303.  
  304.       Int_t binsX = h_2d->GetXaxis()->GetNbins();
  305.       Int_t minX  = h_2d->GetXaxis()->GetFirst();
  306.       Int_t maxX  = h_2d->GetXaxis()->GetLast()+1;
  307.       Int_t binsY = h_2d->GetYaxis()->GetNbins();
  308.       Int_t minY  = h_2d->GetYaxis()->GetFirst();
  309.       Int_t maxY  = h_2d->GetYaxis()->GetLast()+1;
  310.       Double_t xLowUser  = h_2d->GetXaxis()->GetBinLowEdge(minX);
  311.       Double_t xUpUser   = h_2d->GetXaxis()->GetBinUpEdge(maxX);
  312.       Double_t yLowUser  = h_2d->GetYaxis()->GetBinLowEdge(minY);
  313.       Double_t yUpUser =   h_2d->GetYaxis()->GetBinUpEdge(maxY);
  314.       if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
  315.       TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);
  316.      
  317.       double noise = getNoise(h_2d, 1, 89);
  318.       if(debug) printf("Noise = %f\n", noise);
  319.       for(int k=minX; k<=maxX; k++) {
  320.         for(int j=minY; j<=maxY; j++) {
  321.           double signal = h_2d->GetBinContent(k,j); // detected
  322.           //double p = ((signal - noise) > 1) ? (signal-noise) : 1;
  323.           double p = signal - noise;
  324.           p /= 10000.;
  325.           double p0 = 1.0 - p; // events with zero photons
  326.           //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
  327.           double eta = -log(p0); // constant of the poissonian statistics
  328.           if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
  329.           double x = k*(xUpUser-xLowUser)/double(binsX);
  330.           double y = j*(yUpUser-yLowUser)/double(binsY);
  331.           h_corrected->Fill(x+3,y+8, eta);
  332.         }
  333.        }
  334.       h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");
  335.       h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);
  336.       h_corrected->Draw("colz");
  337.      
  338.       // collection efficiency
  339.       int nPoints =0;
  340.       double efficiency=0;
  341.       for (int i=18; i<=58; i++) {
  342.         for (int j=19; j<=59; j++) {
  343.           double signal = h_corrected->GetBinContent(i,j);
  344.           if(debug) printf("signal %f\n",signal);
  345.           efficiency += signal;
  346.           nPoints++;
  347.         }
  348.       }  
  349.     printf("Signal sum = %f\n   # of points = %d\n",efficiency,nPoints);
  350.   }
  351.  
  352.   /** Draws the sum of channel signals
  353.     * Each channel is a 2d ('h2d') histogram
  354.     * Suitable for 8x8 chs scan
  355.     */
  356.   if( strstr(plopt, "sum") != NULL ) {
  357.     int nX = chXend - chXstart + 1;
  358.     int nY = chYend - chYstart + 1;
  359.           TCanvas* canvas12 = new TCanvas("canvas12","c2",8*200, 8*200);
  360.           canvas12->cd();
  361.           gStyle->SetOptStat(0);
  362.  
  363.     // final histogram parameters
  364.     TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
  365.           Int_t binsX = h0->GetXaxis()->GetNbins();
  366.     Int_t minX  = h0->GetXaxis()->GetFirst();
  367.     Int_t maxX  = h0->GetXaxis()->GetLast()+1;
  368.     Int_t binsY = h0->GetYaxis()->GetNbins();
  369.     Int_t minY  = h0->GetYaxis()->GetFirst();
  370.     Int_t maxY  = h0->GetYaxis()->GetLast()+1;
  371.     Double_t xLowUser  = h0->GetXaxis()->GetBinLowEdge(minX);
  372.     Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
  373.     Double_t yLowUser  = h0->GetYaxis()->GetBinLowEdge(minY);
  374.     Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
  375.     if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
  376.     TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
  377.     TH2F* h[9];
  378.    
  379.     // 2d histogram noise subtraction and poisson scaling
  380.     for(int i=chYstart; i<=chYend; i++) {
  381.       for(int j=chXstart; j<=chXend; j++) {
  382.         int chPosition = map[j][i];
  383.         char hname[128];
  384.               sprintf(hname, "h2d%d", chPosition);  
  385.               int histogram = nX*(i-chYstart)+(j-chXstart);
  386.               h[histogram] = (TH2F *) rootfile->Get(hname);
  387.               int noise = getNoise(h[histogram], 1, 170);
  388.               if (debug) printf("noise: %d\n",noise);
  389.               for(int k=minX; k<=maxX; k++) {
  390.           for(int l=minY; l<=maxY; l++) {
  391.             int signal = h[histogram]->GetBinContent(k,l); // detected
  392.             //double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;
  393.             double p = signal - noise;
  394.             p /= 10000.;
  395.             double p0 = 1.0 - p; // events with zero photons
  396.             //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
  397.             double eta = -log(p0);
  398.             //printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
  399.             double x = k*(xUpUser-xLowUser)/double(binsX);
  400.             double y = l*(yUpUser-yLowUser)/double(binsY);
  401.             h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);
  402.           }
  403.         }    
  404.      
  405.       }
  406.     }
  407.  
  408.     h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");
  409.     //h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);
  410.     h_corrected->Draw("colz");
  411.    
  412.     TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);
  413.     canvas13->Divide(2);
  414.     canvas13->cd(1);
  415.     h[16]->Draw("colz");
  416.     canvas13->cd(2);
  417.     h[8]->Draw("colz");
  418.   }
  419.  
  420.   /** Draws the beam profile and fits it with error function
  421.     * on some background function
  422.     */
  423.   if (strstr(plopt, "beam") != NULL) {
  424.  
  425.   TCanvas* canvas9 = new TCanvas("canvas9","canvas9", 500,500);
  426.   canvas9->cd();
  427.   char hname[128];
  428.   sprintf(hname, "hnhitsx%d", 36);
  429.   TH1F* h_laser = (TH1F*)rootfile->Get(hname);
  430.   h_laser->Draw();
  431.   h_laser->SetStats(1);
  432.  
  433.   TF1* err = new TF1("err","[0]+[1]*TMath::Erf((x-[2])/[3])",17.2,17.50);
  434.   err->SetParameter(0,2500);
  435.   err->SetParameter(1, h_laser->GetMaximum());
  436.   err->SetParameter(2, h_laser->GetBinCenter(h_laser->GetMaximumBin()));
  437.   err->SetParameter(3, 0.001);
  438.   h_laser->Fit(err,"qr");
  439.   h_laser->Fit(err,"lr");
  440.   double sigma = err->GetParameter(3);
  441.   printf("sigma = %2.0f um,     FWHM = %2.0f um\n", sigma*1000, 2.35*sigma*1000);
  442.   }
  443.  
  444.   if (strstr(plopt, "map") != NULL) {
  445.     for (int i=7; i>=0; i--) {
  446.       //for (int j=7; j>=0; j--) printf("(%d, %d) ", j,i);
  447.       for (int j=7; j>=0; j--) printf("%2d ", map[j][i]);
  448.       printf("\n");
  449.     }
  450.      
  451.   }
  452.    
  453.   return(0);
  454. }
  455.  
  456. /** Function calculates the noise from one channel
  457.   * it runs through the bins along x and returns the average value
  458.   */
  459. double getNoise(TH2F* histogram, int yStart, int yEnd)
  460. {
  461.   double noise=0;
  462.   int count=0;
  463.   for(int j=yStart; j<yEnd; j++) {
  464.     double value = histogram->GetBinContent(j,2);
  465.     //if (noise < value) noise = value;
  466.     noise += value;
  467.     count++;
  468.     }
  469.   return (noise/double(count));
  470. }
  471.  
  472.