Rev 1 |
Blame |
Compare with Previous |
Last modification |
View Log
| RSS feed
/*
Skripta za kalibracijo pozicije kanalov
Vhod: histogrami utezenih pozicij
Navodilo:
1.klikni na GCUT toolbar
2.poklikaj vrhove vseh kristalov po vrstnem redu zgoraj levo->zgoraj desno ->spodaj desno
3. zakljuci na zadnjem vrhu tako, da se enkrat kliknes na isto tocko
Izhod: Datoteka s kalibracijskimi histogrami, kjer je st. entrijev enako stevilki kanala
.x calibration.cxx("/tmp/pet.root","calibration.root")
*/
int calibration(char *fname, char* fout){
char cmd[256];
TFile *f= new TFile(fname);
TFile *fnew= new TFile(fout,"RECREATE");
for (int ii=0;ii<4;ii++){
sprintf(cmd,"c%d",ii);
TCanvas *c= new TCanvas(cmd,cmd,1500,1000);
sprintf(cmd,"pmt1%d",ii);
TH2 *h = (TH2 *) f->Get(cmd);
if (!h) {
cout << "Histogram not found " << cmd << endl;
continue;
} else {
}
c->ToggleToolBar();
c->Divide(2,1);
c->cd(1);
h->Draw("colz");
sprintf(cmd,"pmt1%d_calib",ii);
TH2F *h1= h->Clone(cmd);
h1->Reset();
c->WaitPrimitive("CUTG");
TCutG *mycutg = (TCutG*)gROOT->GetListOfSpecials()->FindObject("CUTG");
if (mycutg) {
sprintf(cmd,"cutg%d",ii);
mycutg->SetName(cmd);
mycutg->Print();
int n= mycutg->GetN();
double x,y;
TAxis *xaxis = h1->GetXaxis();
TAxis *yaxis = h1->GetYaxis();
for (int i=0; i < xaxis->GetNbins();i++){
double x0=yaxis->GetBinCenter(i+1);
for (int j=0; j < yaxis->GetNbins();j++){
double y0=yaxis->GetBinCenter(j+1);
double rmin=1e10;
int w=-1;
for (int k=0;k<n;k++){
mycutg->GetPoint(k, x, y);
double r2 = (x0-x)*(x0-x)+(y0-y)*(y0-y);
if (r2<rmin || k==0 ) {
w=k;
rmin=r2;
}
}
h1->Fill(x0,y0,w);
}
}
c->cd(2);
h1->Draw("colz");
h->Write();
mycutg->Write();
}
}
fnew->Write();
return 0;
}