Subversion Repositories f9daq

Rev

Rev 1 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1 f9daq 1
/*
2
 
3
Skripta za kalibracijo pozicije kanalov
4
Vhod: histogrami utezenih pozicij
5
 
6
Navodilo:
7
1.klikni na GCUT toolbar
8
2.poklikaj vrhove vseh kristalov po vrstnem redu zgoraj levo->zgoraj desno ->spodaj desno
9
3. zakljuci na zadnjem vrhu tako, da se enkrat kliknes na isto tocko
10
 
11
Izhod: Datoteka s kalibracijskimi histogrami, kjer je st. entrijev enako stevilki kanala
12
 
13
 .x calibration.cxx("/tmp/pet.root","calibration.root")
14
 
15
*/
16
 
17
int calibration(char *fname, char* fout){
18
 
19
char cmd[256];
20
TFile *f= new TFile(fname);
21
TFile *fnew= new TFile(fout,"RECREATE");
22
for (int ii=0;ii<4;ii++){
23
  sprintf(cmd,"c%d",ii);
24
  TCanvas *c= new TCanvas(cmd,cmd,1500,1000);
25
  sprintf(cmd,"pmt1%d",ii);
26
  TH2 *h = (TH2 *) f->Get(cmd);
27
  if (!h) {
28
     cout << "Histogram not found " << cmd << endl;
29
     continue;
30
  } else {
31
 
32
  }
33
  c->ToggleToolBar();
34
  c->Divide(2,1);
35
  c->cd(1);
36
  h->Draw("colz");
37
  sprintf(cmd,"pmt1%d_calib",ii);
38
  TH2F *h1= h->Clone(cmd);
39
  h1->Reset();
40
 
41
  c->WaitPrimitive("CUTG");
42
  TCutG *mycutg  = (TCutG*)gROOT->GetListOfSpecials()->FindObject("CUTG");
43
  if (mycutg) {
44
   sprintf(cmd,"cutg%d",ii);
45
   mycutg->SetName(cmd);
46
   mycutg->Print();
47
   int n= mycutg->GetN();
48
   double x,y;
49
   TAxis *xaxis =  h1->GetXaxis();
50
   TAxis *yaxis =  h1->GetYaxis();
51
   for (int i=0; i < xaxis->GetNbins();i++){
52
     double x0=yaxis->GetBinCenter(i+1);
53
     for (int j=0; j < yaxis->GetNbins();j++){
54
       double y0=yaxis->GetBinCenter(j+1);
55
       double rmin=1e10;
56
       int w=-1;
57
       for (int k=0;k<n;k++){
58
         mycutg->GetPoint(k, x, y);
59
         double r2 = (x0-x)*(x0-x)+(y0-y)*(y0-y);
60
         if (r2<rmin || k==0 ) {
61
           w=k;
62
           rmin=r2;
63
         }
64
 
65
       }
66
       h1->Fill(x0,y0,w);
67
     }
68
   }
69
   c->cd(2);
70
   h1->Draw("colz");
71
   h->Write();
72
   mycutg->Write();
73
  }
74
}
75
fnew->Write();
76
return 0;
77
}