//*-- AUTHOR : J.Kempter
//*-- Modified :

//_HADES_CLASS_DESCRIPTION 
///////////////////////////////////////////////////////////////////////////////
//
// HMdcOffset
//
// Defines the offset  parameter for MDC calibration. Uses HMdcCalParRaw container
// as  input/output for calibration data
//
///////////////////////////////////////////////////////////////////////////////
using namespace std;
#include <math.h>
#include <stdlib.h>
#include <fstream> 
#include <iomanip> 
#include <iostream> 

#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TStopwatch.h>
#include <TDirectory.h>
#include <TNtuple.h>
#include <TMath.h>

#include "heventheader.h"
#include "hmdcoffset.h"
#include "hmdcdef.h"
#include "hmdctrackddef.h"
#include "hmdcsizescells.h"
#include "hmdcraw.h"
#include "hmdcclus.h"
#include "hdebug.h"
#include "hades.h"
#include "hiterator.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hmdcdetector.h"
#include "hdetector.h"
#include "hevent.h"
#include "hcategory.h"
#include "hlocation.h"
#include "hmdccal1.h"
#include "hmdclookupgeom.h"
#include "hmdclookupraw.h"
#include "hmdccalpar.h"
#include "hmdccalparraw.h"
#include "hmdclookupraw.h"
#include "hstarthit.h"
#include "hstartdef.h"
#include "hmdctimecut.h"

ClassImp(HMdcOffset)

const Int_t HMdcOffset::nbin        = 2048;
const Int_t HMdcOffset::nbinp1      = HMdcOffset::nbin+1;
const Int_t HMdcOffset::nbinm1      = HMdcOffset::nbin-1;
const Int_t HMdcOffset::nSubEvents  = (6*192)*3;


HMdcOffset::HMdcOffset(void)
{
    setDefault();
    initVariables();
    //initMemory();
}

HMdcOffset::HMdcOffset(Text_t* name,Text_t* title) : HReconstructor(name,title)
{
    setDefault();
    initVariables();
    //initMemory();
}
HMdcOffset::~HMdcOffset(void)
{
    // Destructor.
    if (iter)             delete iter;
    if (iter_start)       delete iter_start;
    if (iter_clus)        delete iter_clus;
    if (fNameAsciiOffset) delete fNameAsciiOffset;
    if (fNameRootOffset)  delete fNameRootOffset;
    if (hreverse)         delete[] hreverse;
    iter = NULL;
    iter_start= NULL;
    iter_clus = NULL;

}
void HMdcOffset::initVariables()
{
    // pointer to iterators and categories
    rawCat       =0;
    hitStartCat  =0;
    clusCat      =0;
    iter         =0;
    iter_start   =0;
    iter_clus    =0;
    locraw.set(4,0,0,0,0);

    // switch options
    isPulserFile =kFALSE;
    noStart      =kFALSE;
    useTimeCuts  =kTRUE;
    useClusters  =kTRUE;
    useWireOffset=kTRUE;
    useTof       =kFALSE;
    debug        =kFALSE;
    // pointers to fileds and hists
    hreverse     = (MyField*) new MyField;
    hint         =0;
    hinv         =0;

    // container pointers
    calparraw    =0;
    lookupgeom   =0;
    timecut      =0;
    sizescells   =0;

    // constants
    signalspeed=0.004; // [ns/mm]

    // varables in finalize
    validRange   =1000.;
    meanhOffset  =0;
    myoffset=0;
    myerror=0;

    // counters
    eventcounter =0;
    skipcounter  =0;
    nSkipEvents  =0;

    // output file names
    fNameAsciiOffset =0;
    fNameRootOffset  =0;
    ferrorlog=0;
    offsetTuple      =0;
    offsetPulserTuple=0;

    // fit variables
    yequalzero       =0;
    crosspointX      =0;
    fitpar0          =0;
    fitpar0error     =0;
    fitpar1          =0;
    fitpar1error     =0;

    fitparNoise0     =0;
    fitparNoise0error=0;
    fitparNoise1     =0;
    fitparNoise1error=0;
    totalsigma       =0;

    fitGaussMean     =0;
    fitGaussSigma    =0;

    
}
void HMdcOffset::setDefault()
{
    // Sets the the default values for the fits.
    setNoise(100, 50);
    setThreshold(0.15, 0.5);
    setRangeGauss(50);
}

void HMdcOffset::setOutputAscii(Char_t *c)
{
    // Sets ascii output of HMdcOffset for debugging output.
    if (fNameAsciiOffset) delete fNameAsciiOffset;
    fNameAsciiOffset = new Char_t[strlen(c)+1];
    strcpy(fNameAsciiOffset, c);
}

