// dichiaro la funzione, la definiro' poi
double myFunc(double *x, double *par);

// il mio "main"
void fitDoppio() {
  // istogramma
  TH1F* h=new TH1F("h","Double particle spectrum",100,0,20);

  // prima particella
  const int nPart1 = 1000;
  const float mass1 = 4.2;
  const float sigma1 = 0.6;

  // seconda particella
  const int nPart2 = 500;
  const float mass2 = 5.1;
  const float sigma2 = 0.3;

  // fondo
  const int nBkgn= 1700;

  // riempiamo isto con prima particella
  for (int i = 0; i < nPart1; ++i) {
    float rnd = gRandom->Gaus(mass1, sigma1);
    h->Fill(rnd);
  }

  // riempiamo isto con seconda particella
  for (int i = 0; i < nPart2; ++i) {
    float rnd = gRandom->Gaus(mass2, sigma2);
    h->Fill(rnd);
  }

  // mettiamo il fondo
  for (int i = 0; i<nBkgn; ++i) {
    float rnd=gRandom->Rndm()*20.;
    h->Fill(rnd);
  }

  // intanto guardiamo che faccia ha...
  h->DrawCopy();

  TF1 func("func",myFunc,0.,20.,7);
  func.SetParLimits(2,0,100);
  func.SetParameter(1,700);
  func.SetParameter(2,4);
  func.SetParameter(3,1);
  func.SetParLimits(5,0,100);
  func.SetParameter(4,500);
  func.SetParameter(5,16);
  func.SetParameter(6,0.9);
  func.SetParameter(7,15);

  h->Fit("func","r");
}

// definisco la mia funzione
double myFunc(double *x, double *par) {
  // mi sembra di vedere due picchi e fondo uniforme, facci funzione
  // gaus1+gaus2+flat

  return 
    par[0]*exp(-((x[0]-par[1])/(2*par[2]))*((x[0]-par[1])/(2*par[2]))) + // gaus1
    par[3]*exp(-((x[0]-par[4])/(2*par[5]))*((x[0]-par[4])/(2*par[5]))) + // gaus2
    par[6]; // flat
}