/***********************************************
funtion to fit an arbitrary top threshlold
 **********************************************/
void topGraphs(TH1D *prob,Double_t x[],Double_t y[],Double_t y1[],const Int_t nPoints){

  //arrays that determine graph size
  //const Int_t nPoints=sizeof(x);

  //graphs

  //beginning of graph
  const Double_t xMin=347;
  //end of graph
  const Double_t xMax=352;
  //graph step size
  const Double_t xStep=(xMax-xMin)/nPoints;

  for (Int_t i=0;i<nPoints;i++){

    //x that we calculate cross section from
    Double_t myx=xMin+i*xStep;
    //lower limit of the integral
    Double_t myxMin=0.9*myx;
    //integral step size
    Double_t myxStep=(myx-myxMin)/100;

    //fill the simple graph array
    x[i]=myx;
    y[i]=topSection(myx);
    //cleanSec->SetPoint(i,myx,topSection(myx));

    //make histo to integrate
    TH1D *smear=new TH1D("smear","smear histogram",100,xMin,xMax);
    //fill it with product of spectrum*cross-section

    for (Int_t jLoop=1;jLoop<101;jLoop++){
      smear->SetBinContent(jLoop,prob->GetBinContent(jLoop)*
            topSection(myxMin+jLoop*myxStep));
    }
    //integrate
    y1[i]=smear->Integral();

    //delete loop histo
    delete smear;
  }



  return;
}


//a simple functional form for
//the top cross section


Double_t g1(Double_t x,Double_t s){
  Double_t val;

  if (x < 0)
    val=  exp(-pow( x-s ,2));
  else
    val= 1- exp(-pow(- x-s ,2));

  return val +0.1 +0.5*exp(-pow(x,2));

}

Double_t topSection(Double_t x){
  Double_t mean=349;
  Double_t s=0.8325546111576978;
  return g1(x-mean,s);
}




//flat top cross section
Double_t topSectionLinear(Double_t x){

  Double_t Wtop=3;
  Double_t Mtop=349.5;
  Double_t MaxSection=1.0;

  if (x<(Mtop-Wtop/2)){
    return 0;
  }
  else if ((x>(Mtop-Wtop/2)) && (x<(Mtop+Wtop/2)) ){
    Double_t f=MaxSection/Wtop*(x-(Mtop-Wtop/2));
    return f;
  }
  else if (x>=(Mtop+Wtop/2)){
    return MaxSection;
  }

}