Subversion Repositories f9daq

#include <stdio.h>
#include <stdlib.h>

#include <TROOT.h>
#include <TH1D.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TH2F.h>
#include <TH2I.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TSystem.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TPaveText.h>
#include <TSpectrum.h>
#include <TF1.h>
#include <algorithm>

#include "base.h"
#include "fit.h"

int fit(int kanal) {
  GRID m = mapping();
  char name[128];
  sprintf(name,"hcharge0_0;%d",m.pozicijaPixla[kanal]+1);

  TF1 * g0,  * g1, * g2,  * g3;
  float gain;
  float mean1, sigma1, mean2;

  TH1D * h = (TH1D *) gDirectory->Get(name);
  h->GetXaxis()->SetRangeUser(-50,300);
  h->Draw();
 
  TSpectrum * s = new TSpectrum(4);
  int nfound = s->Search(h);

  float * xpeaks = s->GetPositionX();

  //SORTIRANJE
  printf("Peaks found:\t");
  for(int i=0; i<nfound;i++){
    printf("%f\t",xpeaks[i]);
    if(i==nfound-1) printf("\n");
  }    
  float* first(&xpeaks[0]);
  float* last(first + nfound);
  std::sort(first, last);
  printf("Sorted:\t\t");
  for(int i=0; i<nfound;i++){
    printf("%f\t",xpeaks[i]);
    if(i==nfound-1) printf("\n");
  }

  for (int p=0;p<nfound;p++) {
    float xp = xpeaks[p];
    printf("Vrh številka %d je na poziciji %f\n",p,xp);
    if(p==0) {
      g0 = new TF1("0ph","gaus",xp-14,xp+14);
      h->Fit(g0,"QR");
      //g0->GetParameters(&par[0]);
    } else if(p==1) {
      g1 = new TF1("1ph","gaus",xp-14,xp+14);
      h->Fit(g1,"QR+");
      //g1->GetParameters(&par[3]);
     } else if(p==2) {
      g2 = new TF1("2ph","gaus",xp-14,xp+14);
      h->Fit(g2,"QR+");
      //g2->GetParameters(&par[6]);
    } else {
      g3 = new TF1("3ph","gaus",xp-14,xp+14);
      h->Fit(g3,"QR+");
      //g3->GetParameters(&par[9]);
    }
  }

  if(nfound >= 1){
    mean1 = g0->GetParameter(1);
    sigma1 = g0->GetParameter(2);
    if (sigma1<7) gain = 0;
    else gain = 30;
  } else {
    gain = 0;
  }
  if(nfound > 1){
    mean2 = g1->GetParameter(1);
    printf("%f \t %f\n", mean1, mean2);
    if (mean1<mean2) gain = mean2 - mean1;
    else gain = mean1 - mean2;
    printf("nfound: %d \t gain: %.1f \t mean1: %.1f \t sigma1: %.1f \t mean2: %.1f\n", nfound, gain, mean1, sigma1, mean2);
  }

  delete s;
  return gain;

}

  /*
  for (int p=0;p<nfound;p++) {
    float xp = xpeaks[p];
    printf("Vrh številka %d je na poziciji %f\n",p,xp);
    if(p==0) {
      TF1 * g0 = new TF1("0ph","gaus",xp-14,xp+14);
      h->Fit(g0,"QR");
      g0->GetParameters(&par[0]);
      delete g0;
    } else if(p==1) {
      TF1 * g1 = new TF1("1ph","gaus",xp-14,xp+14);
      h->Fit(g1,"QR+");
      g1->GetParameters(&par[3]);
      delete g1;
     } else if(p==2) {
      TF1 * g2 = new TF1("2ph","gaus",xp-14,xp+14);
      h->Fit(g2,"QR+");
      g2->GetParameters(&par[6]);
      delete g2;
    } else {
      TF1 * g3 = new TF1("3ph","gaus",xp-14,xp+14);
      h->Fit(g3,"QR+");
      g3->GetParameters(&par[9]);
      delete g3;
    }
  }
  if(nfound==1) total = new TF1("mstotal","gaus(0)",xpeaks[0]-14,xpeaks[0]+14);
  else if(nfound==2) total = new TF1("mstotal","gaus(0)+gaus(3)",xpeaks[0]-14,xpeaks[1]+14);
  else if(nfound==3) total = new TF1("mstotal","gaus(0)+gaus(3)+gaus(6)",xpeaks[0]-14,xpeaks[2]+14);
  else total = new TF1("mstotal","gaus(0)+gaus(3)+gaus(6)+gaus(9)",xpeaks[0]-14,xpeaks[3]+14);
  total->SetParameters(par);
  h->Fit(total,"QR0+");

  float mean1 = total->GetParameter(1);
  //float error1 = total->GetParError(1);
  float sigma1 = total->GetParameter(2);

  float mean2 = total->GetParameter(4);
  //float error2 = total->GetParError(4);
  */


/*
int test(int channel) {
  char name[128];
  sprintf(name,"hcharge2_0;%d",channel);
  TH1D *h=(TH1D *) gDirectory->Get(name);
  //h->Draw("");
  double par[12];
  TSpectrum *s = new TSpectrum(4);
  int nfound = s->Search(h,10);
  //printf("Found %d candidate peaks to fit\n",nfound);
  float *xpeaks = s->GetPositionX();

  TF1 *g0 = new TF1("0ph","gaus",xpeaks[0]-3,xpeaks[0]+3);
  h->Fit(g0,"QR");
  g0->GetParameters(&par[0]);

  TF1 *g1 = new TF1("1ph","gaus",xpeaks[1]-3,xpeaks[1]+3);
  h->Fit(g1,"QR+");
  g1->GetParameters(&par[3]);
 

 
  //Float_t gainErr = sRt((error1/mean1)*(error1/mean1)+(error2/mean2)*(error2/mean2))
  delete s;
  return nfound;
}*/