Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3.  
  4. #include <TROOT.h>
  5. #include <TH1D.h>
  6. #include <TH1F.h>
  7. #include <TH2D.h>
  8. #include <TH2F.h>
  9. #include <TH2I.h>
  10. #include <TCanvas.h>
  11. #include <TStyle.h>
  12. #include <TSystem.h>
  13. #include <TFile.h>
  14. #include <TDirectory.h>
  15. #include <TPaveText.h>
  16. #include <TSpectrum.h>
  17. #include <TF1.h>
  18. #include <algorithm>
  19.  
  20. #include "base.h"
  21. #include "fit.h"
  22.  
  23. int fit(int kanal) {
  24.   GRID m = mapping();
  25.   char name[128];
  26.   sprintf(name,"hcharge0_0;%d",m.pozicijaPixla[kanal]+1);
  27.  
  28.   TF1 * g0,  * g1, * g2,  * g3;
  29.   float gain;
  30.   float mean1, sigma1, mean2;
  31.  
  32.   TH1D * h = (TH1D *) gDirectory->Get(name);
  33.   h->GetXaxis()->SetRangeUser(-50,300);
  34.   h->Draw();
  35.  
  36.   TSpectrum * s = new TSpectrum(4);
  37.   int nfound = s->Search(h);
  38.  
  39.   float * xpeaks = s->GetPositionX();
  40.  
  41.   //SORTIRANJE
  42.   printf("Peaks found:\t");
  43.   for(int i=0; i<nfound;i++){
  44.     printf("%f\t",xpeaks[i]);
  45.     if(i==nfound-1) printf("\n");
  46.   }    
  47.   float* first(&xpeaks[0]);
  48.   float* last(first + nfound);
  49.   std::sort(first, last);
  50.   printf("Sorted:\t\t");
  51.   for(int i=0; i<nfound;i++){
  52.     printf("%f\t",xpeaks[i]);
  53.     if(i==nfound-1) printf("\n");
  54.   }
  55.  
  56.   for (int p=0;p<nfound;p++) {
  57.     float xp = xpeaks[p];
  58.     printf("Vrh številka %d je na poziciji %f\n",p,xp);
  59.     if(p==0) {
  60.       g0 = new TF1("0ph","gaus",xp-14,xp+14);
  61.       h->Fit(g0,"QR");
  62.       //g0->GetParameters(&par[0]);
  63.     } else if(p==1) {
  64.       g1 = new TF1("1ph","gaus",xp-14,xp+14);
  65.       h->Fit(g1,"QR+");
  66.       //g1->GetParameters(&par[3]);
  67.      } else if(p==2) {
  68.       g2 = new TF1("2ph","gaus",xp-14,xp+14);
  69.       h->Fit(g2,"QR+");
  70.       //g2->GetParameters(&par[6]);
  71.     } else {
  72.       g3 = new TF1("3ph","gaus",xp-14,xp+14);
  73.       h->Fit(g3,"QR+");
  74.       //g3->GetParameters(&par[9]);
  75.     }
  76.   }
  77.  
  78.   if(nfound >= 1){
  79.     mean1 = g0->GetParameter(1);
  80.     sigma1 = g0->GetParameter(2);
  81.     if (sigma1<7) gain = 0;
  82.     else gain = 30;
  83.   } else {
  84.     gain = 0;
  85.   }
  86.   if(nfound > 1){
  87.     mean2 = g1->GetParameter(1);
  88.     printf("%f \t %f\n", mean1, mean2);
  89.     if (mean1<mean2) gain = mean2 - mean1;
  90.     else gain = mean1 - mean2;
  91.     printf("nfound: %d \t gain: %.1f \t mean1: %.1f \t sigma1: %.1f \t mean2: %.1f\n", nfound, gain, mean1, sigma1, mean2);
  92.   }
  93.  
  94.   delete s;
  95.   return gain;
  96.  
  97. }
  98.  
  99.   /*
  100.   for (int p=0;p<nfound;p++) {
  101.     float xp = xpeaks[p];
  102.     printf("Vrh številka %d je na poziciji %f\n",p,xp);
  103.     if(p==0) {
  104.       TF1 * g0 = new TF1("0ph","gaus",xp-14,xp+14);
  105.       h->Fit(g0,"QR");
  106.       g0->GetParameters(&par[0]);
  107.       delete g0;
  108.     } else if(p==1) {
  109.       TF1 * g1 = new TF1("1ph","gaus",xp-14,xp+14);
  110.       h->Fit(g1,"QR+");
  111.       g1->GetParameters(&par[3]);
  112.       delete g1;
  113.      } else if(p==2) {
  114.       TF1 * g2 = new TF1("2ph","gaus",xp-14,xp+14);
  115.       h->Fit(g2,"QR+");
  116.       g2->GetParameters(&par[6]);
  117.       delete g2;
  118.     } else {
  119.       TF1 * g3 = new TF1("3ph","gaus",xp-14,xp+14);
  120.       h->Fit(g3,"QR+");
  121.       g3->GetParameters(&par[9]);
  122.       delete g3;
  123.     }
  124.   }
  125.   if(nfound==1) total = new TF1("mstotal","gaus(0)",xpeaks[0]-14,xpeaks[0]+14);
  126.   else if(nfound==2) total = new TF1("mstotal","gaus(0)+gaus(3)",xpeaks[0]-14,xpeaks[1]+14);
  127.   else if(nfound==3) total = new TF1("mstotal","gaus(0)+gaus(3)+gaus(6)",xpeaks[0]-14,xpeaks[2]+14);
  128.   else total = new TF1("mstotal","gaus(0)+gaus(3)+gaus(6)+gaus(9)",xpeaks[0]-14,xpeaks[3]+14);
  129.   total->SetParameters(par);
  130.   h->Fit(total,"QR0+");
  131.  
  132.   float mean1 = total->GetParameter(1);
  133.   //float error1 = total->GetParError(1);
  134.   float sigma1 = total->GetParameter(2);
  135.  
  136.   float mean2 = total->GetParameter(4);
  137.   //float error2 = total->GetParError(4);
  138.   */
  139.  
  140.  
  141. /*
  142. int test(int channel) {
  143.   char name[128];
  144.   sprintf(name,"hcharge2_0;%d",channel);
  145.   TH1D *h=(TH1D *) gDirectory->Get(name);
  146.   //h->Draw("");
  147.   double par[12];
  148.   TSpectrum *s = new TSpectrum(4);
  149.   int nfound = s->Search(h,10);
  150.   //printf("Found %d candidate peaks to fit\n",nfound);
  151.   float *xpeaks = s->GetPositionX();
  152.  
  153.   TF1 *g0 = new TF1("0ph","gaus",xpeaks[0]-3,xpeaks[0]+3);
  154.   h->Fit(g0,"QR");
  155.   g0->GetParameters(&par[0]);
  156.  
  157.   TF1 *g1 = new TF1("1ph","gaus",xpeaks[1]-3,xpeaks[1]+3);
  158.   h->Fit(g1,"QR+");
  159.   g1->GetParameters(&par[3]);
  160.  
  161.  
  162.  
  163.   //Float_t gainErr = sRt((error1/mean1)*(error1/mean1)+(error2/mean2)*(error2/mean2))
  164.   delete s;
  165.   return nfound;
  166. }*/