Subversion Repositories f9daq

/*

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