Subversion Repositories f9daq

Rev

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

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