void HMdcOffset::setOutputRoot(Char_t *c)
{
    // Sets rootfile output of HMdcOffset where all created histograms were written.
    //
    if (fNameRootOffset) delete fNameRootOffset;
    fNameRootOffset = new Char_t[strlen(c)+1];
    strcpy(fNameRootOffset, c);
}
void HMdcOffset::setUseTof(TString filename)
{
    // Retrieves TF1's for min tof substraction from root file
    // and switches the usesage of tof substraction kTRUE.
    if (filename.CompareTo("")==0)
    {   // check for empty string
	useTof=kFALSE;
	Warning("setUseTof()","Empty String received for inputfile name of tof parameters!");
        exit(1);
    }
    else
    {   // if string is not empty
	TFile* inp=new TFile(filename,"READ");
	if(inp==0)
	{   // if file does not exist
	    Warning("setUseTof()","Input file not existent!");
	    exit(1);
	}
	else
	{   // if file exists
	    cout<<"HMdcOffset::setUseTof(): Reading from "<<filename.Data()<<endl;
	    Char_t fitname[300];
	    for(Int_t mdc=0;mdc<4;mdc++){
		for(Int_t lay=0;lay<6;lay++){
		    sprintf(fitname,"fit[%i][%i]",mdc,lay);
                    toffunc[mdc][lay]=(TF1*)inp->Get(fitname);
		    if(toffunc[mdc][lay]==0)
		    {   // if fit object has not been found in file
			Warning("setUseTof()","Error in reading input file!");
                        exit(1);
		    }
		}
	    }
	    useTof=kTRUE;
	}
    }
}
void HMdcOffset::printStatus()
{
    printf("**************** HMdcOffset ***********************\n");
    printf("* ACTUAL SETTINGS                                 *\n");
    if(useTimeCuts)
    {
	printf("*    TIME CUTS are used                           *\n");
    }
    if(!useTimeCuts)
    {
	printf("*    NO TIME CUTS are used                        *\n");
    }
    if(useClusters)
    {
	printf("*    TIMES FILLED FROM CLUSTERS                   *\n");
    }
    if(!useClusters)
    {
	printf("*    TIMES FILLED FROM RAW                        *\n");
    }
    if(useWireOffset&&useClusters)
    {
	printf("*    TIME OF SIGNAL ON WIRE SUBTRACTED            *\n");
    }
    if(!useWireOffset)
    {
	printf("*    TIME OF SIGNAL ON WIRE NOT SUBTRACTED        *\n");
    }
    if(useTof)
    {
	printf("*    MINIMUM TOF WILL BE SUBSTRACTED              *\n");
    }
    if(!useTof)
    {
	printf("*    MINIMUM TOF WILL NOT BE SUBSTRACTED          *\n");
    }
    if(isPulserFile)
    {
	printf("*    Mode set to PULSER FILE                      *\n");
    }
    if(!isPulserFile)
    {
	printf("*    Mode set to REAL DATA FILE                   *\n");
    }

    if(!noStart)
    {
	printf("*    START TIMES are used                         *\n");
    }
    if(noStart)
    {
	printf("*    NO START TIMES are used                      *\n");
    }
    printf("*    SIGNALSPEED = %5.3f [ns/mm]                 *\n",signalspeed);
    printf("*    VALID OFFSET RANGE = %5.1f                  *\n",validRange);
    printf("*    NOISE OFFSET = %4i                          *\n",offsetfitNoise);
    printf("*    NOISE WIDTH  = %4i                          *\n",widthfitNoise);
    printf("*    MIN Threshold   = %1.2f MAX Threshold = %1.2f  *\n",minfitthreshold,maxfitthreshold);
    printf("*    Range GAUSS FIT = %4i                       *\n",rangeGauss);
    printf("**************** HMdcOffset ***********************\n");
}
Bool_t HMdcOffset::init(void)
{
    //  Inits HMdcOffset and the needed HMdcCalParRaw, if this container does not exists
    if(!noStart)
    {
	hitStartCat=gHades->getCurrentEvent()->getCategory(catStartHit);       //????
	if (!hitStartCat)
	{
          Error("init()","NO START HIT CATEGORY IN INPUT!");
          return kFALSE;
	}
	else iter_start=(HIterator *)hitStartCat->MakeIterator("native");
    }
    if(useClusters)
    {
	clusCat=gHades->getCurrentEvent()->getCategory(catMdcClus);       //????
	if (!clusCat)
	{
          Error("init()","NO CLUS CATEGORY IN INPUT!");
          return kFALSE;
	}
	else iter_clus=(HIterator *)clusCat->MakeIterator("native");

	lookupraw =(HMdcLookupRaw*) gHades->getRuntimeDb()->getContainer("MdcLookupRaw");
	if(useWireOffset)
	{
	    sizescells=HMdcSizesCells::getExObject();
	    if(!sizescells)
	    {
		Error("init()","NO HMDCSIZESCELLS CONTAINER IN INPUT!");
		return kFALSE;
	    }
	}
    }

    calparraw =(HMdcCalParRaw*) gHades->getRuntimeDb()->getContainer("MdcCalParRaw");
    timecut   =(HMdcTimeCut*)   gHades->getRuntimeDb()->getContainer("MdcTimeCut");
    lookupgeom=(HMdcLookupGeom*)gHades->getRuntimeDb()->getContainer("MdcLookupGeom");

    rawCat=gHades->getCurrentEvent()->getCategory(catMdcRaw);
    if (!rawCat)
    {
	Error("init()","NO MDC RAW CATEGORY IN INPUT!");
	return kFALSE;
    }
    else {
	// creates an iterator which loops over all fired cells
	iter=(HIterator *)rawCat->MakeIterator();
    }
    fActive=kTRUE;
    eventcounter=0;
    skipcounter=0;
    printStatus();
    return kTRUE;
}

Bool_t HMdcOffset::reinit(void)
{
    //
    cout<<"skipcount "<<skipcounter<<endl;
    skipcounter=0;

    fActive=kTRUE;
    return kTRUE;
}



void HMdcOffset::createHist(TFile* file,Int_t s, Int_t m, Int_t l, Int_t c)
{
    // Histograms for inverted Time1 and integrated Time1 per Tdc-Channel are created
    // in a subdirectory structure.
    static Char_t title[50], tmp[30];

    sprintf(title,"%s%i%s%i%s%i%s%i","Sector ",s,"  Module ",m,"  Mbo ", l,"  Tdc ", c);
    file->TDirectory::Cd("hinv");
    sprintf(tmp, "%s%i%s%i%s%02i%s%02i%s", "hinv[",s,"][",m,"][",l,"][",c,"]");
    hinv = new MyHist(tmp, title, nbin,-nbin+0.5,0.5); //inverted hist.
    file->TDirectory::Cd("../hint");
    sprintf(tmp, "%s%i%s%i%s%02i%s%02i%s", "hint[",s,"][",m,"][",l,"][",c,"]");
    hint = new MyHist(tmp, title, nbin,-nbin+0.5,0.5); // integrated hist.
    file->TDirectory::Cd("..");
}
void HMdcOffset::createHist_2D(Int_t s, Int_t m)
{
    // Histograms for inverted Time1 and integrated Time1 per Tdc-Channel are created
    // in a subdirectory structure.
    static Char_t title[50], tmp[30];
    for(Int_t mb=0;mb<16;mb++)
    {
	sprintf(title,"%s%i%s%i%s%i","Sector ",s,"  Module ",m,"  Mbo ", mb);
	sprintf(tmp, "%s%i%s%i%s%02i%s", "htime1_mbo[",s,"][",m,"][",mb,"]");
	htime1_mbo[mb] = new TH2F(tmp, title,700,-200,500,96,1,97); //hist time vers tdc.
    }
    for(Int_t lay=0;lay<6;lay++)
    {
	sprintf(title,"%s%i%s%i%s%i","Sector ",s,"  Module ",m,"  Lay ", lay);
	sprintf(tmp, "%s%i%s%i%s%02i%s", "htime1_lay[",s,"][",m,"][",lay,"]");
	htime1_lay[lay] = new TH2F(tmp, title,700,-200,500,210,1,211); //hist time vers tdc.
    }
    for(Int_t lay=0;lay<6;lay++)
    {
	sprintf(title,"%s%i%s%i%s%i","Sector ",s,"  Module ",m,"  Lay ", lay);
	sprintf(tmp, "%s%i%s%i%s%02i%s", "htime1_lay_inv_norm[",s,"][",m,"][",lay,"]");
	htime1_lay_inv_norm[lay] = new TH2F(tmp, title,700,-200,500,210,1,211); //hist time vers tdc.
    }
    for(Int_t lay=0;lay<6;lay++)
    {
	sprintf(title,"%s%i%s%i%s%i","Sector ",s,"  Module ",m,"  Lay ", lay);
	sprintf(tmp, "%s%i%s%i%s%02i%s", "htime1_lay_int[",s,"][",m,"][",lay,"]");
	htime1_lay_int[lay] = new TH2F(tmp, title,700,-200,500,210,1,211); //hist time vers tdc.
    }

    for(Int_t lay=0;lay<6;lay++)
    {
	sprintf(title,"%s%i%s%i%s%i","Sector ",s,"  Module ",m,"  Lay ", lay);
	sprintf(tmp, "%s%i%s%i%s%02i%s", "htime1_lay_int_norm[",s,"][",m,"][",lay,"]");
	htime1_lay_int_norm[lay] = new TH2F(tmp, title,700,-200,500,210,1,211); //hist time vers tdc.
    }

}
void HMdcOffset::deleteHist()
{
    // Created histograms are deleted
    delete hinv;
    delete hint;
}
void HMdcOffset::deleteHist_2D()
{
    // Created histograms are deleted
    for(Int_t mb=0;mb<16;mb++)
    {
      delete htime1_mbo[mb];
    }
    for(Int_t lay=0;lay<6;lay++)
    {
	delete htime1_lay[lay];
	delete htime1_lay_inv_norm[lay];
	delete htime1_lay_int[lay];
	delete htime1_lay_int_norm[lay];
    }

}

