Subversion Repositories f9daq

Rev

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