Subversion Repositories f9daq

Rev

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