void HMdcOffset::fillHist(Int_t s, Int_t m, Int_t l, Int_t c)
{
    // Histograms for inverted Time1 and integrated Time1 per Tdc-Channel are filled

    hint->SetBinContent(1,(Double_t)(*hreverse)[s][m][l][c][0]);
    hinv->SetBinContent(1,(Double_t)(*hreverse)[s][m][l][c][0]);

    for(Int_t k=2; k<nbinp1; k++)
    {
	hinv->SetBinContent(k,(Double_t)(*hreverse)[s][m][l][c][k-1]);
	hint->SetBinContent(k,(hint->GetBinContent(k-1)+hinv->GetBinContent(k)));
    }
    integral[s][m][l][c]=(Int_t)hint->GetBinContent(nbin-2);
}
void HMdcOffset::fillHist_2D(Int_t s, Int_t m, Int_t l, Int_t c)
{
    // Histograms for Time1 vers Tdc-Channel and cell are filled
    HMdcLookupChan& chan=(*lookupgeom)[s][m][l][c];
    Int_t layer= chan.getNLayer();
    Int_t cell= chan.getNCell();
    Int_t myoffbin=(Int_t)(2048-offsets[s][m][l][c]); // convert offset x to offset bin
    Float_t temp;
    for(Int_t k=1; k<nbinp1; k++) // loop over all bins of hinv  (x from -2048 to 0)
    {
	if(k-myoffbin<500&&k-myoffbin>-200) // check if bin in range
	{
	    Int_t mybin=200+k-myoffbin;

	    temp=htime1_mbo[l]->GetCellContent(mybin,c+1);
	    htime1_mbo[l]    ->SetCellContent(mybin,c+1   ,hinv->GetBinContent(k)+temp);

	    temp=htime1_lay[layer]->GetCellContent(mybin,cell+1);
	    htime1_lay[layer]->SetCellContent(mybin,cell+1,hinv->GetBinContent(k)+temp);

	    temp=htime1_lay_int[layer]->GetCellContent(mybin,cell+1);
	    htime1_lay_int[layer]->SetCellContent(mybin,cell+1,hint->GetBinContent(k)+temp);

	    temp=htime1_lay_inv_norm[layer]->GetCellContent(mybin,cell+1);
	    htime1_lay_inv_norm[layer]->SetCellContent(mybin,cell+1,hinv->GetBinContent(k)+temp);

	    temp=htime1_lay_int_norm[layer]->GetCellContent(mybin,cell+1);
	    htime1_lay_int_norm[layer]->SetCellContent(mybin,cell+1,hint->GetBinContent(k)+temp);

	}
    }

    Float_t max_inv=hinv->GetMaximum();
    Float_t max_int=hint->GetMaximum();

    for(Int_t k=1; k<nbinp1; k++) // loop over all bins of hinv  (x from -2048 to 0)
    {
	if(max_inv!=0)
	{
	    temp=htime1_lay_inv_norm[layer]->GetCellContent(k,cell+1);
	    htime1_lay_inv_norm[layer]->SetCellContent(k,cell+1,temp/max_inv);
	}
	if(max_int!=0)
	{
	    temp=htime1_lay_int_norm[layer]->GetCellContent(k,cell+1);
	    htime1_lay_int_norm[layer]->SetCellContent(k,cell+1,temp/max_int);
	}
    }

}

