//*-- AUTHOR : J.Kempter
//*-- Modified :
//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////////////////////////
//
// HMdcOffsetCheckTimeShift
//
// 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 <TStopwatch.h>
#include <TDirectory.h>
#include "heventheader.h"
#include "hmdcoffsetchecktimeshift.h"
#include "hmdcdef.h"
#include "hmdcraw.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 "hmdccalpar.h"
#include "hmdccalparraw.h"
#include "hmdclookupraw.h"
#include "hstartcal.h"
#include "hmdctimecut.h"
ClassImp(HMdcOffsetCheckTimeShift)
//const Int_t HMdcOffsetCheckTimeShift::cells_Num[] ={184, 192, 172, 172, 192, 192};
const Int_t HMdcOffsetCheckTimeShift::nbin = 2048;
const Int_t HMdcOffsetCheckTimeShift::nbinp1 = HMdcOffsetCheckTimeShift::nbin+1;
const Int_t HMdcOffsetCheckTimeShift::nbinm1 = HMdcOffsetCheckTimeShift::nbin-1;
const Int_t HMdcOffsetCheckTimeShift::nSubEvents = (6*192)*3;
HMdcOffsetCheckTimeShift::HMdcOffsetCheckTimeShift(void)
{
setDefault();
rawCat = NULL;
calStartCat=NULL;
iter = NULL;
iter_start= NULL;
isPulserFile = kFALSE;
noStart=kFALSE;
useTimeCut=kTRUE;
hintegral = (MyFieldMdc*) new MyFieldMdc;
hreverse = (MyFieldMdc*) new MyFieldMdc;
lookupgeom=0;
timecut=0;
filecounter=0;
skipcount=0;
stepsize=1000;
step=0;
numberofevents=1000;
}
HMdcOffsetCheckTimeShift::HMdcOffsetCheckTimeShift(Text_t* name,Text_t* title) : HReconstructor(name,title)
{
setDefault();
rawCat = NULL;
calStartCat=NULL;
iter = NULL;
iter_start= NULL;
isPulserFile = kFALSE;
noStart=kFALSE;
useTimeCut=kTRUE;
hintegral = (MyFieldMdc*) new MyFieldMdc;
hreverse = (MyFieldMdc*) new MyFieldMdc;
lookupgeom=0;
timecut=0;
filecounter=0;
skipcount=0;
stepsize=1000;
step=0;
numberofevents=1000;
}
HMdcOffsetCheckTimeShift::~HMdcOffsetCheckTimeShift(void)
{
// Destructor.
if (iter) delete iter;
if (fNameAsciiOffset) delete fNameAsciiOffset;
if (fNameRootOffset) delete fNameRootOffset;
if (hintegral) delete[] hintegral;
if (hreverse) delete[] hreverse;
iter = NULL;
iter_start= NULL;
}
void HMdcOffsetCheckTimeShift::setDefault()
{
// Sets the the default values for the fits.
setNoise(100, 50);
setThreshold(0.15, 0.5);
setRangeGauss(50);
}
void HMdcOffsetCheckTimeShift::setOutputAscii(Char_t *c)
{
// Sets ascii output of HMdcOffsetCheckTimeShift for debugging output.
if (fNameAsciiOffset) delete fNameAsciiOffset;
fNameAsciiOffset = new Char_t[strlen(c)+1];
strcpy(fNameAsciiOffset, c);
}
void HMdcOffsetCheckTimeShift::setOutputRoot(Char_t *c)
{
// Sets rootfile output of HMdcOffsetCheckTimeShift where all created histograms were written.
//
if (fNameRootOffset) delete fNameRootOffset;
fNameRootOffset = new Char_t[strlen(c)+1];
strcpy(fNameRootOffset, c);
}
Bool_t HMdcOffsetCheckTimeShift::init(void)
{
// Inits HMdcOffsetCheckTimeShift and the needed HMdcCalParRaw, if this container does not exists
if(!noStart){
calStartCat=gHades->getCurrentEvent()->getCategory(catStartCal); //????
if (!calStartCat) {
calStartCat=gHades->getSetup()->getDetector("Start")->buildCategory(catStartCal);
if (!calStartCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catStartCal,calStartCat,"Start");
}
iter_start=(HIterator *)calStartCat->MakeIterator("native");
}
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) {
rawCat=gHades->getSetup()->getDetector("Mdc")->buildCategory(catMdcRaw);
if (!rawCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catMdcRaw,rawCat,"Mdc");
}
// creates an iterator which loops over all fired cells
iter=(HIterator *)rawCat->MakeIterator();
fActive=kTRUE;
eventcounter=0;
step=0;
return kTRUE;
}
Bool_t HMdcOffsetCheckTimeShift::reinit(void)
{
// Inits HMdcOffsetCheckTimeShift and the needed HMdcCalParRaw, if this container does not exists
filecounter++;
cout<<"file "<<filecounter<<endl;
strcpy(filenames[filecounter-1],
strrchr((Text_t*)gHades->getRuntimeDb()->getCurrentFileName(),47)+1);
cout<<filenames[filecounter-1]<<endl;
skipcount=0;
step=0;
fActive=kTRUE;
return kTRUE;
}
void HMdcOffsetCheckTimeShift::createHist(TFile* file,Int_t s, Int_t m, Int_t step,Int_t f)
{
// 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%02i","Sector ",s," Module ",m," Step ",step," File ",f);
file->TDirectory::Cd("hinv");
sprintf(tmp, "%s%i%s%i%s%02i%s%02i%s", "hinv[",s,"][",m,"][",step,"][",f,"]");
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,"][",step,"][",f,"]");
hint = new MyHist(tmp, title, nbin,-nbin+0.5,0.5); // integrated hist.
file->TDirectory::Cd("..");
}
void HMdcOffsetCheckTimeShift::deleteHist()
{
// Created histograms are deleted
delete hinv;
delete hint;
}
void HMdcOffsetCheckTimeShift::fillHist(Int_t s, Int_t m, Int_t step,Int_t f)
{
// Histograms for inverted Time1 and integrated Time1 per Tdc-Channel are filled
MyInt *hi = &((*hintegral)[s][m][step][f][0]);
MyInt *hr = &((*hreverse) [s][m][step][f][0]);
*hi = *hr;
hint->SetBinContent(1, *hi++);
hinv->SetBinContent(1, *hr++);
for(Int_t k=2; k<nbinp1; k++) {
*hi=*(hi-1) + *hr;
hint->SetBinContent(k, *hi++);
hinv->SetBinContent(k, *hr++);
}
}
Int_t HMdcOffsetCheckTimeShift::fitHist(Int_t s, Int_t m, Int_t step,Int_t file)
{
// The offset is calculated from two linear fits. Some pictures of the fitted histograms
// can be found on
MDC calibration page
.
// 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 calculated
// offset.The mean and sigma of the fit is written to the ascii debugging output.
fitGaussMean = 0;
fitGaussSigma = 0;
crosspointX=0;
yequalzero=0;
totalsigma=0;
fitpar0 =0;
fitpar0error=0;
fitpar1 =0;
fitpar1error=0;
fitparNoise0 =0;
fitparNoise1 =0;
fitparNoise0error=0;
fitparNoise1error=0;
MyInt *hI = &((*hintegral)[s][m][step][file][0]);
Int_t xmaxfitbin=nbinm1;
Float_t yminfit=minfitthreshold*hI[nbinm1-1]; // setting y thresholds for fitting
Float_t ymaxfit=maxfitthreshold*hI[nbinm1-1]; // relative to maximum y-value
while (hI[--xmaxfitbin] > ymaxfit && xmaxfitbin); // find to ymaxfit corresponding bin
if (!xmaxfitbin) return 1;
Int_t xminfitbin=xmaxfitbin;
while (hI[--xminfitbin] > yminfit && xminfitbin); // find to yminfit corresponding bin
if (!xminfitbin) return 2;
TF1 *f = new TF1("f1","pol1",hint->GetBinCenter(xminfitbin),hint->GetBinCenter(xmaxfitbin)); // call fit function poynom 1.order in x-range
hint->Fit("f1","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]
delete f;
if (fitpar1==0) return 3;
yequalzero=0;
yequalzero =-fitpar0/fitpar1;
if (-yequalzero < (offsetfitNoise+widthfitNoise) || -yequalzero > nbin) return 4;
f = new TF1("f2","pol1",
yequalzero-(offsetfitNoise+widthfitNoise),
yequalzero-offsetfitNoise);
hint->Fit("f2","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);
delete f;
if (fitpar1==fitparNoise1) return 5;
crosspointX=0;
crosspointX=-(fitparNoise0-fitpar0)/(fitpar1-fitparNoise1);
if (crosspointX<0 || crosspointX>nbinm1 || isNaN(crosspointX)) return 6;
// calculating total sigma
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))
);
if (!isPulserFile) return 0;
f = new TF1("f3","gaus", -crosspointX-rangeGauss, -crosspointX+rangeGauss);
hinv->Fit("f3","RNQ");
fitGaussMean =f->GetParameter(1); //mean
fitGaussSigma=f->GetParameter(2); //sigma
delete f;
// FIXME: Error handling for gauss fit ???
return 0;
}
void HMdcOffsetCheckTimeShift::writeAscii(ofstream &fout, Int_t s, Int_t m, Int_t step,Int_t f)
{
// 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) << step << " "
<< setw(5) << f << " "
<< setw(10) << -yequalzero << " "
<< setw(10) << crosspointX << " "
<< setw(15) << -fitGaussMean << " "
<< setw(15) << fitGaussSigma << " "
<< setw(25) <<filenames[f] << " "
<< endl;
else
fout
<< setw(3) << s << " "
<< setw(4) << m << " "
<< setw(4) << step << " "
<< setw(5) << f << " "
<< setw(10) << -yequalzero << " "
<< setw(10) << crosspointX << " "
<< setw(10) << totalsigma << " "
<< setw(25) <<filenames[f] << " "
<< endl;
}
void HMdcOffsetCheckTimeShift::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("..");
}
ofstream *HMdcOffsetCheckTimeShift::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) << "step" << " "
<< setw(5) << "file" << " "
<< setw(10) << "yequalzero" << " "
<< setw(10) << "Offset" << " "
<< setw(15) << "GaussMean" << " "
<< setw(15) << "GaussSigma" << " "
<< setw(25) << "filename" << endl;
*fout << setw(14) << "[MdcCalParRaw]" <<endl;
}
else
{
*fout //<< "Layer Cell crosspoint" << endl;
<< setw(3) << "s" << " "
<< setw(4) << "mo" << " "
<< setw(4) << "step" << " "
<< setw(5) << "file" << " "
<< setw(10) << "yequalzero" << " "
<< setw(10) << "Offset" << " "
<< setw(10) << "Error" << " "
<< setw(25) << "filename" <<endl;
*fout << setw(14) << "[MdcCalParRaw]" <<endl;
}
return fout;
}
Bool_t HMdcOffsetCheckTimeShift::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 HMdcOffsetCheckTimeShift::finalize" << endl;
// Start Stop Watch
TStopwatch timer;
timer.Start();
// open Root file for writing
TFile *file = NULL;
if (fNameRootOffset)
file = new TFile(fNameRootOffset,"RECREATE");
// open ascii file for writing
ofstream *fout = openAsciiFile();
Int_t err;
Float_t offsets[6][4][20][50];
Float_t offsetErr[6][4][20][50];
HDetector *mdcDet = gHades->getSetup()->getDetector("Mdc");
if (!mdcDet)
cout << "Error in HMdcOffsetCheckTimeShift::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);
if (dirMod)
{
dirMod->mkdir("hinv");
dirMod->mkdir("hint");
}
for(Int_t f=0;f<50;f++){ // loop over files
for(Int_t step=0;step<20;step++){ // loop over steps per file
createHist(file, s, mo,step,f);
fillHist (s, mo,step,f);
if ((err=fitHist(s, mo,step,f))) crosspointX = fitpar1 = -err;
offsets[s][mo][step][f]=crosspointX;
offsetErr[s][mo][step][f]=totalsigma;
if (file) writeHist (file);
if (*fout) writeAscii(*fout, s, mo,step,f);
deleteHist();
}
cout<<"."<<flush;
}
cout<<" "<<endl;
if (dirMod) dirMod->TDirectory::Cd("..");
}
cout<<" "<<endl;
if (dirSec) dirSec->TDirectory::Cd("..");
}
cout << endl;
// write offsets to calparraw
// stop watch
cout << "Time of cell loop (offset calculation time):" << endl;
timer.Print();
timer.Stop();
// close root file
file->Save();
file->Close();
delete file;
// close ascii file
fout->close();
delete fout;
return kTRUE;
}
TDirectory *HMdcOffsetCheckTimeShift::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 HMdcOffsetCheckTimeShift::getstarttime(){
Int_t i=0;
HStartCal* startcal=0;
iter_start->Reset();
Float_t starttime=0;
while ((startcal=(HStartCal *)iter_start->Next())!=0) {
switch (startcal->getStrip()){
case 0: starttime=(startcal->getTime());
break;
case 1: starttime=(startcal->getTime());
break;
case 2: starttime=(startcal->getTime());
break;
case 3: starttime=(startcal->getTime());
break;
case 4: starttime=(startcal->getTime());
break;
case 5: starttime=(startcal->getTime());
break;
case 6: starttime=(startcal->getTime());
break;
case 7: starttime=(startcal->getTime());
break;
}
i++;
}
if (i!=1)cout<<"Multiplicity in Start > 1"<<endl;
return starttime;
}
Int_t HMdcOffsetCheckTimeShift::execute()
{
// 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;
skipcount++; // starts with 0 for each file
step=(Int_t)(skipcount/stepsize);
if(eventcounter%10000==0){
cout<<eventcounter<<" events analyzed"<<endl;
}
eventcounter++; // counts number of calculated events of all files
if(step>19)return 0;
if(skipcount>(step*stepsize+numberofevents))return 0;
if(skipcount%1000==0)cout<<skipcount<<endl;
iter->Reset();
if(!noStart)
{
Float_t 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(useTimeCut)
{
if(timecut->cutTime1(s,mo,testTime1)&&
timecut->cutTime2(s,mo,testTime2)&&
timecut->cutTimesDif(s,mo,testTime2,testTime1))
{
(*hreverse)[s][mo][step][filecounter-1]
[
(Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope())+starttime)
]++;
}
}
else
{
(*hreverse)[s][mo][step][filecounter-1]
[
(Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope())+starttime)
]++;
}
}
}
else
{
while ((raw=(HMdcRaw*)iter->Next()))
{
s=mo=mb=t=-99;
raw->getAddress(s, mo, mb, t);
testTime1=raw->getTime(1);
testTime2=raw->getTime(2);
if(useTimeCut)
{
if(timecut->cutTime1(s,mo,testTime1)&&
timecut->cutTime2(s,mo,testTime2)&&
timecut->cutTimesDif(s,mo,testTime2,testTime1))
{
// we have to use Float_ts here
(*hreverse)[s][mo][step][filecounter-1]
[
(Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope()))
]++;
}
}
else
{
(*hreverse)[s][mo][step][filecounter-1]
[
(Int_t)(2048-((Float_t)(raw->getTime(1))*(*calparraw)[s][mo][mb][t].getSlope()))
]++;
}
}
}
return 0;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.