Subversion Repositories f9daq

Rev

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