Int_t HMdcOffset::fitHist(Int_t s, Int_t m, Int_t l, Int_t c)
{
    // The offset is calculated from two linear fits. Some pictures of the fitted histograms
    // can be found on Begin_Html <a href="http://www-hades.gsi.de/docs/mdc/mdcana/#calibration">MDC calibration page</a>End_Html .
    // The first linear fit (red) is done to find the rising edge of the integrated spectra.
    // The second linear fit (blue) is done to substract the noise.
    // The ranges for both fits can be set through the functions setNoise()(x-range) for the
    // second fit (default values: 100,50) and setThreshold()(y-range) for the first fit
    // (default values :0.15,0.50). The y-range is calculated by the percentage of the maximum
    // height of the spectra.
    // The offset ist calculated from the intersection point of both linear fits.
    // The sigma of the offset is calculated from the sigmas of the two linear fits.
    //
    // If the analyzed file is a pulser file for external calibration a gaussian fit is done.
    // To set this fit the function setPulserFile() has to be used in the macro. The range of
    // the fit is set through the function setRangeGauss()(default value: 50) around the maximum.
    // The mean and sigma of the fit is written to the ascii debugging output.
    // The hist is checked for multiple peaks.The peaks which have been found are written to an array.

    fitGaussMean  = 0;
    fitGaussSigma = 0;
    fitpar0     =0;
    fitpar0error=0;
    fitpar1     =0;
    fitpar1error=0;
    fitparNoise0     =0;
    fitparNoise1     =0;
    fitparNoise0error=0;
    fitparNoise1error=0;
    totalsigma=0;
    yequalzero=0;
    crosspointX=0;
    //-----------------------------------------------------------
    // if hist is empty do nothing
    if(hint->Integral()==0)
    {   
        if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: Integral==0 "<<endl;
	yequalzero=1;
	crosspointX=1;
	return 1;

    }
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // Check for overflows. In case of an overflow bins
    // might be negative
    if(hint->GetMinimum()<0)
    {
        cout<<" "<<endl;
	cout<<"WARNING: IN HINT "<<s<<" "<<m<<" "<<l<<" "<<c<<" NEGATIV VALUES DETECTED!"<<endl;
        cout<<" "<<endl;
	fprintf(ferrorlog,"%i %i %i %2i %s\n",s,m,l,c,"fit: in hint negative values detected!");
    }
    //-----------------------------------------------------------


    //-----------------------------------------------------------
    // finding upper fitting range for signal fit
    Int_t   xmaxfitbin=nbinm1;

    Float_t yminfit=minfitthreshold*hint->GetBinContent(nbinm1-1);     // setting y thresholds for fitting
    Float_t ymaxfit=maxfitthreshold*hint->GetBinContent(nbinm1-1);     // relative to maximum y-value
    while (hint->GetBinContent(--xmaxfitbin) > ymaxfit && xmaxfitbin); // find to ymaxfit corresponding bin
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // if no match for the upper fit range has been found in
    // the desired Y range
    if (!xmaxfitbin)
    {
	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: xmaxfitbin=0, interal "<<hint->Integral()<<endl;
	yequalzero=1;
	crosspointX=1;
	return 1;
    }
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // finding lower fitting range for signal fit
    Int_t xminfitbin=xmaxfitbin;
    while (hint->GetBinContent(--xminfitbin) > yminfit && xminfitbin); // find to yminfit corresponding bin
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // if no match for the lower fit range has been found in
    // the desired Y range
    if (!xminfitbin)
    {
	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: xminfitbin=0, interal "<<hint->Integral()<<endl;
	yequalzero=2;
	crosspointX=2;
	return 2;
    }
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // if fit range is just 1 bin the fit will give an error message.
    // this is just to suppress the error message.
    if(xmaxfitbin-xminfitbin==1)
    {   // fit range is just 1 bin!
	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: fit range 1 bin, interal "<<hint->Integral()<<endl;
	yequalzero=2;
	crosspointX=2;
	return 2;

    }
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // create pol1 fit function in desired range and fit the hist.
    TF1 *f = new TF1("fitsignal","pol1",hint->GetBinCenter(xminfitbin),hint->GetBinCenter(xmaxfitbin)); // call fit function poynom 1.order in x-range
    hint->Fit("fitsignal","RNQ");
    fitpar0     =0;
    fitpar0error=0;
    fitpar1     =0;
    fitpar1error=0;
    fitpar0     =f->GetParameter(0);  // get fitparameters
    fitpar0error=f->GetParError (0);
    fitpar1     =f->GetParameter(1);  //(f1->GetParameter(1)<0.000001)?-1:(f1->GetParameter(1));
    fitpar1error=f->GetParError (1);  // y=fitpar[1]* x +  fitpar[0]
    //-----------------------------------------------------------
    

    //-----------------------------------------------------------
    // if result of fit is nan or inf.
    if(TMath::Finite(f->GetParameter(0))==0 ||
       TMath::Finite(f->GetParameter(1))==0
      )
    {   // bullshit !!
	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: fit result strange, integral "
	        <<hint->Integral()
		<<" xminbin "<<xminfitbin
		<<" xmaxbin "<<xmaxfitbin
		<<" fitpar0 "<<f->GetParameter(0)
		<<" fitpar1 "<<f->GetParameter(1)
		<<endl;
	yequalzero=3;
	crosspointX=3;
	delete f;
	return 3;
    }
    //-----------------------------------------------------------


    //-----------------------------------------------------------
    // remove fit function
    delete f;
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // if fitpar1 is 0 we can not do anything usefull
    if (fitpar1==0)
    {
	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: fitpar1==0 , integral "
	        <<hint->Integral()
		<<" xminbin "<<xminfitbin
		<<" xmaxbin "<<xmaxfitbin
		<<" fitpar0 "<<fitpar0
		<<" fitpar1 "<<fitpar1
		<<endl;
	yequalzero=3;
	crosspointX=3;
	return 3;
    }
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // calculate cross point of fit with x-axis
    // and check for good range
    yequalzero=0;
    yequalzero =-fitpar0/fitpar1;

    if (-yequalzero < (offsetfitNoise+widthfitNoise) || -yequalzero > nbin)
    {
	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: yequalzero out of range, integral "
	        <<hint->Integral()
                <<" -yequalzero "<<-yequalzero
		<<" xminbin "<<xminfitbin
		<<" xmaxbin "<<xmaxfitbin
		<<" fitpar0 "<<fitpar0
		<<" fitpar1 "<<fitpar1
		<<endl;
	crosspointX=4;
	return 4;
    }
    //-----------------------------------------------------------

    //-----------------------------------------------------------
    // check weather the range for the noise fit contains any
    // counts. This is just to suppress error messages from the
    // fit. If the is no signal yequalzero will be taken as offset.
    if(hint->Integral(hint->FindBin(yequalzero-(offsetfitNoise+widthfitNoise)),
		      hint->FindBin(yequalzero-offsetfitNoise))<10
      )
    {  // almost no counts in range of fit!

	if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: no counts in noise region , integral "
	        <<hint->Integral()
                <<" -yequalzero "<<-yequalzero
		<<" minbinNoise "<<hint->FindBin(yequalzero-(offsetfitNoise+widthfitNoise))
		<<" maxbinNoise "<<hint->FindBin(yequalzero-offsetfitNoise)
		<<endl;
	crosspointX=-yequalzero;
	fitparNoise0     =0;
	fitparNoise1     =0;
	fitparNoise0error=0;
	fitparNoise1error=0;
    }
    else
    {
	//-----------------------------------------------------------
	// create pol1 fit function for the noise fit in
        // the desired range and fit the hist
	f = new TF1("fitnoise","pol1",
		    yequalzero-(offsetfitNoise+widthfitNoise),
		    yequalzero-offsetfitNoise);
	hint->Fit("fitnoise","RNQ");
	fitparNoise0     =0;
	fitparNoise1     =0;
	fitparNoise0error=0;
	fitparNoise1error=0;
	fitparNoise0     =f->GetParameter(0); // get fitparameters
	fitparNoise1     =f->GetParameter(1);
	fitparNoise0error=f->GetParError (0); // get fitparameters
	fitparNoise1error=f->GetParError (1);
	//-----------------------------------------------------------

	//-----------------------------------------------------------
	// if result of fit is nan or inf.
        // if so yequalzero is taken as offset
	if(TMath::Finite(f->GetParameter(0))==0 ||
	   TMath::Finite(f->GetParameter(1))==0
	  )
	{   // bullshit !!
	    if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: noise fit result strange , integral "
	             <<hint->Integral()
		     <<" -yequalzero "<<-yequalzero
		     <<" minbinNoise "<<hint->FindBin(yequalzero-(offsetfitNoise+widthfitNoise))
		     <<" maxbinNoise "<<hint->FindBin(yequalzero-offsetfitNoise)
		     <<" fitparNoise0 "<<f->GetParameter(0)
		     <<" fitparNoise1 "<<f->GetParameter(1)
		     <<endl;
	    crosspointX=-yequalzero;
	    fitparNoise0     =0;
	    fitparNoise1     =0;
	    fitparNoise0error=0;
	    fitparNoise1error=0;
       }

	//-----------------------------------------------------------
	// if fitpar1==fitparNoise1 we can not calculate
        // any cross point. if so yequalzero is taken as offset
	if (fitpar1==fitparNoise1)
	{
	    if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: no intersection for fit functions , integral "
	             <<hint->Integral()
		     <<" -yequalzero "<<-yequalzero
		     <<" minbinNoise "<<hint->FindBin(yequalzero-(offsetfitNoise+widthfitNoise))
		     <<" maxbinNoise "<<hint->FindBin(yequalzero-offsetfitNoise)
                     <<" fitpar0     "<<fitpar0
                     <<" fitpar1     "<<fitpar1
                     <<" fitparNoise0 "<<f->GetParameter(0)
		     <<" fitparNoise1 "<<f->GetParameter(1)
		     <<endl;
	    crosspointX=-yequalzero;
	    fitparNoise0     =0;
	    fitparNoise1     =0;
	    fitparNoise0error=0;
	    fitparNoise1error=0;
	}

	
	//-----------------------------------------------------------
	// finally calculate cross point between noise and signal
	// fit to define the offset
	crosspointX=-(fitparNoise0-fitpar0)/(fitpar1-fitparNoise1);
	//-----------------------------------------------------------

	//-----------------------------------------------------------
	// check if value is ok and in meaningfull range
	if (TMath::Finite(crosspointX)==0 || crosspointX<0 || crosspointX>nbinm1)
	{
	    if(debug)cout<<s<<" "<<m<<" "<<l<<" "<<c<<" skip: crospointX out of range or strange , integral "
	             <<hint->Integral()
		     <<" -yequalzero "<<-yequalzero
                     <<" crosspaintX "<<crosspointX
             	     <<" xminbin "<<xminfitbin
		     <<" xmaxbin "<<xmaxfitbin
		     <<" minbinNoise "<<hint->FindBin(yequalzero-(offsetfitNoise+widthfitNoise))
		     <<" maxbinNoise "<<hint->FindBin(yequalzero-offsetfitNoise)
                     <<" fitpar0     "<<fitpar0
                     <<" fitpar1     "<<fitpar1
                     <<" fitparNoise0 "<<f->GetParameter(0)
		     <<" fitparNoise1 "<<f->GetParameter(1)
		     <<endl;
	    crosspointX=-yequalzero;
	    fitparNoise0     =0;
	    fitparNoise1     =0;
	    fitparNoise0error=0;
	    fitparNoise1error=0;
	}
 	//-----------------------------------------------------------

	//-----------------------------------------------------------
	// remove fit function
	delete f;
	//-----------------------------------------------------------
    }

    //-----------------------------------------------------------
    // calculating total sigma of the found offset
    totalsigma= (Float_t)sqrt((pow(fitparNoise0error,2)/pow((fitparNoise1-fitpar1),2))
			      +(pow(fitpar0error,2)/pow((fitparNoise1-fitpar1),2))
			      +(pow((fitparNoise0-fitpar0),2)*pow(fitpar1error,2)/pow((fitparNoise1-fitpar1),4))
			      +(pow((fitparNoise0-fitpar0),2)*pow(fitparNoise1error,2)/pow((fitparNoise1-fitpar1),4))
			     );
    //-----------------------------------------------------------

    // check if fitpars are meaningful
    if(fitpar1<100000&&fitpar1>-100000)
    {
	fitslope1[s][m][l][c]=fitpar1;
    }
    else
    {
	fitslope1[s][m][l][c]=0;
    }
    if(fitparNoise1<100000&&fitpar1>-100000)
    {
	fitslope2[s][m][l][c]=fitparNoise1;
    }
    else
    {
	fitslope2[s][m][l][c]=0;
    }

    if (!isPulserFile) return 0;


    //f = new TF1("f3","gaus", -crosspointX-rangeGauss, -crosspointX+rangeGauss);
    // If there are more than 1 peak only the most prominant is fitted
    Float_t localmax=(2048-(hinv->GetMaximumBin()));
    f = new TF1("f3","gaus",-localmax-rangeGauss, -crosspointX+rangeGauss);

    hinv->Fit("f3","RNQ");
    fitGaussMean =f->GetParameter(1);  //mean
    fitGaussSigma=f->GetParameter(2);  //sigma

    findMultiplePeaks(s,m,l,c);

    delete f;

    // FIXME: Error handling for gauss fit ???

    return 0;
}
void HMdcOffset::findMultiplePeaks(Int_t s,Int_t m,Int_t l,Int_t c)
{
    for(Int_t j=0;j<2047;j++)
    {
	htime1temp->SetBinContent(j,hinv->GetBinContent(j));// copy hinv to temp hist
    }
    htime1temp->Smooth(10);// make it smooth
    Char_t address[50];
    Float_t binlowold[5];
    Float_t binhighold[5];
    Float_t binlow[5];
    Float_t binhigh[5];
    Int_t binmax[5];
    Float_t offout[5];

    Float_t max=0;
    Float_t min=0;
    Float_t threshold=0.05;

    for(Int_t j=0;j<5;j++)// reset helper arrays
    {
	binlowold[j]=0;
	binhighold[j]=0;
	binlow[j]=0;
	binhigh[j]=0;
	binmax[j]=0;
	offout[j]=0;
    }
    max=min=0;
    max=htime1temp->GetMaximum(); // get max height
    min=threshold*max;
    Int_t count=0;
    for(Int_t bin=0;bin<2047;bin++)
    {
	if(count<5)
	{
	    binlow [count]=htime1temp->GetBinContent(bin);
	    binhigh[count]=htime1temp->GetBinContent(bin+1);

	    if(bin>40 && binlow[count]>20  // skip the first 40 tdc bins and take into account only if > 20 counts
	       && fabs(binhigh[count]-binlow[count])>1
	       && (binhighold[count]-binlowold[count]>0)
	       && (binhigh[count]-binlow[count]<0) ) // if zero was crossed
	    {
		binmax[count]=bin;
		count++;
	    }
	    if(count<5)
	    {
		binlowold [count]=binlow [count];
		binhighold[count]=binhigh[count];
	    }
            else
	    {
		cout<<"MORE THAN 5 PEAKS!"<<endl;
	    }
	}
    }
    if(count>1)
    {
        cout<<" "<<endl;
	cout<<"sector "<<s<<" module "<<m<<" mbo "<<l<<" tdc "<<c<<endl;
	for(Int_t i=0;i<count;i++)
	{
	    if(i<5)
	    {
		cout<<"peak found: "<<i+1<<" bin="<<binmax[i]<<" ,x="<<(2048-binmax[i])<<" max="<<binlow[i]<<endl;
                sprintf(address,"%i %i %i %2i",s,m,l,c);
		fprintf(ferrorlog,"%s %i %s %4i %s %4i %s %8.0f\n",address,i+1,"multiple peaks: bin=",binmax[i]," ,x=",(2048-binmax[i])," max=",binlow[i]);
	    }
	    else
	    {
		cout<<"MORE THAN 5 PEAKS!"<<endl;
	    }
	}
	cout<<" "<<endl;
    }

    for(Int_t j=0;j<5;j++)// coppy to global array
    {
	if(binmax[j]>0)
	{
	    offsetpulser[s][m][l][c][j]=(2048-binmax[j]);
	}
	else
	{
	    offsetpulser[s][m][l][c][j]=0;
	}
    }
}
void HMdcOffset::writeAscii(ofstream &fout, Int_t s, Int_t m, Int_t l, Int_t c)
{
    // The adresses of the channel(sector,module,mbo,tdc),offset, two fitparameters for the
    // two linear fits, the sigma of the offset ,offset and the integral(number of entries for this channel)
    // are written to the ascii debugging output.In the case of an pulser file (external calibration)
    // the mean and the sigma of the gaussian fit are also written to the output.
    if (isPulserFile)
        fout
            << setw(3)  << s                 << " "
            << setw(4)  << m                 << " "
            << setw(4)  << l                 << " "
            << setw(5)  << c                 << " "
	    << setw(10) << -yequalzero       << " "
            << setw(10) << yequalzero        << " "
	    << setw(10) << crosspointX       << " "
            << setw(15) << -fitGaussMean     << " "
            << setw(15) << fitGaussSigma     << " "
            << setw(15) << fitpar0           << " "
            << setw(15) << fitpar1           << " "
            << setw(15) << fitparNoise0      << " "
            << setw(15) << fitparNoise1      << " "
            << setw(15) << totalsigma        << " "
            << setw(10) << integral[s][m][l][c]

            << endl;
    else
        fout
            << setw(3)  << s                  << " "
            << setw(4)  << m                  << " "
            << setw(4)  << l                  << " "
            << setw(5)  << c                  << " "
	    << setw(10) << -yequalzero        << " "
	    << setw(10) << crosspointX        << " "
            << setw(15) << fitpar0            << " "
            << setw(15) << fitpar1            << " "
            << setw(15) << fitparNoise0       << " "
            << setw(15) << fitparNoise1       << " "
            << setw(15) << totalsigma         << " "
            << setw(10) << integral[s][m][l][c]
            << endl;
}

