Subversion Repositories f9daq

Rev

Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

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