Subversion Repositories f9daq

Rev

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