void HMdcOffset::writeHist(TFile* file)
{
    // All created histograms are written to a rootfile.The file is structured for sector,module,mbo.
    
    file->TDirectory::Cd("hinv");
    hinv->Write();
    file->TDirectory::Cd("../hint");
    hint->Write();
    file->TDirectory::Cd("..");
}
void HMdcOffset::writeHist_2D()
{
    // All created histograms are written to a rootfile.
    for(Int_t mb=0;mb<16;mb++)
    {
	htime1_mbo[mb]->SetEntries(1);
        htime1_mbo[mb]->SetOption("colz");
        htime1_mbo[mb]->Write();

    }
    for(Int_t lay=0;lay<6;lay++)
    {
	htime1_lay[lay]->SetEntries(1);
        htime1_lay[lay]->SetOption("colz");
        htime1_lay[lay]->Write();

	htime1_lay_inv_norm[lay]->SetEntries(1);
        htime1_lay_inv_norm[lay]->SetOption("colz");
        htime1_lay_inv_norm[lay]->Write();

	htime1_lay_int[lay]->SetEntries(1);
        htime1_lay_int[lay]->SetOption("colz");
        htime1_lay_int[lay]->Write();

	htime1_lay_int_norm[lay]->SetEntries(1);
        htime1_lay_int_norm[lay]->SetOption("colz");
        htime1_lay_int_norm[lay]->Write();
    }
}

