// the following sourcecode is an excerpt of the actual program I use
// but one problem at a time...
// I run this from within root via .L fitterHappier.c and then calling
// the function with the arguments I actually need.

void swap(Float_t *pa, Float_t *pb)
{
	Float_t tmp=*pa;
	*pa=*pb;
	*pb=tmp;
}

void bubbleSort(Float_t *array, Int_t size)
{
	for( int obergrenze = size - 1; obergrenze>0; --obergrenze )
	{
		for( int pos = 0; pos < obergrenze; ++pos )
		{
			if( array[ pos ] > array[ pos + 1 ] )
			{ 
				swap( &array[ pos ], &array[ pos + 1 ] );
			}
		}
	}
}

void fitterHappier(char filename[1024], Int_t numPeaks=5, Float_t rangeSpace=5, Float_t pPos=0)
{
	char* title = "current fit";
	gROOT->Reset();
	ifstream infile;
	Float_t x,min,max;
	Int_t nlines	= 0;
	Int_t sum	= 0;

	infile.open(filename);
	infile >> min >> max;

	char nuTit[200];
	sprintf(nuTit,"%s (%s)", title, filename);
	TH1F *h1 = new TH1F("h1", nuTit, max+1, min, max);

	while (1) {
		infile >> x;
		if (!infile.good()) break;
		h1->SetBinContent(nlines, x);
		// need this later in the actual program - just found out today
		// that TH1F does not only provide the number of entries (which
		// simply returns the number of bins in this case) has this
		// already implemented...
		sum	+= x;
		nlines++;
	}

	if(nlines==0){
		cout << filename << " does not contain data!" << endl;
		break;
	}

	infile.close();
	h1->GetXaxis()->SetTitle("ADC-Channels");
	h1->GetYaxis()->SetTitle("Number of Events");

	TSpectrum *mySpec	= new TSpectrum(numPeaks);
	Int_t myNumPeaks	= mySpec->Search(h1,3,"new", 0.01);
	Float_t *myPeakPositions = mySpec->GetPositionX();

	// TSpectrum::Search() does somehow provide the peaks unordered
	// on the 5.10-Mac so I sort them first
	bubbleSort( myPeakPositions, myNumPeaks );

	char funcDef[300] = "gaus(0)";
	for(Int_t i = 1; i < myNumPeaks; i++){
		sprintf(funcDef,"%s + gaus(%i)",funcDef,i*3);
	}

	// coarse approximation to the spacing between the peaks
	// due to the nature of the signal they will be equidistant
	Float_t dist	= myPeakPositions[1] - myPeakPositions[0];
	TF1 *myFunc		= new TF1("myFunc", funcDef, 332, myPeakPositions[myNumPeaks-1]+dist/2-rangeSpace);

	for(Int_t i=0; i<myNumPeaks; ++i){
		myFunc->SetParameter(i*3, ( h1->GetBinContent(myPeakPositions[i]) )*3);
		myFunc->SetParameter(i*3+1,myPeakPositions[i]);
		myFunc->SetParameter(i*3+2,(dist - rangeSpace)/2);
		if (pPos!=0){
			myFunc->SetParLimits( (myNumPeaks-1)*3+1, myPeakPositions[i]-pPos, myPeakPositions[i]+pPos);
		}
	}
	

	h1->Fit("myFunc","RLQ");
	myFunc->SetRange(0,500);
	myFunc->SetLineStyle(2);
	myFunc->DrawCopy("SAME");
	h1->GetXaxis()->SetRangeUser(myPeakPositions[0]-dist,myPeakPositions[0]+dist*(myNumPeaks)*2);
}