WebSVN - f9daq - Rev 326 - /cvi/apps/sa02

#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;
}

Subversion Repositories f9daq

(root)/cvi/apps/sa02_daq/analysis/projection.c @ 354 - Rev 326