Subversion Repositories f9daq

#include <stdio.h>

#include <stdlib.h>

#include <TROOT.h>

#include <TH1D.h>

#include <TH2D.h>

#include <TH3D.h>

#include <TCanvas.h>

#include <TStyle.h>

#include <TSystem.h>

#include <TFile.h>

#include <TDirectory.h>

#include <TPaveText.h>

#include <iostream>

#include <string>

#include "base.h"

#include "projection.h"

#include "savetoroot.h"

int projection(int runNumber, int save) {

  int HAPDnumber; //Number of HAPD: 0, 1, 2 or 3!

  int direction, nx, ny;

  char name[128];

  char pdfname[128];

  char buf[256];

  float pNoise;

  int color[6] = {1,2,4,6,8,9};

  char hname[0xF];

  const char * serialNumber = getSN(runNumber);

  std::string serialNumberTemp = serialNumber;

  std::string HAPDserialNumber,FEBserialNumber;

  GRID m = mapping();

  int vrstaPisave = 82;

  TCanvas * c, * c2;

  TPad * pad1, * pad2, * pad3, * pad4, * c2pad1, * c2pad2;

  TPad * VirtualPad2[12];

  TLine * l;

  TPaveText * info, * axisName;

  TAxis * axis;

  TH3D * h;

  TH2D * slice[12], * sum; //take the right slice from 3D histogram

  TH1D * p[12], * zamenjan[12];

  /********** GLOBALNE **********/

    gStyle->SetOptStat(0);

    //gStyle->SetOptTitle(1);

    gStyle->SetTitleFontSize(.15);

  TDirectory *CurrentDirectory = gDirectory->CurrentDirectory();

  for(HAPDnumber=0;HAPDnumber<4;HAPDnumber++){

    // Change code here for runs where one or more modules were disabled during measurement

    //if(HAPDnumber==0) continue;

    //if(HAPDnumber==1) continue;

    //if(HAPDnumber==2) continue;

    //if(HAPDnumber==3) continue;

    SaveToRootInit();

    //std::cout << serialNumberTemp << std::endl;

    FEBserialNumber = serialNumberTemp.substr(0,serialNumberTemp.find(","));

    serialNumberTemp.erase(0,FEBserialNumber.length()+1);

    if(!FEBserialNumber.compare("noserial")) continue;

    HAPDserialNumber = FEBserialNumber.substr(0,FEBserialNumber.find("/"));

    FEBserialNumber.erase(0,FEBserialNumber.find("/")+1);

    sprintf(hname,"hxy%d_0;1",HAPDnumber);

    printf("%s \t %s \t %s \n",FEBserialNumber.c_str(),HAPDserialNumber.c_str(),hname);

    //~ h = (TH3D *) gDirectory->Get(hname);

    h = (TH3D *) CurrentDirectory->Get(hname);

    /********** SMER MERITVE **********/

    h->GetZaxis()->SetRange(1,1);

    nx = h->GetNbinsX();

    ny = h->GetNbinsY();

    if (nx>=ny) direction = 0;

    else direction = 1;

    printf("Direction = %d\n", direction);

    /********** CANVAS **********/

    if (!direction) {

      sprintf(buf,"%04d_%s_px",runNumber,HAPDserialNumber.c_str());

      c = new TCanvas(buf,buf,0,0,1200,1200);

    } else {

      sprintf(buf,"%04d_%s_py",runNumber,HAPDserialNumber.c_str());

      c = new TCanvas(buf,buf,0,0,1200,1200);

    }

    pad1 = new TPad("pad1","Title",0.05,0.965,0.95,0.99,0,0);

      pad1->SetFillStyle(4000);

      pad1->Draw();

    pad2 = new TPad("pad2","Graphs",0.05,.05,.95,0.97,0,0);

      pad2->SetFillStyle(4000);

      pad2->Draw();

    pad3 = new TPad("pad3","AxisTitle",0.57,0.03,0.92,0.0470,0,0);

      pad3->SetFillStyle(4000);

      pad3->Draw();

    pad4 = new TPad("pad4","Lines",0.07,0.05,1,0.97,0,0);

      pad4->SetFillStyle(4000);

      pad4->Draw();

      pad4->cd();

      l = new TLine();

      l->SetLineColor(kBlue);

      l->DrawLineNDC(0,.5003,.085,.5003);

      l->DrawLineNDC(0.84,.5003,.9,.5003);

      l->DrawLineNDC(0.452,0,.452,0.992);

    /********** PAD 1 **********/

    pad1->cd();

      info = new TPaveText(0.2,0,0.8,1,"ndc");

      info->SetBorderSize(0);

      info->SetFillColor(4000);

      info->SetTextSize(0.8);

      info->SetTextAlign(22);

      sprintf(name,"%s,  Nx: %d,  Ny: %d", serialNumber, nx, ny);

      info->AddText(name);

      //info->Draw();

    /********** PAD3 **********/

    pad3->cd();

      axisName = new TPaveText(0.2,0,0.8,1,"ndc");

      axisName->SetBorderSize(0);

      axisName->SetFillColor(4000);

      axisName->SetTextSize(.7);

      axisName->SetTextFont(vrstaPisave);

      axisName->SetTextAlign(22);

      sprintf(name,"Incident light position (mm)");

      axisName->AddText(name);

      axisName->Draw();

    /********** PAD2 **********/

    pad2->cd();

    pad2->Divide(1,12,0,0);

    for (int i=0;i<12;i++) {

      if (!direction) {

        pad2->cd(11 - i+1);

        VirtualPad2[11 - i] = (TPad *)(pad2->cd(11 - i+1)); //Tole nastaviš zato, da lahko daš Log skalo

        VirtualPad2[11 - i]->SetLogy(1);

        //VirtualPad2[11 - i]->SetLeftMargin(.05);

        //VirtualPad2[11 - i]->SetRightMargin(.05);

        //VirtualPad2[11 - i]->SetGrid();

      } else {

        pad2->cd(i+1);

        VirtualPad2[i] = (TPad *)(pad2->cd(11 - i+1)); //Tole nastaviš zato, da lahko daš Log skalo

        VirtualPad2[i]->SetLogy(1);

        //VirtualPad2[i]->SetLeftMargin(.05);

        //VirtualPad2[i]->SetRightMargin(.05);

        //VirtualPad2[i]->SetGrid();

      }

      for (int j=0;j<12;j++) {

        if (!direction) h->GetZaxis()->SetRange(m.koordinatniSistem[j][i]+1,m.koordinatniSistem[j][i]+1);

        else h->GetZaxis()->SetRange(m.koordinatniSistem[i][j]+1,m.koordinatniSistem[i][j]+1);

        slice[j] = (TH2D *)h->Project3D("pyx");

        if (!direction) {

          sprintf(name,"PX: (%d,%d)",j,i);

          p[j] = slice[j]->ProjectionX(name,i+1,i+1);

          sprintf(name,"Projection X of row %d",i);

        } else {

          sprintf(name,"PY: (%d,%d)",i,j);

          p[j] = slice[j]->ProjectionY(name,i+1,i+1);

          sprintf(name,"Projection Y of column %d",i);

        }

      axis = p[j]->GetXaxis();

      if(!direction) zamenjan[j] = new TH1D(name,p[j]->GetTitle(),axis->GetNbins(),-32.9,33.6);

      else zamenjan[j] = new TH1D(name,p[j]->GetTitle(),axis->GetNbins(),-33.97,32.55);

      for (int k = 0; k < axis->GetNbins(); k++) zamenjan[j]->SetBinContent(k+1,p[j]->GetBinContent(k+1));

      if(!direction) VirtualPad2[11]->SetBottomMargin(.17);

      else VirtualPad2[0]->SetBottomMargin(.17);

      zamenjan[j]->SetLineWidth(2);

      zamenjan[j]->SetTitle(name);

      zamenjan[j]->SetTitleOffset(1.5);

      zamenjan[j]->GetYaxis()->SetRangeUser(1,10000);

      zamenjan[j]->GetYaxis()->SetTickLength(0.01);

      zamenjan[j]->GetYaxis()->SetTitleFont(vrstaPisave);

      zamenjan[j]->GetYaxis()->SetLabelFont(vrstaPisave);

      zamenjan[j]->GetYaxis()->SetLabelSize(0.05);

      zamenjan[j]->GetYaxis()->SetLabelColor(kBlack);

      zamenjan[j]->GetYaxis()->SetTitle("Number of Events");

      zamenjan[j]->GetYaxis()->SetNdivisions(0);

      zamenjan[j]->GetYaxis()->CenterTitle();

      zamenjan[j]->GetYaxis()->SetTitleSize(0.05);

      zamenjan[j]->GetYaxis()->SetTitleOffset(.5);

      zamenjan[j]->GetYaxis()->SetTitleColor(kWhite);

      zamenjan[j]->GetXaxis()->SetTickLength(0.1);

      zamenjan[j]->GetXaxis()->SetTitleFont(vrstaPisave);

      zamenjan[j]->GetXaxis()->SetLabelFont(vrstaPisave);

      zamenjan[j]->GetXaxis()->SetLabelSize(0.15);

      zamenjan[j]->GetXaxis()->SetLabelColor(kBlack);

      zamenjan[j]->GetXaxis()->SetLabelOffset(.05);

      zamenjan[j]->GetXaxis()->SetTitle("Incident light position (mm)");

      //zamenjan[j]->GetXaxis()->SetNdivisions(0);

      //zamenjan[j]->GetXaxis()->CenterTitle();

      zamenjan[j]->GetXaxis()->SetTitleSize(0.01);

      zamenjan[j]->GetXaxis()->SetTitleColor(kWhite); //zato ker ga ne moreš pozicionirat

      zamenjan[j]->GetXaxis()->SetTitleOffset(.5);

      pNoise = noise(zamenjan[j]);

      for (int k=1;k<=zamenjan[j]->GetSize();k++) zamenjan[j]->SetBinContent(k,zamenjan[j]->GetBinContent(k)-pNoise);

      if (j<6) zamenjan[j]->SetLineColor(kWhite+color[j]);

      else zamenjan[j]->SetLineColor(kWhite+color[j-6]);

      zamenjan[j]->SetLineWidth(1);

      if (j==0) zamenjan[j]->DrawCopy();

      else zamenjan[j]->DrawCopy("SAME");

         if (!direction) sprintf(text_for_STR, "Scan_1D_X_%s_i_%d_j_%d", HAPDserialNumber.c_str(),i,j);

         else            sprintf(text_for_STR, "Scan_1D_Y_%s_i_%d_j_%d", HAPDserialNumber.c_str(),i,j);

        SaveToRootAddTH1DCopy(zamenjan[j], text_for_STR);

      /*

        //Comment the "zamenjan[j]" if you want to se scale in position steps! 1000 steps = 0.3959mm

        p[j]->SetLineWidth(2);

        p[j]->SetTitle("");

        p[j]->SetTitleOffset(1.5);

        p[j]->GetYaxis()->SetRangeUser(1,10000);

        p[j]->GetYaxis()->SetTickLength(0.01);

        p[j]->GetYaxis()->SetTitleFont(vrstaPisave);

        p[j]->GetYaxis()->SetLabelFont(vrstaPisave);

        p[j]->GetYaxis()->SetLabelSize(0.05);

        p[j]->GetYaxis()->SetLabelColor(kBlack);

        p[j]->GetYaxis()->SetTitle("Stevilo dogodkov");

        p[j]->GetYaxis()->SetNdivisions(0);

        p[j]->GetYaxis()->CenterTitle();

        p[j]->GetYaxis()->SetTitleSize(0.05);

        p[j]->GetYaxis()->SetTitleOffset(.5);

        p[j]->GetYaxis()->SetTitleColor(kWhite);

        p[j]->GetXaxis()->SetTickLength(0.1);

        p[j]->GetXaxis()->SetTitleFont(vrstaPisave);

        p[j]->GetXaxis()->SetLabelFont(vrstaPisave);

        p[j]->GetXaxis()->SetLabelSize(0.15);

        p[j]->GetXaxis()->SetLabelColor(kBlack);

        p[j]->GetXaxis()->SetLabelOffset(.05);

        p[j]->GetXaxis()->SetTitle("Pozicija vpadne svetlobe (korak mizice)");

        //p[j]->GetXaxis()->SetNdivisions(0);

        //p[j]->GetXaxis()->CenterTitle();

        p[j]->GetXaxis()->SetTitleSize(0.01);

        p[j]->GetXaxis()->SetTitleColor(kWhite); //zato ker ga ne moreš pozicionirat

        p[j]->GetXaxis()->SetTitleOffset(.5);

        pNoise = noise(p[j]);

        for (int k=1;k<=p[j]->GetSize();k++) p[j]->SetBinContent(k,p[j]->GetBinContent(k)-pNoise);

        if (j<6) p[j]->SetLineColor(kWhite+color[j]);

        else p[j]->SetLineColor(kWhite+color[j-6]);

        p[j]->SetLineWidth(1);

        if (j==0) p[j]->DrawCopy();

        else p[j]->DrawCopy("SAME");

      */

      }

    }

    c->Modified();

    c->Update();

    sprintf(buf,"%04d_%s_linearScan",runNumber,HAPDserialNumber.c_str());

    c2 = new TCanvas(buf,buf,0,0,2000,2000);

    c2->cd();

    c2pad2 = new TPad("c2pad2","c2pad2",0.125,0,.875,0.75,0,0);

    c2pad2->Draw();

    c2pad2->cd();

      c->DrawClonePad();

    c2->cd();

    c2pad1 = new TPad("c2pad1","c2pad1",0.2,0.74,.8,1,0,0);

    c2pad1->Draw();

    c2pad1->cd();

      c2pad1->SetRightMargin(.15);

      c2pad1->SetLogz();

      gStyle->SetTitleFontSize(0.07);

      sprintf(buf,"hxy%d_sum_0",HAPDnumber);

      //~ sum = (TH2D *) gDirectory->Get(buf);

      sum = (TH2D *) CurrentDirectory->Get(buf);

      if(!direction) sprintf(buf,"Scan in X direction");

      else sprintf(buf,"Scan in Y direction");

      sum->SetTitle(buf);

      sum->GetZaxis()->SetRangeUser(1,10000);

      sum->Draw("COLZ");

         if (!direction) sprintf(text_for_STR, "Scan_2D_X_%s", HAPDserialNumber.c_str());

         else            sprintf(text_for_STR, "Scan_2D_Y_%s", HAPDserialNumber.c_str());

        SaveToRootAddTH2DCopy(sum, text_for_STR);

    if (save==1){

      switch (direction) {

        case 0: sprintf(pdfname,"../modules/%04d/%04d_%s_2_2DX.pdf",runNumber,runNumber,HAPDserialNumber.c_str());break;

        case 1: sprintf(pdfname,"../modules/%04d/%04d_%s_3_2DY.pdf",runNumber,runNumber,HAPDserialNumber.c_str());break;

      }

      c->SaveAs(pdfname,"pdf");

    } else if (save==2){

      sprintf(pdfname,"../modules/%04d/%04d_%s.pdf",runNumber,runNumber,HAPDserialNumber.c_str());

      /*switch (direction) {

        case 0: sprintf(pdfname,"../modules/%04d/%04d_%s.pdf",runNumber,runNumber,HAPDserialNumber.c_str());break;

        case 1: sprintf(pdfname,"../modules/%04d/%04d_%s.pdf)",runNumber,runNumber,HAPDserialNumber.c_str());break;

      }*/

      c2->SaveAs(pdfname,"pdf");

    }

    sprintf(text_for_STR,"../modules/%04d/%04d_%s_out.root",runNumber,runNumber,HAPDserialNumber.c_str());

    SaveToRootWrite(text_for_STR, HAPDserialNumber.c_str());

  }

  return 0;

}

float noise(TH1D * proj) {

  float nBin=proj->GetSize(); //ugotovi koliko binov je v histogramu

  float wholeAverage=0;   //sem da povprečje celotnega kanala

  float noise=0;    //sem da povprečn šum

  int j=0;

  for (int i=1;i<=nBin;i++) wholeAverage+=proj->GetBinContent(i)/nBin;

  for (int i=1;i<=nBin;i++) {

    if (proj->GetBinContent(i)<1.1*wholeAverage) {

      noise+=proj->GetBinContent(i);

      j++;

    }

  }

  if(j!=0) return noise/j;

  else return 0;

}