Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3.  
  4. #include <TROOT.h>
  5. #include <TH1D.h>
  6. #include <TH2D.h>
  7. #include <TH3D.h>
  8. #include <TCanvas.h>
  9. #include <TStyle.h>
  10. #include <TSystem.h>
  11. #include <TFile.h>
  12. #include <TDirectory.h>
  13. #include <TPaveText.h>
  14.  
  15. #include <iostream>
  16. #include <string>
  17.  
  18. #include "base.h"
  19. #include "projection.h"
  20. #include "savetoroot.h"
  21.  
  22. int projection(int runNumber, int save) {
  23.  
  24.   int HAPDnumber; //Number of HAPD: 0, 1, 2 or 3!
  25.  
  26.   int direction, nx, ny;
  27.   char name[128];
  28.   char pdfname[128];
  29.   char buf[256];
  30.   float pNoise;
  31.   int color[6] = {1,2,4,6,8,9};
  32.  
  33.   char hname[0xF];
  34.  
  35.   const char * serialNumber = getSN(runNumber);
  36.   std::string serialNumberTemp = serialNumber;
  37.   std::string HAPDserialNumber,FEBserialNumber;
  38.  
  39.   GRID m = mapping();
  40.   int vrstaPisave = 82;
  41.  
  42.   TCanvas * c, * c2;
  43.   TPad * pad1, * pad2, * pad3, * pad4, * c2pad1, * c2pad2;
  44.   TPad * VirtualPad2[12];
  45.   TLine * l;
  46.   TPaveText * info, * axisName;
  47.   TAxis * axis;
  48.   TH3D * h;
  49.   TH2D * slice[12], * sum; //take the right slice from 3D histogram
  50.   TH1D * p[12], * zamenjan[12];
  51.  
  52.   /********** GLOBALNE **********/
  53.     gStyle->SetOptStat(0);
  54.     //gStyle->SetOptTitle(1);
  55.     gStyle->SetTitleFontSize(.15);
  56.    
  57.   TDirectory *CurrentDirectory = gDirectory->CurrentDirectory();
  58.  
  59.   for(HAPDnumber=0;HAPDnumber<4;HAPDnumber++){
  60.  
  61.     // Change code here for runs where one or more modules were disabled during measurement
  62.     //if(HAPDnumber==0) continue;
  63.     //if(HAPDnumber==1) continue;
  64.     //if(HAPDnumber==2) continue;
  65.     //if(HAPDnumber==3) continue;
  66.  
  67.     SaveToRootInit();
  68.            
  69.     //std::cout << serialNumberTemp << std::endl;
  70.     FEBserialNumber = serialNumberTemp.substr(0,serialNumberTemp.find(","));
  71.     serialNumberTemp.erase(0,FEBserialNumber.length()+1);
  72.     if(!FEBserialNumber.compare("noserial")) continue;
  73.     HAPDserialNumber = FEBserialNumber.substr(0,FEBserialNumber.find("/"));
  74.     FEBserialNumber.erase(0,FEBserialNumber.find("/")+1);
  75.  
  76.     sprintf(hname,"hxy%d_0;1",HAPDnumber);
  77.     printf("%s \t %s \t %s \n",FEBserialNumber.c_str(),HAPDserialNumber.c_str(),hname);
  78.  
  79.     //~ h = (TH3D *) gDirectory->Get(hname);
  80.     h = (TH3D *) CurrentDirectory->Get(hname);
  81.  
  82.     /********** SMER MERITVE **********/
  83.     h->GetZaxis()->SetRange(1,1);
  84.     nx = h->GetNbinsX();
  85.     ny = h->GetNbinsY();
  86.     if (nx>=ny) direction = 0;
  87.     else direction = 1;
  88.  
  89.     printf("Direction = %d\n", direction);
  90.  
  91.     /********** CANVAS **********/
  92.     if (!direction) {
  93.       sprintf(buf,"%04d_%s_px",runNumber,HAPDserialNumber.c_str());
  94.       c = new TCanvas(buf,buf,0,0,1200,1200);
  95.     } else {
  96.       sprintf(buf,"%04d_%s_py",runNumber,HAPDserialNumber.c_str());
  97.       c = new TCanvas(buf,buf,0,0,1200,1200);
  98.     }
  99.  
  100.     pad1 = new TPad("pad1","Title",0.05,0.965,0.95,0.99,0,0);
  101.       pad1->SetFillStyle(4000);
  102.       pad1->Draw();
  103.  
  104.     pad2 = new TPad("pad2","Graphs",0.05,.05,.95,0.97,0,0);
  105.       pad2->SetFillStyle(4000);
  106.       pad2->Draw();
  107.  
  108.     pad3 = new TPad("pad3","AxisTitle",0.57,0.03,0.92,0.0470,0,0);
  109.       pad3->SetFillStyle(4000);
  110.       pad3->Draw();
  111.  
  112.    
  113.     pad4 = new TPad("pad4","Lines",0.07,0.05,1,0.97,0,0);
  114.       pad4->SetFillStyle(4000);
  115.       pad4->Draw();
  116.       pad4->cd();
  117.  
  118.       l = new TLine();
  119.       l->SetLineColor(kBlue);
  120.       l->DrawLineNDC(0,.5003,.085,.5003);
  121.       l->DrawLineNDC(0.84,.5003,.9,.5003);
  122.       l->DrawLineNDC(0.452,0,.452,0.992);
  123.    
  124.     /********** PAD 1 **********/
  125.     pad1->cd();
  126.       info = new TPaveText(0.2,0,0.8,1,"ndc");
  127.       info->SetBorderSize(0);
  128.       info->SetFillColor(4000);
  129.       info->SetTextSize(0.8);
  130.       info->SetTextAlign(22);
  131.    
  132.       sprintf(name,"%s,  Nx: %d,  Ny: %d", serialNumber, nx, ny);
  133.       info->AddText(name);
  134.       //info->Draw();
  135.  
  136.     /********** PAD3 **********/
  137.     pad3->cd();
  138.       axisName = new TPaveText(0.2,0,0.8,1,"ndc");
  139.       axisName->SetBorderSize(0);
  140.       axisName->SetFillColor(4000);
  141.       axisName->SetTextSize(.7);
  142.       axisName->SetTextFont(vrstaPisave);
  143.       axisName->SetTextAlign(22);
  144.      
  145.       sprintf(name,"Incident light position (mm)");
  146.       axisName->AddText(name);
  147.       axisName->Draw();
  148.    
  149.     /********** PAD2 **********/
  150.     pad2->cd();
  151.     pad2->Divide(1,12,0,0);
  152.  
  153.     for (int i=0;i<12;i++) {
  154.       if (!direction) {
  155.         pad2->cd(11 - i+1);
  156.         VirtualPad2[11 - i] = (TPad *)(pad2->cd(11 - i+1)); //Tole nastaviš zato, da lahko daš Log skalo
  157.         VirtualPad2[11 - i]->SetLogy(1);
  158.         //VirtualPad2[11 - i]->SetLeftMargin(.05);
  159.         //VirtualPad2[11 - i]->SetRightMargin(.05);
  160.         //VirtualPad2[11 - i]->SetGrid();
  161.       } else {
  162.         pad2->cd(i+1);
  163.         VirtualPad2[i] = (TPad *)(pad2->cd(11 - i+1)); //Tole nastaviš zato, da lahko daš Log skalo
  164.         VirtualPad2[i]->SetLogy(1);
  165.         //VirtualPad2[i]->SetLeftMargin(.05);
  166.         //VirtualPad2[i]->SetRightMargin(.05);
  167.         //VirtualPad2[i]->SetGrid();
  168.       }
  169.       for (int j=0;j<12;j++) {
  170.         if (!direction) h->GetZaxis()->SetRange(m.koordinatniSistem[j][i]+1,m.koordinatniSistem[j][i]+1);
  171.         else h->GetZaxis()->SetRange(m.koordinatniSistem[i][j]+1,m.koordinatniSistem[i][j]+1);
  172.  
  173.         slice[j] = (TH2D *)h->Project3D("pyx");
  174.  
  175.         if (!direction) {
  176.           sprintf(name,"PX: (%d,%d)",j,i);
  177.           p[j] = slice[j]->ProjectionX(name,i+1,i+1);
  178.           sprintf(name,"Projection X of row %d",i);
  179.         } else {
  180.           sprintf(name,"PY: (%d,%d)",i,j);
  181.           p[j] = slice[j]->ProjectionY(name,i+1,i+1);
  182.           sprintf(name,"Projection Y of column %d",i);
  183.         }
  184.       axis = p[j]->GetXaxis();
  185.  
  186.       if(!direction) zamenjan[j] = new TH1D(name,p[j]->GetTitle(),axis->GetNbins(),-32.9,33.6);
  187.       else zamenjan[j] = new TH1D(name,p[j]->GetTitle(),axis->GetNbins(),-33.97,32.55);
  188.  
  189.       for (int k = 0; k < axis->GetNbins(); k++) zamenjan[j]->SetBinContent(k+1,p[j]->GetBinContent(k+1));
  190.  
  191.       if(!direction) VirtualPad2[11]->SetBottomMargin(.17);
  192.       else VirtualPad2[0]->SetBottomMargin(.17);
  193.  
  194.       zamenjan[j]->SetLineWidth(2);
  195.  
  196.       zamenjan[j]->SetTitle(name);
  197.       zamenjan[j]->SetTitleOffset(1.5);
  198.  
  199.       zamenjan[j]->GetYaxis()->SetRangeUser(1,10000);
  200.       zamenjan[j]->GetYaxis()->SetTickLength(0.01);
  201.       zamenjan[j]->GetYaxis()->SetTitleFont(vrstaPisave);
  202.       zamenjan[j]->GetYaxis()->SetLabelFont(vrstaPisave);
  203.       zamenjan[j]->GetYaxis()->SetLabelSize(0.05);
  204.       zamenjan[j]->GetYaxis()->SetLabelColor(kBlack);
  205.       zamenjan[j]->GetYaxis()->SetTitle("Number of Events");
  206.       zamenjan[j]->GetYaxis()->SetNdivisions(0);
  207.       zamenjan[j]->GetYaxis()->CenterTitle();
  208.       zamenjan[j]->GetYaxis()->SetTitleSize(0.05);
  209.       zamenjan[j]->GetYaxis()->SetTitleOffset(.5);
  210.       zamenjan[j]->GetYaxis()->SetTitleColor(kWhite);
  211.  
  212.       zamenjan[j]->GetXaxis()->SetTickLength(0.1);
  213.       zamenjan[j]->GetXaxis()->SetTitleFont(vrstaPisave);
  214.       zamenjan[j]->GetXaxis()->SetLabelFont(vrstaPisave);
  215.       zamenjan[j]->GetXaxis()->SetLabelSize(0.15);
  216.       zamenjan[j]->GetXaxis()->SetLabelColor(kBlack);
  217.       zamenjan[j]->GetXaxis()->SetLabelOffset(.05);
  218.       zamenjan[j]->GetXaxis()->SetTitle("Incident light position (mm)");
  219.       //zamenjan[j]->GetXaxis()->SetNdivisions(0);
  220.       //zamenjan[j]->GetXaxis()->CenterTitle();
  221.       zamenjan[j]->GetXaxis()->SetTitleSize(0.01);
  222.       zamenjan[j]->GetXaxis()->SetTitleColor(kWhite); //zato ker ga ne moreš pozicionirat
  223.       zamenjan[j]->GetXaxis()->SetTitleOffset(.5);
  224.  
  225.       pNoise = noise(zamenjan[j]);
  226.       for (int k=1;k<=zamenjan[j]->GetSize();k++) zamenjan[j]->SetBinContent(k,zamenjan[j]->GetBinContent(k)-pNoise);
  227.  
  228.       if (j<6) zamenjan[j]->SetLineColor(kWhite+color[j]);
  229.       else zamenjan[j]->SetLineColor(kWhite+color[j-6]);
  230.       zamenjan[j]->SetLineWidth(1);
  231.       if (j==0) zamenjan[j]->DrawCopy();
  232.       else zamenjan[j]->DrawCopy("SAME");
  233.          
  234.          if (!direction) sprintf(text_for_STR, "Scan_1D_X_%s_i_%d_j_%d", HAPDserialNumber.c_str(),i,j);
  235.          else            sprintf(text_for_STR, "Scan_1D_Y_%s_i_%d_j_%d", HAPDserialNumber.c_str(),i,j);
  236.         SaveToRootAddTH1DCopy(zamenjan[j], text_for_STR);
  237.  
  238.       /*
  239.         //Comment the "zamenjan[j]" if you want to se scale in position steps! 1000 steps = 0.3959mm
  240.         p[j]->SetLineWidth(2);
  241.  
  242.         p[j]->SetTitle("");
  243.         p[j]->SetTitleOffset(1.5);
  244.  
  245.         p[j]->GetYaxis()->SetRangeUser(1,10000);
  246.         p[j]->GetYaxis()->SetTickLength(0.01);
  247.         p[j]->GetYaxis()->SetTitleFont(vrstaPisave);
  248.         p[j]->GetYaxis()->SetLabelFont(vrstaPisave);
  249.         p[j]->GetYaxis()->SetLabelSize(0.05);
  250.         p[j]->GetYaxis()->SetLabelColor(kBlack);
  251.         p[j]->GetYaxis()->SetTitle("Stevilo dogodkov");
  252.         p[j]->GetYaxis()->SetNdivisions(0);
  253.         p[j]->GetYaxis()->CenterTitle();
  254.         p[j]->GetYaxis()->SetTitleSize(0.05);
  255.         p[j]->GetYaxis()->SetTitleOffset(.5);
  256.         p[j]->GetYaxis()->SetTitleColor(kWhite);
  257.  
  258.         p[j]->GetXaxis()->SetTickLength(0.1);
  259.         p[j]->GetXaxis()->SetTitleFont(vrstaPisave);
  260.         p[j]->GetXaxis()->SetLabelFont(vrstaPisave);
  261.         p[j]->GetXaxis()->SetLabelSize(0.15);
  262.         p[j]->GetXaxis()->SetLabelColor(kBlack);
  263.         p[j]->GetXaxis()->SetLabelOffset(.05);
  264.         p[j]->GetXaxis()->SetTitle("Pozicija vpadne svetlobe (korak mizice)");
  265.         //p[j]->GetXaxis()->SetNdivisions(0);
  266.         //p[j]->GetXaxis()->CenterTitle();
  267.         p[j]->GetXaxis()->SetTitleSize(0.01);
  268.         p[j]->GetXaxis()->SetTitleColor(kWhite); //zato ker ga ne moreš pozicionirat
  269.         p[j]->GetXaxis()->SetTitleOffset(.5);
  270.  
  271.         pNoise = noise(p[j]);
  272.         for (int k=1;k<=p[j]->GetSize();k++) p[j]->SetBinContent(k,p[j]->GetBinContent(k)-pNoise);
  273.  
  274.         if (j<6) p[j]->SetLineColor(kWhite+color[j]);
  275.         else p[j]->SetLineColor(kWhite+color[j-6]);
  276.         p[j]->SetLineWidth(1);
  277.         if (j==0) p[j]->DrawCopy();
  278.         else p[j]->DrawCopy("SAME");
  279.       */
  280.       }
  281.     }
  282.     c->Modified();
  283.     c->Update();
  284.  
  285.     sprintf(buf,"%04d_%s_linearScan",runNumber,HAPDserialNumber.c_str());
  286.     c2 = new TCanvas(buf,buf,0,0,2000,2000);
  287.     c2->cd();
  288.     c2pad2 = new TPad("c2pad2","c2pad2",0.125,0,.875,0.75,0,0);
  289.     c2pad2->Draw();
  290.     c2pad2->cd();
  291.       c->DrawClonePad();
  292.     c2->cd();
  293.     c2pad1 = new TPad("c2pad1","c2pad1",0.2,0.74,.8,1,0,0);
  294.     c2pad1->Draw();
  295.     c2pad1->cd();
  296.       c2pad1->SetRightMargin(.15);
  297.       c2pad1->SetLogz();
  298.       gStyle->SetTitleFontSize(0.07);
  299.       sprintf(buf,"hxy%d_sum_0",HAPDnumber);
  300.       //~ sum = (TH2D *) gDirectory->Get(buf);
  301.       sum = (TH2D *) CurrentDirectory->Get(buf);
  302.       if(!direction) sprintf(buf,"Scan in X direction");
  303.       else sprintf(buf,"Scan in Y direction");
  304.       sum->SetTitle(buf);
  305.       sum->GetZaxis()->SetRangeUser(1,10000);
  306.       sum->Draw("COLZ");
  307.  
  308.        
  309.          if (!direction) sprintf(text_for_STR, "Scan_2D_X_%s", HAPDserialNumber.c_str());
  310.          else            sprintf(text_for_STR, "Scan_2D_Y_%s", HAPDserialNumber.c_str());
  311.         SaveToRootAddTH2DCopy(sum, text_for_STR);
  312.  
  313.     if (save==1){
  314.       switch (direction) {
  315.         case 0: sprintf(pdfname,"../modules/%04d/%04d_%s_2_2DX.pdf",runNumber,runNumber,HAPDserialNumber.c_str());break;
  316.         case 1: sprintf(pdfname,"../modules/%04d/%04d_%s_3_2DY.pdf",runNumber,runNumber,HAPDserialNumber.c_str());break;
  317.       }
  318.       c->SaveAs(pdfname,"pdf");
  319.     } else if (save==2){
  320.       sprintf(pdfname,"../modules/%04d/%04d_%s.pdf",runNumber,runNumber,HAPDserialNumber.c_str());
  321.       /*switch (direction) {
  322.         case 0: sprintf(pdfname,"../modules/%04d/%04d_%s.pdf",runNumber,runNumber,HAPDserialNumber.c_str());break;
  323.         case 1: sprintf(pdfname,"../modules/%04d/%04d_%s.pdf)",runNumber,runNumber,HAPDserialNumber.c_str());break;
  324.       }*/
  325.       c2->SaveAs(pdfname,"pdf");
  326.     }
  327.    
  328.     sprintf(text_for_STR,"../modules/%04d/%04d_%s_out.root",runNumber,runNumber,HAPDserialNumber.c_str());
  329.     SaveToRootWrite(text_for_STR, HAPDserialNumber.c_str());
  330.   }
  331.  
  332.   return 0;
  333. }
  334.  
  335. float noise(TH1D * proj) {
  336.   float nBin=proj->GetSize(); //ugotovi koliko binov je v histogramu
  337.   float wholeAverage=0;   //sem da povprečje celotnega kanala
  338.   float noise=0;    //sem da povprečn šum
  339.   int j=0;
  340.   for (int i=1;i<=nBin;i++) wholeAverage+=proj->GetBinContent(i)/nBin;
  341.   for (int i=1;i<=nBin;i++) {
  342.     if (proj->GetBinContent(i)<1.1*wholeAverage) {
  343.       noise+=proj->GetBinContent(i);
  344.       j++;
  345.     }
  346.   }
  347.   if(j!=0) return noise/j;
  348.   else return 0;
  349. }
  350.