ofstream *HMdcOffset::openAsciiFile()
{
    // Opens the ascii debugging output and writes the headerline.
    ofstream *fout=NULL;
    if (fNameAsciiOffset)
        fout = new ofstream(fNameAsciiOffset, ios::out);
    if (fout)             // changed *fout
        if (isPulserFile)
        {
            *fout //<< "Layer     Cell      crosspoint" << endl;
                << setw(3)  << "s"              << " "
                << setw(4)  << "mo"             << " "
                << setw(4)  << "mb"             << " "
                << setw(5)  << "t"              << " "
		<< setw(10) << "yequalzero"     << " "
		<< setw(10) << "Offset"         << " "
                << setw(15) << "GaussMean"      << " "
                << setw(15) << "GaussSigma"     << " "
                << setw(15) << "fitpar0"        << " "
                << setw(15) << "fitpar1"        << " "
                << setw(15) << "fitparNoise0"   << " "
                << setw(15) << "fitparNoise1"   << " "
                << setw(15) << "totalsigma"     << " "
                << setw(10) << "Integral"       << " "<< endl;
          *fout << setw(14) << "[MdcCalParRaw]" <<endl;
        }
        else
        {
            *fout //<< "Layer     Cell      crosspoint" << endl;
                << setw(3)  << "s"              << " "
                << setw(4)  << "mo"             << " "
                << setw(4)  << "mb"             << " "
                << setw(5)  << "t"              << " "
		<< setw(10) << "yequalzero"     << " "
		<< setw(10) << "Offset"         << " "
                << setw(15) << "fitpar0"        << " "
                << setw(15) << "fitpar1"        << " "
                << setw(15) << "fitparNoise0"   << " "
                << setw(15) << "fitparNoise1"   << " "
                << setw(15) << "totalsigma"     << " "
                << setw(10) << "Integral"       << " "<< endl;
          *fout << setw(14) << "[MdcCalParRaw]" <<endl;
        }
    return fout;
}
void HMdcOffset::fillNTuples(Int_t s,Int_t mo,Int_t mb,Int_t t)
{
    // fill NTuples for offsets and multiple peaks
    HMdcLookupChan& chan=(*lookupgeom)[s][mo][mb][t];
    offsetTuple->Fill(s,mo,mb,t,
		      chan.getNLayer(),
		      chan.getNCell(),
		      (*calparraw)[s][mo][mb][t].getSlope(),
		      (*calparraw)[s][mo][mb][t].getSlopeErr(),
		      offsets[s][mo][mb][t],
		      offsetErr[s][mo][mb][t],
		      integral[s][mo][mb][t],
		      offset1[s][mo][mb][t],
		      fitslope1[s][mo][mb][t],
		      fitslope2[s][mo][mb][t]);

    offsetPulserTuple->Fill(s,mo,mb,t,
			    chan.getNLayer(),
			    chan.getNCell(),
			    offsets[s][mo][mb][t],
			    offsetErr[s][mo][mb][t],
			    integral[s][mo][mb][t],
			    offsetpulser[s][mo][mb][t][0],
			    offsetpulser[s][mo][mb][t][1],
			    offsetpulser[s][mo][mb][t][2],
			    offsetpulser[s][mo][mb][t][3],
			    offsetpulser[s][mo][mb][t][4]);

}
void HMdcOffset::fillArrays(TH1F *hOffset, Int_t s,Int_t mo,Int_t mb,Int_t t)
{
    if(!isPulserFile)
    {
	hOffset->Fill(crosspointX);
	offsets[s][mo][mb][t]=crosspointX;
	offsetErr[s][mo][mb][t]=totalsigma;
	offset1[s][mo][mb][t]=-yequalzero;
    }
    else
    {
	hOffset->Fill(fitGaussMean);
	offsets[s][mo][mb][t]=-fitGaussMean;
	offsetErr[s][mo][mb][t]=fitGaussSigma;
	offset1[s][mo][mb][t]=-yequalzero;
    }
}
void HMdcOffset::fillCalParRaw(TH1F *hOffsetcorr, Int_t s,Int_t mo,Int_t mb,Int_t t)
{
    // fills CalParRaw with offsets, offseterrors and Methods
    // Checks if offset is in valid range arround mean value of offsets.
    // If it is out of range the value is replaced by the mean value.
    // Not connected channels are marked with -7.
    myoffset=offsets  [s][mo][mb][t];
    myerror =offsetErr[s][mo][mb][t];
    HMdcLookupChan& chan=(*lookupgeom)[s][mo][mb][t];
    Int_t connected=0;
    connected=chan.getNCell();
    if(!myoffset){myoffset=0;}

    if(myoffset>=0)   // all offsets where no error showed up
    {
	if((myoffset<meanhOffset-validRange) ||
	   (myoffset>meanhOffset+validRange))
	{   // if offset out of alowed range replace it by mean

	    (*calparraw)[s][mo][mb][t].setOffset(meanhOffset);
	    (*calparraw)[s][mo][mb][t].setOffsetMethod(0);
	    (*calparraw)[s][mo][mb][t].setOffsetErr(100);
	    offsets[s][mo][mb][t]=meanhOffset;

	    hOffsetcorr->Fill(meanhOffset);
	}
	else  // if offsets in allowed range and no error flag
	{
	    hOffsetcorr->Fill(myoffset);
	    (*calparraw)[s][mo][mb][t].setOffset(myoffset);
	    (*calparraw)[s][mo][mb][t].setOffsetMethod(2);
	    if(myerror>99) // if some nonsense value shows up in Err
	    {
		(*calparraw)[s][mo][mb][t].setOffsetErr(100);
		offsetErr[s][mo][mb][t]=100;
	    }
	    else // if Err gives nomal values
	    {
		(*calparraw)[s][mo][mb][t].setOffsetErr(myerror);
	    }
	}
    }
    else
    {  // if offsets contain error flags
	if(connected==-1) // if no wire is connected
	{
	    (*calparraw)[s][mo][mb][t].setOffset(-7);
	    offsets[s][mo][mb][t]=-7;
	    hOffsetcorr->Fill(-7);
	}
	else // if wire is connected but error showed up
	{
	    (*calparraw)[s][mo][mb][t].setOffset(myoffset);
	    hOffsetcorr->Fill(myoffset);
	}
	(*calparraw)[s][mo][mb][t].setOffsetMethod(2);
	(*calparraw)[s][mo][mb][t].setOffsetErr(100);
	offsetErr[s][mo][mb][t]=100;
    }
}

