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