Bool_t HMdcOffset::finalize(void)
{
    // This function is called after the execute function is finished. At this point
    // the arrays for the drift-time are filled. Froms this arrays the histograms for
    // the drift-time and the integrated drift-time are filled. The fits for the calcutation
    // of the offsets are done and the offsets and the sigma of the offsets are
    // calculated .All histograms are written to a rootfile output and the information
    // of the fits to an ascii debugging output. The offset and the sigma of the offset
    // are filled into the container (HMdcCalParRaw) for the calibration parameters.
    // To characterize the method of calculation of the offsets the function setOffsetMethod()
    // of HMdcCalParRaw is used. As default value for the method 2 is used for  automatic
    // calibration. If no calibration is done this value is set to 0.
    cout << "From HMdcOffset::finalize" << endl;

    // Start Stop Watch
    TStopwatch timer;
    timer.Start();

    // open Root file for writing
    htime1temp = new MyHist("time1temp","time1temp", nbin,-nbin+0.5,0.5); //temp time1 hist.

    TFile *file = NULL;
    offsetTuple=NULL;
    offsetPulserTuple=NULL;
    ferrorlog=fopen("offset_error_log.txt","w"); // open file for log of errors/warnings

    if (fNameRootOffset)
    {
	file = new TFile(fNameRootOffset,"RECREATE");
	file->cd();
	offsetTuple=new TNtuple("offset","offset","s:m:mb:tdc:lay:cell:slope:slopeErr:off:offErr:integral:off1:fit1slope:fit2slope");
	offsetPulserTuple=new TNtuple("multipeaks","multipeaks","s:m:mb:tdc:lay:cell:off:offErr:integral:p1:p2:p3:p4:p5");


	// open ascii file for writing
        ofstream *fout=NULL;
	if (fNameAsciiOffset)
        {
	    fout = openAsciiFile();
        }
	TH1F *hOffset = new TH1F("Offset", "Offset", 512, 0, 2048);
	TH1F *hOffsetcorr = new TH1F("Offsetcorr", "Offsetcorr", 512, 0, 2048);

	Int_t err=-3;

	meanhOffset=0;

        initArrays();

	HDetector *mdcDet = gHades->getSetup()->getDetector("Mdc");
	if (!mdcDet)
	    cout << "Error in HMdcOffset::finalize: Mdc-Detector setup (gHades->getSetup()->getDetector(\"Mdc\")) missing." << endl;
	else
	    for(Int_t s=0; s<6; s++) {  //loop over sectors
		cout<<"Sector "<<s<<endl;
		TDirectory *dirSec=NULL;
		if (file) dirSec=Mkdir(file, "sector", s);
		for(Int_t mo=0; mo<4; mo++) {  //loop over modules
		    cout<<"Module "<<mo<<endl;
		    if (!mdcDet->getModule(s, mo)) continue;
		    TDirectory *dirMod=NULL;
		    if (dirSec) dirMod=Mkdir(dirSec, "module", mo);
                    createHist_2D(s,mo);
		    Int_t nMb=(*calparraw)[s][mo].getSize();
		    for(Int_t mb=0; mb<nMb; mb++) {  //loop over layer
			TDirectory *dirMbo=NULL;
			if (dirMod)
			{
			    dirMbo=Mkdir(dirMod, "mbo", mb, 2);
			    dirMbo->mkdir("hinv");
			    dirMbo->mkdir("hint");
			}
			Int_t nTdc=(*calparraw)[s][mo][mb].getSize();
			for(Int_t t=0; t<nTdc; t++) {  //loop over cells
			    createHist(file, s, mo, mb, t);
			    fillHist  (s, mo, mb, t);
			    if ((err=fitHist(s, mo, mb, t))) crosspointX = fitpar1 = -err;
                            fillArrays(hOffset,s,mo,mb,t);
                            fillHist_2D(s,mo,mb,t);
			    if (file)  writeHist (file);
			    if (fNameAsciiOffset){writeAscii(*fout, s, mo, mb, t);}
			    deleteHist();
			}
			if (dirMbo) dirMbo->TDirectory::Cd("..");
			cout << "." << flush;
		    }
		    writeHist_2D();
		    deleteHist_2D();
		    cout<<" "<<endl;
		    if (dirMod) dirMod->TDirectory::Cd("..");
		}
		cout<<" "<<endl;
		if (dirSec) dirSec->TDirectory::Cd("..");
	    }
	cout << endl;

	// write offsets to calparraw
	myoffset=0;
	myerror=0;

	meanhOffset=(Float_t)((Int_t)(hOffset->GetMean()));

	for(Int_t s=0; s<6; s++)
	{  //loop over sectors
	    for(Int_t mo=0; mo<4; mo++)
	    {  //loop over modules
		if (!mdcDet->getModule(s, mo)) continue;
		Int_t nMb=(*calparraw)[s][mo].getSize();
		for(Int_t mb=0; mb<nMb; mb++)
		{  //loop over layer
		    Int_t nTdc=(*calparraw)[s][mo][mb].getSize();
		    for(Int_t t=0; t<nTdc; t++)
		    {  //loop over cells
                        fillCalParRaw(hOffsetcorr,s,mo,mb,t);
			fillNTuples(s,mo,mb,t);
		    }
		}
	    }
	}

	// stop watch
	cout << "Time of cell loop (offset calculation time):" << endl;
	timer.Print();
	timer.Stop();

	hOffset->Write();
	hOffsetcorr->Write();
	offsetTuple->Write();
        offsetPulserTuple->Write();

	delete hOffset;
	delete hOffsetcorr;

	// close root file
	file->Save();
	file->Close();
	delete file;

	if (fNameAsciiOffset)
	{ // close ascii file
	    fout->close();
	    delete fout;
	}
	else
	{
	    cout<<"WARNING: NO ASCII OUTPUT SET, NO FILE WRITTEN TO DISK!"<<endl;
            fprintf(ferrorlog,"%s\n","WARNING: NO ASCII OUTPUT SET, NO FILE WRITTEN TO DISK!");
	}
    }
    else
    {
	cout<<"WARNING: NO ROOT OUTPUT SET, NO FILE WRITTEN TO DISK!"<<endl;
        fprintf(ferrorlog,"%s\n","WARNING: NO ROOT OUTPUT SET, NO FILE WRITTEN TO DISK!");
    }
    
    fclose(ferrorlog);
    return kTRUE;
}

TDirectory *HMdcOffset::Mkdir(TDirectory *dirOld, Char_t *c, Int_t i, Int_t p)    //! Makes new Dir, changes to Dir, returns pointer to new Dir
{
    // Function to create subdirectories
    static Char_t sDir[255];// = new Char_t[strlen(c)+3];
    static Char_t sFmt[10];
    sprintf(sFmt, "%%s %%0%1ii", p);
    sprintf(sDir, sFmt, c, i);
    TDirectory *dirNew = dirOld->mkdir(sDir);
    dirOld->cd(sDir);
    //  free (sDir);
    return dirNew;
}

Float_t HMdcOffset::getstarttime(){
    // Need some work for multiple hists in start detector
    // Better select multiplicity 1 in start.
    Int_t i=0;
    Int_t startmod=-1;
    HStartHit* starthit=0;
    iter_start->Reset();
    Float_t starttime=0;
    while ((starthit=(HStartHit *)iter_start->Next())!=0) {
	startmod=starthit->getModule();
	if(startmod==0)
	{
	    i++;
	    if(starthit->getFlag()) starttime=starthit->getTime();
	}
    }
    if (i!=1) Error("getstarttime(int)","Multiplicity in Start > 1 or 0");
    return starttime;
}

Int_t HMdcOffset::execute()
{
    Int_t result=0;
    if(useClusters)result=executeClus();
    else           result=executeRaw();
    return result;
}
Int_t HMdcOffset::executeClus()
{
    // Fired cells are take from the list of Cells inside one HMdcClus.
    // Raw data of Time1  multiplied by the slope of the channel taken from the
    // container of the calibration parameters are filled into the array for the
    // drift-time . This array contains 2048 (corresponding to the resulution of
    // the tdc chip) entries for each tdc channel.
    static HMdcRaw *raw  =NULL;
    static HMdcClus *clus=NULL;
    Int_t s, mo, mb, t;
    Int_t l,c;

    Float_t testTime1=0;
    Float_t testTime2=0;

    Double_t xp1,yp1,zp1,xp2,yp2,zp2,signalpath;
    Double_t wireoff =0.;
    Float_t starttime=0;
    Float_t tof      =0;

    iter_clus->Reset();
    if(skipcounter>nSkipEvents)
    {
	if(!noStart) // Start Time is used
	{
	    starttime=getstarttime();
	}
	while ((clus=(HMdcClus*)iter_clus->Next()))
	{
	    s=mo=mb=t==-99;
	    l=c=-99;
	    xp1=clus->getX();
	    yp1=clus->getY();
	    zp1=clus->getZ();
	    xp2=clus->getXTarg();
	    yp2=clus->getYTarg();
	    zp2=clus->getZTarg();

	    s =clus->getSec();
	    locraw[0]=s;

	    Int_t ioseg=clus->getIOSeg();
	    Int_t mod=clus->getMod();

	    for(Int_t lay=0;lay<12;lay++)
	    {  // loop over layers
		//---------------finding module-----------
		if(mod<0)
		{ 	// combined clusters
		    if(ioseg==0){
			(lay<6)? locraw[1]=0 : locraw[1]=1;
		    }else if(ioseg==1){
			(lay<6)? locraw[1]=2 : locraw[1]=3;
		    }
		}
		else
		{  // chamber clusters
		    locraw[1]=mod;
		}
		mo=locraw[1];
		//---------------finding layer------------
		(lay<6)? l=lay : l=lay-6;
		//----------------------------------------

		Int_t nCells=clus->getNCells(lay);
		for(Int_t cell=0;cell<nCells;cell++)
		{  // loop over cells in layer
		    c=clus->getCell(lay,cell);

		    if(useWireOffset)
		    {
			signalpath=(*sizescells)[s][mo][l].getSignPath(xp1,yp1,zp1,xp2,yp2,zp2,c);
			if(signalpath>0) wireoff=signalpath*signalspeed;
		    }

		    HMdcLookupCell& mycell=(*lookupraw)[s][mo][l][c];
		    mb=mycell.getNMoth();
		    t =mycell.getNChan();
		    locraw[2]=mb;
		    locraw[3]=t;
		    raw=(HMdcRaw*)rawCat->getObject(locraw);

		    if(((*calparraw)[s][mo][mb][t].getSlope())!=0)
		    {
			testTime1=raw->getTime(1)-(starttime/(*calparraw)[s][mo][mb][t].getSlope());
			testTime2=raw->getTime(2)-(starttime/(*calparraw)[s][mo][mb][t].getSlope());
		    }
		    else
		    {
			testTime1=raw->getTime(1);
			testTime2=raw->getTime(2);
		    }

                    if(useTof) tof=toffunc[mo][l]->Eval(c); // if useTof is switched on

		    if(useTimeCuts)  // if Time Cuts and Start Time
		    {
			if(testTimeCuts(s,mo,testTime1,testTime2))
			{
			    // we have to use Float_ts here
			    (*hreverse)[s][mo][mb][t]
				[
				 (Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope())+starttime-wireoff-tof)
				]++;
			}
		    }
		    else
		    {   // No Time Cuts and Start Time
			(*hreverse)[s][mo][mb][t]
			    [
			     (Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope())+starttime-wireoff-tof)
			    ]++;
		    }
		}
	    }
	}
	if(eventcounter%nStep==0){
	    cout<<eventcounter<<" events analyzed"<<endl;
	}
    }
    skipcounter++;
    eventcounter++;
    return 0;
}
Int_t HMdcOffset::executeRaw()
{
    // Raw data of Time1  multiplied by the slope of the channel taken from the
    // container of the calibration parameters are filled into the array for the
    // drift-time . This array contains 2048 (corresponding to the resulution of
    // the tdc chip) entries for each tdc channel.
    static HMdcRaw *raw=NULL;
    Int_t s, mo, mb, t;
    Float_t testTime1=0;
    Float_t testTime2=0;
    Float_t starttime=0.;
    Float_t tof      =0.;

    iter->Reset();
    if(skipcounter>nSkipEvents)
    {
	if(!noStart) // Start Time is used
	{
	    starttime=getstarttime();
	}
	while ((raw=(HMdcRaw*)iter->Next()))
	{
	    s=mo=mb=t=-99;
	    raw->getAddress(s, mo, mb, t);
	    if(((*calparraw)[s][mo][mb][t].getSlope())!=0)
	    {
		testTime1=raw->getTime(1)-(starttime/(*calparraw)[s][mo][mb][t].getSlope());
		testTime2=raw->getTime(2)-(starttime/(*calparraw)[s][mo][mb][t].getSlope());
	    }
	    else
	    {
		testTime1=raw->getTime(1);
		testTime2=raw->getTime(2);
	    }

	    if(useTof)
	    {
		HMdcLookupChan& chan=(*lookupgeom)[s][mo][mb][t];
		Int_t l= chan.getNLayer();
		Int_t c= chan.getNCell();
		tof=toffunc[mo][l]->Eval(c);
	    }

	    if(useTimeCuts)  // if Time Cuts and Start Time
	    {
		if(testTimeCuts(s,mo,testTime1,testTime2))
		{
		    // we have to use Float_ts here
		    (*hreverse)[s][mo][mb][t]
			[
			 (Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope())+starttime-tof)
			]++;
		}
	    }
	    else
	    {   // No Time Cuts and Start Time
		(*hreverse)[s][mo][mb][t]
		    [
		     (Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope())+starttime-tof)
		    ]++;
	    }
	}
	if(eventcounter%nStep==0){
	    cout<<eventcounter<<" events analyzed"<<endl;
	}
    }
    skipcounter++;
    eventcounter++;
    return 0;
}