Commit 73b3f26e authored by Vratislav Chudoba's avatar Vratislav Chudoba

Source for libraries upload.

parent b50c7064
//////////////////////////////////////////////////////////
// //
// AculCalibration //
// //
//
//Some description of this very useful class,
//its properties and a short example how to
//use it. This text may be partly used in
//PhD thesis.
//
//Description of the detector itself.
//
// //
//////////////////////////////////////////////////////////
#include "AculCalibration.h"
ClassImp(AculCalibration);
AculCalibration::AculCalibration()
{
//default constructor
fCurrentHStack = NULL;
fCurrentHistList.IsOwner();
kRaNOPEAKS = 0;
fLowerPeakRelativeHight = 0.;
fUpperPeakRelativeHight = 0.;
fPeakPositionTolerance = 0.;
fFitFuncLineWidth = 1;
fFitMinSigma = 0.;
fFitPeakThreshold = 0.;
for(Int_t i = 0; i < DEFAULTNOPEAKS; i++) {
fEnergy[i] = 0;
}
fCalInformation = 0;
Reset();
}
AculCalibration::AculCalibration(const char* parfile)
{
//constructor which fills fA, fB, fC, fD from file parfile
fCurrentHStack = NULL;
fCurrentHistList.IsOwner();
kRaNOPEAKS = 0;
fLowerPeakRelativeHight = 0.;
fUpperPeakRelativeHight = 0.;
fPeakPositionTolerance = 0.;
fFitFuncLineWidth = 1;
fFitMinSigma = 0.;
fFitPeakThreshold = 0.;
for(Int_t i = 0; i < DEFAULTNOPEAKS; i++) {
fEnergy[i] = 0;
}
fCalInformation = 0;
SetCalibrationParameters(parfile);
}
AculCalibration::~AculCalibration()
{
DeleteStacks();
// delete fCalInformation;
// fCalInformation->Close();
}
Int_t AculCalibration::SearchPeaks(const TH1 *hin, Double_t sigma, Option_t *option, Double_t threshold, const Int_t searchedpeaks)
{
//Function searching peaks in inputed TH1 spectrum and selects the peaks in the histogram.
//
// hin:
// sigma:
// option:
// threshold:
// searchedpeaks:
TSpectrum sc; //by default for 100 peaks
Int_t nopeaks = sc.Search(hin, sigma, "goff", threshold);
TString opt = option;
opt.ToLower();
const Double_t tStep = 0.05;
while ( nopeaks > searchedpeaks && threshold <= 1) {
threshold = threshold + tStep;
nopeaks = sc.Search(hin, sigma, "goff", threshold);
}
if (!nopeaks) {
return 0;
}
if (opt.Contains("goff")) {
return nopeaks;
}
TPolyMarker *pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
if (pm) {
hin->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new TPolyMarker(nopeaks, sc.GetPositionX(), sc.GetPositionY());
hin->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1.3);
return nopeaks;
}
Int_t AculCalibration::PeaksFitting(TH1I* hSpectrum, Option_t* option, Double_t sigmamin)
{
//function searching peaks in inputed spectrum
//
// hSpectrum:
// option: posible option are
// GP: explanation needed
// WRITEBAD: explanation needed
// Q: displays on the histogram how the spectrum were fitted and doesn't writes out
// the numbers of channels in which are peaks
// V: displays on the histogram how the spectrum were fitted and writes out the numbers
// of channels in which are peaks
// sigmamin:
if (!hSpectrum) return 1;
Int_t dimension = hSpectrum->GetDimension();
if (dimension > 1) {
Error("PeaksFitting", "Only implemented for 1-d histograms");
return 1;
}
TString opt = option;
opt.ToLower();
if (!kRaNOPEAKS) {
Error("PeaksFitting", "kRaNOPEAKS is set to zero; calibration spectrum must be set");
return 1;
}
Int_t peaksNumber = SearchPeaks(hSpectrum, sigmamin, "", fFitPeakThreshold, kRaNOPEAKS);
if (peaksNumber != kRaNOPEAKS) {
Info("PeaksFitting", "In histogram %s was found %d peaks", hSpectrum->GetName(), peaksNumber);
return 1;
}
//predelat jako volitelny vypis
Info("PeaksFitting", "Number of peaks: %d", peaksNumber);
Double_t peak[peaksNumber]; //pracovni pole pro zapis piku v neusporadanem poradi
Double_t *peakPosition;
Double_t *peakHight;
for (Int_t i = 0; i < peaksNumber; i++) {
TList *functions = hSpectrum->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
peakPosition = pm->GetX();
peakHight = pm->GetY();
Double_t currentHight = peakHight[i];
//posouva se smerem nahoru a urcuje hrubou hranici piku, ktera je urcena 1/10 jeho vysky
Int_t j = 0;
Int_t fitMin = 0;
while ( currentHight > (peakHight[i]*fUpperPeakRelativeHight) ) {
j++;
fitMin = static_cast<Int_t>(peakPosition[i]) + j;
currentHight = hSpectrum->GetBinContent(fitMin);
}
//totez, ale urcujeme dolni hranici piku
currentHight = peakHight[i];
j = 0;
Int_t fitMax = 0;
while ( currentHight > (peakHight[i]*fLowerPeakRelativeHight) ) {
j++;
fitMax = static_cast<Int_t>(peakPosition[i]) - j;
currentHight = hSpectrum->GetBinContent(fitMax);
}
//fitting
if (opt.Contains("gp")) {
Info("PeaksFitting", "Option containing gp");
char fncname[20];
sprintf(fncname, "gaus_aux_%d", i);
TF1 *gausAux = new TF1(fncname, "gaus", fitMin - 10, fitMax + 10); //pomocny gaus
hSpectrum->Fit(fncname, "0 Q", "", fitMin - 15, fitMax + 15); //prvotni fitovani
sprintf(fncname, "auto_gp_%d", i);
TF1 *fitAuto = new TF1(fncname, "gaus(0) + pol0(3)", fitMin - 15, fitMax + 15); //fce pro automaticke fitovani
fitAuto->SetParameter(0, gausAux->GetParameter(0)); //nastavovani parametru fitovaci fce
fitAuto->SetParameter(1, gausAux->GetParameter(1));
fitAuto->SetParameter(2, gausAux->GetParameter(2));
hSpectrum->Fit(fncname, "0 R Q +", "", fitMin - 15, fitMax + 15); //dodelat zapis vsech fci
hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw);
peak[i] = fitAuto->GetParameter(1); //zapis asi pozice v kanalech do pomocneho pole
if (opt.Contains("V")) {
Info("PeaksFitting", "Peak position is\t %4.2f \tresolution is \t %2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1)));
}
}
else {
char fncname[20];
sprintf(fncname, "auto_g%d", i);
TF1 *fitAuto = new TF1(fncname, "gaus", fitMin, fitMax); //fce pro automaticke fitovani
fitAuto->SetLineWidth(fFitFuncLineWidth);
hSpectrum->Fit(fncname, "+ 0 R Q", ""/*, fitMin - 1, fitMax + 1*/);
// hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw);
hSpectrum->GetFunction(fncname)->InvertBit(TF1::kNotDraw);
peak[i] = fitAuto->GetParameter(1); //zapis asi pozice v kanalech do pomocneho pole
if (opt.Contains("v")) {
Info("PeaksFitting", "Peak position is\t%4.2f\tresolution is \t%2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1)));
}
}//else
//end of fitting
}//for over all analyzed peaks
//peaks sorting
Int_t j[peaksNumber];
TMath::Sort(peaksNumber, peak, j, kFALSE);
for (Int_t i = 0; i < 4; i++) {
fPeak[i] = peak[j[i]];
}
if (!opt.Contains("q") || opt.Contains("v")) {
Info("PeaksFitting", "Control output:");
for (Int_t i = 0; i < peaksNumber; i++) {
printf("\tPeak position is\t%f\n", fPeak[i]);
}
}
// provest kontrolu pomerne polohy piku,
// jestli jsou spatne, provest urcita opatreni,
// napr. zapis daneho histogramu do souboru,
// zapis do souboru s chybama, vypis na obrazovku, ...
for (Int_t i = 0; i < peaksNumber; i++) {
if ( !( (((1-fPeakPositionTolerance)*(fEnergy[0]/fEnergy[i])) < (fPeak[0]/fPeak[i])) && (((1+fPeakPositionTolerance)*(fEnergy[0]/fEnergy[i])) > (fPeak[0]/fPeak[i])) ) ) {
if (fCalInformation && opt.Contains("writebad")) {
fCalInformation->cd();
hSpectrum->Write();
}
return 2;
}
}//for
return 0;
}
Bool_t AculCalibration::SetInputParameters(const char* inputparfile)
{
// Function which loads the data file for calibration
//
// -inputparfile: file containing information about calibration source
// -In file with the data must be preserved systematic
// -There can't be pause between the lines
//
//For example
//................................................................
//223Ra
//
//4 nopeaks //number of peaks
//4.415 E1 //in MeV
//5.153 E2 //in MeV
//5.683 E3 //in MeV
//7.419 E4 //in MeV
//100 lowerchannel //in channels
//4096 upperchannel //in channels
//0.1 lowerpeakhight //in relative units
//0.1 upperpeakhight //in relative units
//0.1 peakpositiontolerance //in relative units
//2 fitfunctionlinewidth //integer 1 - 10
//5 minfitsigma //minimal sigma of the peaks to be fitted
//0.4 fithightthreshold //
//................................................................
const Int_t lineLength = 400;
Char_t line[lineLength];
Char_t parameter[100];
Char_t identificator[100];
ifstream fipr;
fipr.open(inputparfile);
if (!fipr.is_open()) {
Error("SetInputsParameters", "File with input parameters was not opened");
return kFALSE;
}
while (!fipr.eof()) {
fipr.getline(line, lineLength);
if (strlen(line) == 0) {
continue;
}
sscanf(line, "%s %s", parameter, identificator);
if ( strcmp(identificator, "nopeaks") == 0 ) {
kRaNOPEAKS = static_cast<Int_t>(atoi(parameter));
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
fipr.getline(line, lineLength);
sscanf(line, "%s", parameter);
fEnergy[i] = static_cast<Double_t>(atof(parameter));
}
}//if
if ( strcmp(identificator, "lowerpeakhight") == 0 ) {
sscanf(line, "%s", parameter);
fLowerPeakRelativeHight = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "upperpeakhight") == 0 ) {
sscanf(line, "%s", parameter);
fUpperPeakRelativeHight = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "peakpositiontolerance") == 0 ) {
sscanf(line, "%s", parameter);
fPeakPositionTolerance = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "fitfunctionlinewidth") == 0 ) {
sscanf(line, "%s", parameter);
fFitFuncLineWidth = static_cast<Width_t>(atoi(parameter));
}
if ( strcmp(identificator, "minfitsigma") == 0 ) {
sscanf(line, "%s", parameter);
fFitMinSigma = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "fithightthreshold") == 0 ) {
sscanf(line, "%s", parameter);
fFitPeakThreshold = static_cast<Double_t>(atof(parameter));
}
}
fipr.close();
return kTRUE;
}
void AculCalibration::PrintInputParameters()
{
//print alpha source parameters
cout << "Number of peaks: " << kRaNOPEAKS << endl
<< endl;
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
cout << "fEnergy[" << i << "] = " << fEnergy[i] << endl;
}
return;
}
Bool_t AculCalibration::SetCalibrationParameters(const char* calparfile)
{
//Loads the file with calibration parameters
//
// calparfile: file containing calibration parameters in format: crate number, address, subaddress, fA, fB, fC, fD
const Int_t lineLength = 200;
char line[lineLength];
// Int_t crate;
char crate[100];
Int_t i, j;
char cA[40], cB[40], cC[40], cD[40];
//open file with calibration parameters
ifstream calFileR;
calFileR.open(calparfile);
if( !calFileR.is_open() ) {
Error("SetCalibrationParameters", "File %s with calibration data was not opened", calparfile);
return kFALSE;
}
Reset();
//read calibration parameters from file
while (!calFileR.eof()) {
calFileR.getline(line, lineLength);
if ( line[0] != '*' && line[0] != '#' && line[0] != '%') { //possible comment characters
sscanf(line, "%s %d %d %s %s %s %s", crate, &i, &j, cA, cB, cC, cD);
fA[i][j] = atof(cA);
fB[i][j] = atof(cB);
fC[i][j] = atof(cC);
fD[i][j] = atof(cD);
}
}
calFileR.close();
return kTRUE;
}
void AculCalibration::PrintCalibrationParameters(const Int_t blockmin, const Int_t blockmax)
{
//Print calibration parameters fA, fB, fC, fD
//
// blockmin: minimum block data are displayed for
// blockmax: maximum block data are displayed for
for (Int_t i = blockmin; i <= blockmax; i++) {
for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
cout << "C3[" << i << "][" << j << "]" << setw(10) << fA[i][j] << setw(10) << fB[i][j] << setw(10) << fC[i][j] << setw(10) << fD[i][j] << endl;
}
}
}
void AculCalibration::ShowRawSpectra(const char* filename, const Int_t block, TCanvas* rawCanvas, Int_t xaxismin, Int_t xaxismax, /*TObjArray* histList,*/ const Int_t subaddress)
{
//Displays the spectrum from a file, divides the canvas into a sufficient number of pads and displays spectrums of each
//block subaddress on the suitable pads.
//
// filename: input .root file containing spectra to be showed
// block: block which will be drawn
// rawCanvas: canvas on which one you will see the spectrum
// xaxismin: minimal channel, which will be displayed
// xaxismax: maximal channel, which will be displayed
// subaddress:
Char_t address[40];
Char_t histName[40];
Char_t fillCommand[40];
Char_t fillCondition[40];
if (!rawCanvas) {
//rawCanvas = new TCanvas("RawSpectra", "Raw spectra in channels", 1);
cout << "You have to assign TCanvas for raw spectra drawing" << endl;
return;
}
rawCanvas->Clear();
// cout << "hovno" << endl;
rawCanvas->SetFillColor(10);
// cout << "hovno" << endl;
TFile *fr = new TFile(filename);
if (fr->IsOpen() == 0) {
cout << endl << "File " << filename << " was not opened and won't be processed" << endl << endl;
return;
}
TH1I *hRead = 0;
TTree *tr = (TTree*)fr->Get("RAW");
// cout << "hovno" << endl;
if (subaddress > 15) {
rawCanvas->Divide(4, 4);
rawCanvas->SetFillColor(10);
// cout << "hovno" << endl;
for (Int_t i = 0; i < 16; i++) {
cout << i << endl;
rawCanvas->cd(i+1);
hRead = new TH1I("name", "title", 4096, 0, 4095);
sprintf(address, "C3[%d][%d]", block, i);
sprintf(histName, "H3[%d][%d]", block, i);
hRead->SetName(histName);
sprintf(fillCommand, "%s >> %s", address, hRead->GetName());
sprintf(fillCondition, "%s > 0", address);
tr->Draw(fillCommand, fillCondition, "");
if (hRead) {
hRead->SetDirectory(0);
// cout << hRead->GetEntries() << endl;
// if (fHRawList) {
// fHRawList->Add(hRead);
fHRawList.Add(hRead);
// }
hRead->SetAxisRange(xaxismin, xaxismax);
}
}//for
}
else {
fr->cd();
hRead = new TH1I("name", "title", 4096, 0, 4095);
sprintf(address, "C3[%d][%d]", block, subaddress);
sprintf(histName, "H3[%d][%d]", block, subaddress);
hRead->SetName(histName);
sprintf(fillCommand, "%s >> %s", address, hRead->GetName());
sprintf(fillCondition, "%s > 0", address);
// cout << fillCommand << setw(20) << fillCondition << endl;
tr->Draw(fillCommand, fillCondition, "goff");
if (hRead) {
hRead->SetDirectory(0);
// if (fHRawList) {
// fHRawList->Add(hRead);
// }
fHRawList.Add(hRead);
hRead->Draw();
hRead->SetAxisRange(xaxismin, xaxismax);
}
}//else
fr->Close();
rawCanvas->Update();
return;
}
void AculCalibration::ShowSpectra(const char* filename, TCanvas* rawCanvas, Option_t *option, Int_t xaxismin, Int_t xaxismax, const Int_t subaddress)
{
//filename: input .root file with saved filled histograms to be showed
//rawCanvas: canvas on which one you will see the spectrum
//option: THStack options
//xaxismin: Minimum channel, which will be displayed
//xaxismax: Maximum channel, which will be displayed
//subaddress:
TString opt = option;
opt.ToLower();
if (!rawCanvas) {
Error("ShowRawSpectra", "You have to assign TCanvas for raw spectra drawing");
return;
}
rawCanvas->Clear();
TFile fr(filename);
if (fr.IsOpen() == 0) {
Error("ShowRawSpectra", "File %s was not opened and won't be processed", filename);
return;
}
TList *histList;
histList = fr.GetListOfKeys();
Int_t listEntries = histList->GetEntries();
TH1 *hDraw = 0;
DeleteStacks();
if (subaddress >= listEntries) {
fCurrentHStack = new THStack();
for (Int_t i = 0; i < listEntries; i++) { //zkontrolovat hranice
Info("ShowRawSpectra", "Histogram with spectrum of subaddress %d is loading", i);
fr.GetObject(histList->At(i)->GetName(), hDraw);
if (hDraw) {
hDraw->SetDirectory(0);
//hDraw->SetAxisRange(xaxismax, xaxismin); //nefunguje
fCurrentHistList.Add(hDraw);
fCurrentHStack->Add(hDraw);
}
}//for
if ( !fCurrentHStack->GetHists()->IsEmpty() ) {
Info("ShowRawSpectra", "Histogram stack drawing");
fCurrentHStack->Draw(opt.Data());
}
}//if all subaddresses
else {
//zkontrolovat
fr.GetObject(histList->At(subaddress)->GetName(), hDraw);
if (hDraw) {
hDraw->SetAxisRange(xaxismin, xaxismax, "X");
hDraw->Draw();
hDraw->SetDirectory(0);
fCurrentHistList.Add(hDraw);
}
}//else
fr.Close();
rawCanvas->Update();
return;
}
void AculCalibration::FillRawSpectraFile(const char* rawdatafile, const char* block, const char* treename, TCanvas* rawCanvas, Option_t *option, Int_t xaxismin, Int_t xaxismax)
{
//filename: input .root file containing spectra to be showed
//block:
//rawCanvas:
//xaxismin:
//xaxismax:
//variables to be became function parameter
TString opt(option);
opt.ToLower();
if (!rawCanvas) {
Error("ShowRawSpectra", "You have to assign TCanvas for raw spectra drawing");
return;
}
rawCanvas->Clear();
TFile fr(rawdatafile);
if (fr.IsOpen() == 0) {
Error("ShowRawSpectra", "File %s was not opened and won't be processed", rawdatafile);
return;
}
TTree *tr = (TTree*)fr.Get(treename);
char outputfile[300];
sprintf(outputfile, "%s[]Raw.root", block);
TFile fw(outputfile, opt.Data());
if (fw.IsOpen() == 0) {
Error("CalculateCalibParameters", "Output file %s was not created.", outputfile);
return;
}
if (fw.IsWritable() == 0) {
Error("CalculateCalibParameters", "Output file %s is not writable. Set option to \"RECREATE\".", outputfile);
return;
}
char address[40];
char histName[40];
char histTitle[40];
char fillCommand[40];
char fillCondition[40];
fw.cd();
TH1I *hRead = 0;
for (Int_t i = 0; i < 16; i++) { //zkontrolovat hranice
cout << i << endl; //predelat na info
hRead = new TH1I("name", "title", 4096, 0, 4095);
sprintf(address, "%s[%d]", block, i);
sprintf(histName, "%s[%d]", block, i);
sprintf(histTitle, "%s : %s", rawdatafile, histName);
hRead->SetName(histName);
hRead->SetTitle(histTitle);
sprintf(fillCommand, "%s >> %s", address, hRead->GetName());
sprintf(fillCondition, "%s > 0", address);
tr->Draw(fillCommand, fillCondition, "goff"); //prozkoumat goff
hRead->Write();
}//for
fw.Close();
fr.cd();
delete tr;
fr.Close();
rawCanvas->Update();
return;
}
Bool_t AculCalibration::CalculateCalibParameters(const char* inputfilename, const Int_t block, Int_t lowerchannel, Int_t upperchannel, Int_t lowersubaddress, Int_t uppersubaddress)
{
//function for calculation of calibrate parameters for DAQ system based on "Black Windows" program
//
// inputfile: root file with calibration spectra
// block: block to be calibrated as number
// lowerchannel: minimal channel from which the spectrum will be analysed
// upperchannel: maximal channel up to which the spectrum will be analysed
// lowersubaddress: block subaddress
// uppersubaddress: block subbaddress
if ( (block > BLOCKSNUMBER) || (block < 1) ) {
Error("CalculateCalibParameters", "Possible block values have to be in range 1 - %d", BLOCKSNUMBER - 1);
return kFALSE;
}
if ( (uppersubaddress - lowersubaddress) >= ADDRESSNUMBER ) {
Error("CalculateCalibParameters", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
return kFALSE;
}
//promenne potrebne pro fitovani
TF1 *calFunction = new TF1("calib", "pol1", 0, 1000); //predelat jako lokalni promennou fce
TGraph *calGraph = new TGraph(kRaNOPEAKS, fPeak, fEnergy); //lokalni promenna, dohodit pocet vstupu pomoci parametru
//auxiliary variables, particularly for text parameter fields
Char_t outputfilename[40];
Char_t address[40];
Char_t histName[40];
Char_t fillCommand[40];
Char_t fillCondition[40];
Int_t fitControl = 0;
Char_t histTitle[40];
//creation of the output text file
if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "C3[%d][].cal", block); }
else {
if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "C3[%d][%d].cal", block, lowersubaddress); }
else { sprintf(outputfilename, "C3[%d][%d-%d].cal", block, lowersubaddress, uppersubaddress); }
}
// fFileName = outputfilename; //doubtful
// cout << "hovno1" << endl;
// cout << outputfilename << endl;
// fOutCalFile.open(outputfilename);
ofstream outcalfile;
outcalfile.open(outputfilename);
// cout << "hovno2" << endl;
// if (!fOutCalFile.is_open()) {
if (!outcalfile.is_open()) {
Error("CalculateCalibParameters", "Output file %s was not opened", outputfilename);
return kFALSE;
}
//creation of the output root file
if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "C3[%d][].root", block); }
else {
if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "C3[%d][%d].root", block, lowersubaddress); }
else { sprintf(outputfilename, "C3[%d][%d-%d].root", block, lowersubaddress, uppersubaddress); }
}
fCalInformation = new TFile(outputfilename, "RECREATE");
if (fCalInformation->IsOpen() == 0) {
Error("CalculateCalibParameters", "File %s was not opened and won't be processed", outputfilename);
return kFALSE;
}
//input file with raw data opening
TFile *fr = new TFile(inputfilename);
if (fr->IsOpen() == 0) {
Error("CalculateCalibParameters", "File %s was not opened and won't be processed", inputfilename);
return kFALSE;
}
TTree *tr = (TTree*)fr->Get("RAW");
if (!tr) {
Error("CalculateCalibParameters", "Tree \"RAW\" was not found in file %s", inputfilename);
return kFALSE;
}
//raw data histogram filling
TH1I *hRaw = 0;
for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
Info("\n\nCalculateCalibParameters", "Calculating calibration parametres for address C3[%d][%d].", block, i);
// cout << endl
// << endl
// << "Calculating calibration parametres for address C3[" << block << "][" << i << "]." << endl;
//TH1I object preparing
hRaw = new TH1I("name", "title", 4096, 0, 4095);
sprintf(address, "C3[%d][%d]", block, i);
sprintf(histName, "H3[%d][%d]", block, i);
hRaw->SetName(histName);
sprintf(fillCommand, "%s >> %s", address, hRaw->GetName());
sprintf(fillCondition, "%s > 0", address);
//filling from the .root raw data file and content arrangement
tr->Draw(fillCommand, fillCondition, "goff");
if (lowerchannel != 0) {
for (Int_t i = 0; i < lowerchannel; i++) {
hRaw->SetBinContent(i, 0);
}
}
if (upperchannel != 0) {
for (Int_t i = upperchannel; i < 4095; i++) {
hRaw->SetBinContent(i, 0);
}
}
//spectrum analysis
fitControl = PeaksFitting(hRaw, "", fFitMinSigma);
// cout << "Value of fitControl is: " << fitControl << endl;
Info("CalculateCalibParameters", "Value of fitControl is: %d", fitControl);
//incorrectly treated spectrum output
if (fitControl != 0) {
// fOutCalFile << setw(39) << fitControl << endl;
outcalfile << setw(39) << fitControl << endl;
fCalInformation->cd();
hRaw->SetLineColor(2); //red
hRaw->Write();
continue;
}
//correctly treated spectrum saving
if (fCalInformation->IsOpen()) {
fCalInformation->cd();
hRaw->SetLineColor(3); //green
hRaw->Write();
}
//calibration parameters calculation
for (Int_t j = 0; j < kRaNOPEAKS; j++) { //delat podle poctu zkoumanych piku
calGraph->SetPoint(j, fPeak[j], fEnergy[j]); //calibration graph filling
// cout << "vypis:\t fPeak = " << fPeak[i] << "\tfEnergy = " << fEnergy[i] << endl;
}//for
calGraph->Fit(calFunction, "", "goff", 0, 4096); //omezit hlasitost fitovani, udelat volitelne
// fOutCalFile << std::right
outcalfile << std::right
<< setw(2) << "3"
<< setw(4) << block
<< setw(4) << i
<< setw(12) << setprecision(4) << calFunction->GetParameter(1)
<< setw(12) << setprecision(4) << calFunction->GetParameter(0)
<< setw(5) << fitControl
// << setw(10) << calFunction->GetParameter(1)*2.35*(fSigma[0] + fSigma[1] + fSigma[2] + fSigma[3])/4.
<< endl;
fA[block][i] = calFunction->GetParameter(1);
fB[block][i] = calFunction->GetParameter(0);
}//for
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//urcite by se to dalo hodit do samostatne fce
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//outputfile with calibrated spectra
if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "C3[%d][]E.root", block); }
else {
if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "C3[%d][%d]E.root", block, lowersubaddress); }
else { sprintf(outputfilename, "C3[%d][%d-%d]E.root", block, lowersubaddress, uppersubaddress); }
}
TFile *fw = new TFile(outputfilename, "RECREATE");
if (fw->IsOpen() == 0) {
Error("\nCalculateCalibParameters", "File %s was not created and won't be processed\n\n", outputfilename);
// cout << endl << "File " << outputfilename << " was not created and won't be processed" << endl << endl;
return kFALSE;
}
AculRaw *eventr = new AculRaw();
Long64_t noEntries = tr->GetEntries(/*getCondition*/);
tr->SetBranchAddress("channels", &eventr);
TH1F *hEnergy = 0;
TRandom3 ranGen(1);
for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
cout << "Calibration spectrum from address C3[" << block << "][" << i << "]." << endl;
sprintf(histName, "C3[%d][%d]", block, i);
sprintf(histTitle, "%s: %s", inputfilename, histName);
hEnergy = new TH1F(histName, histTitle, 10000, 0., 10.);
for (Int_t j = 0; j < noEntries; j++) {
tr->GetEntry(j);
if ((eventr->C3[block][i]) > 0) {
hEnergy->Fill( (eventr->C3[block][i] + ranGen.Uniform(-0.5, 0.5) )*fA[block][i] + fB[block][i] );
}
}
fw->cd();
hEnergy->Write();
}//for
fw->Close();
/////////////////////////////////////////////////////////
fr->Close();
fCalInformation->Close();
delete fCalInformation; //pokusne
// fOutCalFile.close();
outcalfile.close();
return kTRUE;
}
Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const char* block, const Int_t address, const char* treename, Int_t lowerchannel, Int_t upperchannel, Int_t lowersubaddress, Int_t uppersubaddress)
{
//function is not completely ready to use
//
//
//function for calculation of calibrate parameters for DAQ system based on "Go4"
//
// inputfile: root file with calibration spectra
// block: block name to be calibrated
// lowerchannel: minimal channel from which the spectrum will be analysed
// upperchannel: maximal channel up to which the spectrum will be analysed
// lowersubaddress: block subaddress
// uppersubaddress: block subbaddress
if (kRaNOPEAKS == 0) {
Error("CalculateCalibParameters", "Alpha source parameters was not red");
return 0;
}
if ( (address > BLOCKSNUMBER) || (address < 1) ) {
Error("CalculateCalibParameters", "Possible address values have to be in range 1 - %d", BLOCKSNUMBER - 1);
return kFALSE;
}
//muzu nechat
if ( (uppersubaddress - lowersubaddress) >= ADDRESSNUMBER ) {
Error("CalculateCalibParameters", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
return 0;
}
//probrat podle potreby
//auxiliary variables, particularly for text parameter fields
char outputfilename[300];
//creation of the output text file
if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "%s[].cal", block); }
else {
if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "%s[%d].cal", block, lowersubaddress); }
else { sprintf(outputfilename, "%s[%d-%d].cal", block, lowersubaddress, uppersubaddress); }
}//if
ofstream outcalfile;
outcalfile.open(outputfilename);
if (!outcalfile.is_open()) {
Error("CalculateCalibParameters", "Output file %s was not opened", outputfilename);
return 0;
}//if
//predelat podle nazvu bloku
//creation of the output root file
if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "%s[].root", block); }
else {
if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "%s[%d].root", block, lowersubaddress); }
else { sprintf(outputfilename, "%s[%d-%d].root", block, lowersubaddress, uppersubaddress); }
}
fCalInformation = new TFile(outputfilename, "RECREATE");
if ( !fCalInformation->IsOpen() ) {
Error("CalculateCalibParameters", "File %s was not opened and won't be processed", outputfilename);
return 0;
}
//nechat
//input file with raw data opening
TFile *fr = new TFile(inputfile);
if ( !fr->IsOpen() ) {
Error("CalculateCalibParameters", "File %s was not opened and won't be processed", inputfile);
return 0;
}
TTree *tr = (TTree*)fr->Get(treename);
if (!tr) {
Error("CalculateCalibParameters", "Tree %s was not found in file %s", treename, inputfile);
return 0;
}
//promenne potrebne pro fitovani: presunout nize
//pohlidat delete
TF1 *calFunction = new TF1("calib", "pol1", 0, 1000); //predelat jako lokalni promennou fce
TGraph *calGraph = new TGraph(kRaNOPEAKS, fPeak, fEnergy); //lokalni promenna, dohodit pocet vstupu pomoci parametru
Char_t detectorChannel[100];
Char_t histName[100];
Char_t fillCommand[512];
Char_t fillCondition[200];
Int_t fitControl = 0;
//predelat nazvy histogramu
//zrusit cyklus, napsat jako fci
//raw data histogram filling
TH1I *hRaw = 0;
for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
printf("\n\n");
Info("CalculateCalibParameters", "Calculating calibration parametres for detector channel %s[%d].", block, i);
//TH1I object preparing
hRaw = new TH1I("name", "title", 4096, 0, 4095); //nastavovat hranice histogramu podle parametru fce
sprintf(detectorChannel, "%s[%d]", block, i);
sprintf(histName, "Hist%s[%d]", block, i);
hRaw->SetName(histName);
sprintf(fillCommand, "%s >> %s", detectorChannel, hRaw->GetName());
sprintf(fillCondition, "%s > 0", detectorChannel);
//filling from the .root raw data file and content arrangement
tr->Draw(fillCommand, fillCondition, "goff");
if (lowerchannel != 0) { //zbytecna cast //
for (Int_t i = 0; i < lowerchannel; i++) { //
hRaw->SetBinContent(i, 0); //
} //muzu napsat jednoduseji:
} // vytvaret histogram s pozadovanym rozsahem
if (upperchannel != 0) { //zbytecna cast //
for (Int_t i = upperchannel; i < 4095; i++) { // taky ovsem nemusim, chtelo by to jeste zvazit
hRaw->SetBinContent(i, 0); //
} //
} //
//spectrum analysis
fitControl = PeaksFitting(hRaw, "Q", fFitMinSigma);
Info("CalculateCalibParameters", "Value of fitControl is: %d", fitControl); //ok
//incorrectly treated spectrum output
if (fitControl != 0 && fCalInformation->IsOpen()) { //ok
outcalfile << setw(39) << fitControl << endl;
fCalInformation->cd();
hRaw->SetLineColor(2); //red
hRaw->Write();
continue;
}//if
//correctly treated spectrum saving
if (fCalInformation->IsOpen()) { //ok
fCalInformation->cd();
hRaw->SetLineColor(3); //green
hRaw->Write();
}//if
//calibration parameters calculation //ok
for (Int_t j = 0; j < kRaNOPEAKS; j++) { //delat podle poctu zkoumanych piku
calGraph->SetPoint(j, fPeak[j], fEnergy[j]); //calibration graph filling
}//for
calGraph->Fit(calFunction, "Q", "goff", 0, 4096); //omezit hlasitost fitovani, udelat volitelne, dodelat volby rozsahu
outcalfile
<< block << "\t"
<< address << "\t"
<< i << "\t"
<< setprecision(4) << calFunction->GetParameter(1) << "\t"
<< setprecision(4) << calFunction->GetParameter(0) << "\t\t"
<< fitControl
<< endl;
fA[address][i] = calFunction->GetParameter(1);
fB[address][i] = calFunction->GetParameter(0);
}//for
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//zapis histogramu v MeV do souboru .root
//
//urcite by se to dalo hodit do samostatne fce
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Char_t histTitle[40];
//outputfile with calibrated spectra
if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "%s[]E.root", block); }
else {
if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "%s[%d]E.root", block, lowersubaddress); }
else { sprintf(outputfilename, "%s[%d-%d]E.root", block, lowersubaddress, uppersubaddress); }
}
TFile *fw = new TFile(outputfilename, "RECREATE");
if (fw->IsOpen() == 0) {
Error("CalculateCalibParameters", "File %s was not created and won't be processed\n\n", outputfilename);
return 1;
}
TH1F *hEnergy = 0;
TRandom3 ranGen(1);
for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
Info("CalculateCalibParameters", "Calibration spectrum from address %s[%d]", block, i);
sprintf(histName, "Hist%s[%d]E", block, i);
sprintf(histTitle, "%s: %s", inputfile, histName);
hEnergy = new TH1F(histName, histTitle, 10000, 0., 10.);
sprintf(detectorChannel, "%s[%d]", block, i);
hEnergy->SetName(histName);
sprintf(fillCommand, "%f*(%s+%f)+%f >> %s", fA[address][i], detectorChannel, ranGen.Uniform(-0.5, 0.5), fB[address][i], hEnergy->GetName());
sprintf(fillCondition, "%s > 0", detectorChannel);
//filling from the .root raw data file and content arrangement
tr->Draw(fillCommand, fillCondition, "goff");
fw->cd();
hEnergy->Write();
}//for
fw->Close();
fr->Close();
fCalInformation->Close();
delete fCalInformation; //pokusne
outcalfile.close();
return 1;
}
Int_t AculCalibration::CalibrateBlock(const Char_t* inputfilename, const Int_t block, const Char_t* outputfilename, Int_t lowersubaddress, Int_t uppersubaddress)
{
//probably some of the obsolete functions, maybe does not do anything important
//input data root file opening
TFile *fr = new TFile(inputfilename);
if (fr->IsOpen() == 0) {
cout << endl << "File " << inputfilename << " was not opened and won't be processed" << endl << endl;
return 1;
}
TTree *tr = (TTree*)fr->Get("RAW");
if (!tr) {
cout << "Tree \"RAW\" was not found in file " << inputfilename << "." << endl;
return 1;
}
//outputfile with calibrated spectra
TFile *fw = new TFile(outputfilename, "RECREATE");
if (fw->IsOpen() == 0) {
cout << endl << "File " << outputfilename << " was not created." << endl << endl;
return 1;
}
//check for the calibration parameters
for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
if ( fA[block][i] == 0 && fB[block][i] == 0) {
Warning("CalibrateBlock", "Calibration parameters are empty");
return 2;
}
}
//auxiliary variables for histogram names
Char_t histName[40];
Char_t histTitle[40];
//beginning of the tree reading
AculRaw *eventr = new AculRaw();
Long64_t noEntries = tr->GetEntries(/*getCondition*/);
tr->SetBranchAddress("channels", &eventr);
TH1F *hEnergy = 0;
TRandom3 ranGen(1);
for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
cout << "Calibration spectrum from address C3[" << block << "][" << i << "]." << endl;
sprintf(histName, "C3[%d][%d]", block, i);
sprintf(histTitle, "%s: %s", inputfilename, histName);
hEnergy = new TH1F(histName, histTitle, 10000, 0., 10.);
for (Int_t j = 0; j < noEntries; j++) {
tr->GetEntry(j);
if ((eventr->C3[block][i]) > 0) {
//v tomto miste se potom muze predelat jako bud vyplnovani histogramu nebo noveho stromu
hEnergy->Fill( (eventr->C3[block][i] + ranGen.Uniform(-0.5, 0.5) )*fA[block][i] + fB[block][i] );
}
}
fw->cd();
hEnergy->Write();
}//for
fw->Close();
return 0;
}
void AculCalibration::ShowAnalyzedSpectra(const char *filename, TCanvas* fittedRawCanvas, Int_t xaxismin, Int_t xaxismax, Int_t subaddress)
{
//This function displays analyzed spectrum from a file, divides the canvas into a sufficient number of pads and displays
//spectrums of each block subadress on the suitable pads or displays one selected spectrum .
//Selects the peaks in the histogram and displays on the histogram how the spectrum were fitted.
//
// filename: file .root containing analysed spectra
// fittedRawCanvas: canvas on which one you will see the spectrum
// xaxismin: Minimum channel, which will be displayed
// xaxismax: Maximum channel, which will be displayed
// subaddress:
if ( subaddress > ADDRESSNUMBER ) {
Error("ShowAnalyzedSpectra", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
return;
}
if (!fittedRawCanvas) {
Warning("ShowAnalyzedSpectra", "You have to assign TCanvas for fitted raw spectra drawing");
return;
}
TFile *fr = new TFile(filename, "READ");
if (!fr->IsOpen()) {
cout << "File " << filename << " was not opened" << endl;
return;
}
TList *histList;
histList = fr->GetListOfKeys();
Int_t listEntries = histList->GetEntries();
TH1I *hDraw = 0;
fittedRawCanvas->Clear();
// fittedRawCanvas->SetFillColor(10);
if ( (listEntries > 1) && (listEntries <= 8) ) {
fittedRawCanvas->Divide(2, 4);
fittedRawCanvas->SetFillColor(10);
}
if ( (listEntries > 8) && (listEntries <= 16) ) {
fittedRawCanvas->Divide(4, 4);
fittedRawCanvas->SetFillColor(10);
}
if (subaddress >= listEntries) {
for (Int_t i = 0; i < listEntries; i++) {
fittedRawCanvas->cd(i+1);
fr->GetObject(histList->At(i)->GetName(), hDraw);
if (hDraw) {
hDraw->SetAxisRange(xaxismin, xaxismax, "X");
hDraw->Draw();
hDraw->SetDirectory(0);
// if (fHAnalyzedList) {
// fHAnalyzedList->Add(hDraw);
// }
fHAnalyzedList.Add(hDraw);
}
}//for
}
else {
fr->GetObject(histList->At(subaddress)->GetName(), hDraw);
if (hDraw) {
hDraw->SetAxisRange(xaxismin, xaxismax, "X");
hDraw->Draw();
hDraw->SetDirectory(0);
fHAnalyzedList.Add(hDraw);
}
}
fr->Close();
fittedRawCanvas->Update();
return;
}
void AculCalibration::ShowEnergySpectra(const char *filename, TCanvas* energyCanvas, const Int_t subaddress, Option_t* option, Double_t xaxismin, Double_t xaxismax)
{
//Displays the spectrum of the selected subbaddress block in MeV
//
// filename: file .root containing calibrated spectra in MeV
// energyCanvas: : canvas on which one you will see the spectrum
// subaddress: block subaddress which will be drawn
// option: sum ,+ ,c
// xaxismin: Minimum channel, which will be displayed
// xaxismax: Maximum channel, which will be displayed
if ( subaddress > ADDRESSNUMBER ) {
Error("ShowEnergySpectra", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
return;
}
if (!energyCanvas) {
Warning("ShowEnergySpectra", "You have to assign TCanvas for fitted raw spectra drawing");
return;
}
TString opt = option;
opt.ToLower();
TFile *fr = new TFile(filename, "READ");
if (!fr->IsOpen()) {
cout << "File " << filename << " was not opened" << endl;
return;
}
TList *histList;
histList = fr->GetListOfKeys();
Int_t listEntries = histList->GetEntries();
TH1F *hDraw = 0;
energyCanvas->Clear();
energyCanvas->SetFillColor(10);
if ( (listEntries > 1) && (listEntries <= 8) ) {
energyCanvas->Divide(2, 4);
energyCanvas->SetFillColor(10);
}
if ( (listEntries > 8) && (listEntries <= 16) ) {
energyCanvas->Divide(4, 4);
energyCanvas->SetFillColor(10);
}
if (subaddress >= listEntries) {
if (opt.Contains("sum")) {
energyCanvas->cd(0);
for (Int_t i = 0; i < listEntries; i++) {
fr->GetObject(histList->At(i)->GetName(), hDraw);
if (hDraw) {
hDraw->SetDirectory(0);
//hDraw->SetAxisRange(xaxismin, xaxismax, "X");
if (opt.Contains("c")) { hDraw->SetLineColor(i+1); }
if (opt.Contains("c") && opt.Contains("+")) { hDraw->SetFillColor(i+1); }
// fHEnergyStack->Add(hDraw);
fHEnergyStack.Add(hDraw);
}
}
// if (fHEnergyStack) {
// if (opt.Contains("+")) { fHEnergyStack->Draw(); }
// else { fHEnergyStack->Draw("nostack"); }
// }
if (opt.Contains("+")) { fHEnergyStack.Draw(); }
else { fHEnergyStack.Draw("nostack"); }
}
else {
for (Int_t i = 0; i < listEntries; i++) {
energyCanvas->cd(i+1);
fr->GetObject(histList->At(i)->GetName(), hDraw);
if (hDraw) {
hDraw->SetAxisRange(xaxismin, xaxismax, "X");
hDraw->Draw();
hDraw->SetDirectory(0);
// if (fHEnergyList) {
// fHEnergyList->Add(hDraw);
// }
fHEnergyList.Add(hDraw);
}
}//for
}//else
}//if
else {
fr->GetObject(histList->At(subaddress)->GetName(), hDraw);
energyCanvas->cd(0);
if (hDraw) {
hDraw->SetAxisRange(xaxismin, xaxismax, "X");
hDraw->Draw();
hDraw->SetDirectory(0);
// if (fHEnergyList) {
// fHEnergyList->Add(hDraw);
// }
fHEnergyList.Add(hDraw);
}
}//else
fr->Close();
// fFileName = filename;
// fFileName.Resize(fFileName.Length() - 6);
// fFileName.Append(".cal", 4);
energyCanvas->Update();
return;
}
Bool_t AculCalibration::AddCalFileToList(const char* calfilelist)
{
//this function does not work at the moment
//some problem with TString object fFileName
TString fl = calfilelist;
fl.ToLower();
ofstream fw;
fw.open(fl.Data(), ofstream::app);
if (!fw.is_open()) {
cout << "File " << fl.Data() << " was not opened" << endl;
return kFALSE;
}
// fw << fFileName.Data() << endl;
fw.close();
return kTRUE;
}
void AculCalibration::ClearHistograms(Option_t* option)
{
//clear THStack and TObjArray members
//this function will be removed as soon as possible
TString opt = option;
opt.ToLower();
// fHRawList->Add(hRead);
fHRawList.Clear();
fHAnalyzedList.Clear();
fHEnergyList.Clear();
fHEnergyStack.Clear();
return;
}
void AculCalibration::MakeCalibrationFile(Char_t* calibrationfile, Char_t *calfilelist)
{
//calibrationfile: file with calibration parameters to be created
//calfilelist: file containing list of existing text files with calibration parameters
ifstream calListR;
calListR.open(calfilelist);
if( !calListR.is_open() ) {
cout << "File with list of calibration files was not opened" << endl;
return;
}
//asi fce Reset()
for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
fA[i][j] = 0;
fB[i][j] = 0;
}
}
const Int_t lineLength = 100;
char line[lineLength];
char filename[50];
// ifstream calFileR;
Int_t crate, i, j, id;
// int crate, i, j, id;
char cA[40], cB[40], cSigma[40];
while (!calListR.eof()) {
calListR.getline(line, lineLength);
// cout << line << endl;
sscanf(line, "%s", filename);
cout << filename << endl;
ifstream calFileR;
calFileR.open(filename);
if (calFileR.is_open()) {
cout << filename << " processing" << endl;
calFileR.seekg(0);
while (!calFileR.eof()) {
// cout << " in inner while" << endl;
calFileR.getline(line, lineLength);
sscanf(line, "%d %d %d %s %s %d %s", &crate, &i, &j, cA, cB, &id, cSigma);
// cout << line << endl;
if (id == 0) {
fA[i][j] = static_cast<Double_t>(atof(cA));
fB[i][j] = static_cast<Double_t>(atof(cB));
// fMeanSigma[i][j] = static_cast<Double_t>(atof(cSigma));
cout << fA[i][j] << "\t" << fB[i][j] << endl;
}//if
}//while
calFileR.close();
cout << "calFileR was closed" << endl << endl;
}//if
}//while
calListR.close();
ofstream CalibFileW;
CalibFileW.open(calibrationfile);
if (!CalibFileW.is_open()) {
cout << "Calibration file was not opened" << endl;
return;
}
for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
if (fA[i][j] != 0) {
CalibFileW << std::right
<< setw(2) << "3"
<< setw(4) << i
<< setw(4) << j
<< setw(12) << fA[i][j]
<< setw(12) << fB[i][j]
// << setw(10) << fMeanSigma[i][j]
<< endl;
}
}
}
CalibFileW.close();
return;
}
void AculCalibration::DeleteStacks(Option_t* option) {
if (fCurrentHStack) {
delete fCurrentHStack;
fCurrentHStack = NULL;
}
fCurrentHistList.Delete();
return;
}
void AculCalibration::Reset()
{
//reset calibration parameters fA, fB, fC, fD to zero
for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
fA[i][j] = 0;
fB[i][j] = 0;
fC[i][j] = 0;
fD[i][j] = 0;
// fMeanSigma[i][j] = 0;
}
}
return;
}
#pragma once
#include <TObject.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1I.h>
#include <TPolyMarker.h>
#include <TF1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TGraph.h>
#include <TObjArray.h>
#include <TRandom3.h>
#include <THStack.h>
#include <TString.h>
#include <TSpectrum.h>
#include "AculRaw.h" //potreba zaridit nezavislost na teto tride
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#define DEFAULTNOPEAKS 20
#define BLOCKSNUMBER 24
#define ADDRESSNUMBER 16
using std::cout;
using std::endl;
using std::setw;
using std::setprecision;
using std::stringstream;
using std::ostringstream;
class AculCalibration : public TObject
{
public:
//smysl jako verejne globalni promenne maji:
//fNOSpectra - pocet zkalibrovanych spekter
//???? - pocet spravne zkalibrovanych spekter
//???? - pocet nespravne zkalibrovanych spekter
//fEnergy[4] - tabulka s energiemi piku, nacita se zvenci
//private:
// TObjArray *fHRawList; //list of raw histograms, list is set to owner
TObjArray fHRawList; //list of raw histograms, list is set to owner
TObjArray fHAnalyzedList; //list of fitted and analyzed histograms, list is set to owner
TObjArray fHEnergyList; //list of calibrated histograms, list is set to owner
THStack fHEnergyStack; //some stack
THStack *fCurrentHStack;
//dodelat current histograms
TObjArray fCurrentHistList;
// TRandom3 *fRanGen;
// static TRandom3 fRanGen; //!
// THStack *fHStack;
//parameters to be read from file
Int_t kRaNOPEAKS;
Double_t fEnergy[DEFAULTNOPEAKS];
Double_t fLowerPeakRelativeHight; //pouziva se, private
Double_t fUpperPeakRelativeHight; //pouziva se, private, nastavit nenulovou prednastavenou hodnotu
Double_t fPeakPositionTolerance; //pouziva se, private
Width_t fFitFuncLineWidth; //private
Double_t fFitMinSigma; //pouziva se, private
Double_t fFitPeakThreshold; //pouziva se, private, prozkoumat, k cemu vlastne slouzi ve fci ShowPeaks, popremyslet o vhodnem prednastaveni v konstruktoru
//tyto promenne jsou smyslem tridy, mely by tedy byt snad jako jedine verejne
Double_t fA[BLOCKSNUMBER][ADDRESSNUMBER]; //kalibracni parametry, f(x) = fA*x + fB
Double_t fB[BLOCKSNUMBER][ADDRESSNUMBER]; //kalibracni parametry, f(x) = fA*x + fB
Double_t fC[BLOCKSNUMBER][ADDRESSNUMBER]; //treti kalibracni parametr, jine zavislosti nez pol1
Double_t fD[BLOCKSNUMBER][ADDRESSNUMBER]; //ctvrty kalibracni parametr
//smysl je velmi pochybny
TFile *fCalInformation;
// AculCalibratedData *fCalData;
//private:
Double_t fPeak[DEFAULTNOPEAKS]; //v teto promenne je ulozena momentalni hodnota piku v kanalech, zejmena v jedne fci, mozno udelat ji jako lokalni, bude navratovou hodnotou fce PeaksFitting, predelat delku pole
public:
AculCalibration(); //
AculCalibration(const char* parfile); //
virtual ~AculCalibration();
ClassDef(AculCalibration, 1);
//public
Bool_t SetInputParameters(const char* inputparfile = "parforcal.par");
Bool_t SetCalibrationParameters(const char* calparfile);
//public
void PrintInputParameters();
void PrintCalibrationParameters(const Int_t blockmin = 1, const Int_t blockmax = BLOCKSNUMBER - 1);
//public
Bool_t CalculateCalibParameters(const char* inputfilename, const Int_t block, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t lowersubaddress = 0, Int_t uppersubaddress = 15); //calculate calibration parameters for given block in given file
Bool_t CalculateCalibParameters(const char* inputfile, const char* block, const Int_t address, const char* treename, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t lowersubaddress = 0, Int_t uppersubaddress = 15); //calculate calibration parameters for given block in given file
Int_t CalibrateBlock(const Char_t* inputfilename, const Int_t block, const Char_t* outputfilename, Int_t lowersubaddress = 0, Int_t uppersubaddress = 15); //vysvetlit, co je to outputfile
Int_t PeaksFitting(TH1I* hSpectrum, Option_t* option = "", Double_t sigmamin = 2); //possible options: "V", "Q", ""
Int_t SearchPeaks(const TH1 *hin, Double_t sigma = 2, Option_t *option = "", Double_t threshold = 0.05, const Int_t searchedpeaks = 100);
void FillRawSpectraFile(const char* rawdatafile, const char* block, const char* treename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096);
void ShowRawSpectra(const char* filename, const Int_t block, TCanvas* rawCanvas = NULL, Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = 16);
void ShowSpectra(const char* filename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = 16);
void ShowAnalyzedSpectra(const char* filename, TCanvas* fittedRawCanvas = NULL, Int_t xaxismin = 0, Int_t xaxismax = 4096, Int_t subaddress = 16);
void ShowEnergySpectra(const char *filename, TCanvas* energyCanvas = NULL, const Int_t subaddress = 16, Option_t* option = "", Double_t xaxismin = 0., Double_t xaxismax = 10.); //option: "sum", "c", "+",
// void
//dodelat funkce TTree* Get...(...) pro raw, anal i E spectra
Bool_t AddCalFileToList(const char* calfilelist = "CalFileList.log");
void ClearHistograms(Option_t* option = "");
void DeleteStacks(Option_t* option = "");
void MakeCalibrationFile(Char_t* calibrationfile, Char_t* calfilelist);
void Reset();
//private:
};
#include "AculConvert.h"
ClassImp(AculConvert);
AculConvert::AculConvert()
{
fRawEvent = new AculRaw();
}
AculConvert::~AculConvert()
{
delete fRawEvent;
}
Int_t AculConvert::ConvertRun(const char* experimentalrun, const Int_t nofirstfolder, const Int_t nolastfolder/*, const Char_t* runaddress*/)
{
TString fullWayToRun = experimentalrun;
// fullWayToRun.ReplaceAll("/", "\\"); //for win
fOutFaultFile.open("convertFault.log");
Char_t nameOfFolder[260];
// Int_t i = 0;
for (Int_t i = nofirstfolder; i <= nolastfolder; i++) {
sprintf(nameOfFolder, "%s.%03d", fullWayToRun.Data(), i);
ConvertFolder(nameOfFolder);
}
fOutFaultFile.close();
return 0;
}
Int_t AculConvert::ConvertFolder(const char* rawdatafolder)
{
TString fullWayToFolder = rawdatafolder;
// fullWayToFolder.ReplaceAll("/", "\\"); //for win
//necessary names preparing
Char_t nameOfFile[260];
Char_t nameOfRun[40];
Char_t nameOfFileProtocol[260];
// strcpy(nameOfRun, fullWayToFolder(fullWayToFolder.Last('\\') + 1, fullWayToFolder.Length()).Data()); //for win
strcpy(nameOfRun, fullWayToFolder(fullWayToFolder.Last('/') + 1, fullWayToFolder.Length()).Data()); //for linux
strcpy(strstr(nameOfRun, "."), "\0");
// sprintf(nameOfFileProtocol, "%s\\%s.PRO", fullWayToFolder.Data(), nameOfRun); //for win
sprintf(nameOfFileProtocol, "%s/%s.PRO", fullWayToFolder.Data(), nameOfRun); //for linux
// printf("%s\\%s.PRO\n", fullWayToFolder.Data(), nameOfRun); //for windows
printf("%s/%s.PRO\n", fullWayToFolder.Data(), nameOfRun); //for linux
//file protocol reading
if (ReadFileProtocol(nameOfFileProtocol) != 0) {
return -1;
}
//root file creating
Char_t rootFileName[260];
// sprintf(rootFileName, "%s\.root", fullWayToFolder(fullWayToFolder.Last('\\') + 1, fullWayToFolder.Length()).Data()); //for win
// sprintf(rootFileName, "%s\.root", fullWayToFolder(fullWayToFolder.Last('/') + 1, fullWayToFolder.Length()).Data()); //for linux
sprintf(rootFileName, "%s.root", fullWayToFolder(fullWayToFolder.Last('/') + 1, fullWayToFolder.Length()).Data()); //for linux
cout << endl << endl
<< "===========================================================" << endl
<< "== Raw data converting to file " << rootFileName << endl
<< "===========================================================" << endl << endl;
TFile *fw = new TFile(rootFileName, "RECREATE");
TTree *tw = new TTree("RAW", "Transcription");
tw->Bronch("channels", "AculRaw", &fRawEvent);
//folder processing
Int_t i = 0;
do {
// sprintf(nameOfFile, "%s\\%s.%03d", rawdatafolder, nameOfRun, i); //for win
sprintf(nameOfFile, "%s/%s.%03d", rawdatafolder, nameOfRun, i); //for linux
i++;
cout << endl << "=================================================" << endl;
} while ( ConvertFile(nameOfFile, /*fullWayToFolder,*/ fw) == 0 );
fw->Close();
return 0;
}
Int_t AculConvert::ConvertFile(const char* rawdatafile, TFile* rootRawFile)
{
//file with ACCULINNA raw data opening
FILE *dataFile;
TString fullWayToFile = rawdatafile;
// fullWayToFile.ReplaceAll("/", "\\");
// cout << fullWayToFile.Data() << endl;
//file for error output opening
if ((dataFile = fopen(fullWayToFile.Data(), "rb")) == NULL) {
cout << endl << "File " << fullWayToFile.Data() << " was not found\n";
if ( fOutFaultFile.is_open() ) {
fOutFaultFile << "File " << fullWayToFile.Data() << " was not found\n";
}
return -1;
}
cout << endl << endl
<< "Data file " << fullWayToFile.Data() << " was opened" << endl;
Char_t filename[50];
// sprintf(filename, "%s", fullWayToFile(fullWayToFile.Last('\\') + 1, fullWayToFile.Length()).Data()); //for windows
sprintf(filename, "%s", fullWayToFile(fullWayToFile.Last('/') + 1, fullWayToFile.Length()).Data()); //for linux
//ROOT file for transcripted data creating
TTree *treeRaw;
if (rootRawFile == NULL) {
Char_t rootFileName[260];
// sprintf(rootFileName, "%s\.root", filename);
sprintf(rootFileName, "%s.root", filename);
rootRawFile = new TFile(rootFileName, "RECREATE");
treeRaw = new TTree("RAW", "Transcription to .root file");
treeRaw->Bronch("channels", "AculRaw", &fRawEvent);
//file protocol reading
TString protocolName = filename;
protocolName.Replace(protocolName.First('.'), 4, ".PRO");
if (ReadFileProtocol(protocolName.Data()) != 0) {
return -1;
}
// ReadFileProtocol( protocolName.Data() );
}
//TTree with class structure
// AculRaw *rawEvent = new AculRaw();incorrEventsNumber
else {
treeRaw = (TTree*)(rootRawFile->Get("RAW"));
}
// treeRaw->Bronch("channels", "AculRaw", &fRawEvent);
//file reading
const int BUFFSIZE = 8500000;
static unsigned char buffer[BUFFSIZE];
for (Int_t i = 0; i < BUFFSIZE; i++) {
buffer[i] = 0;
}
//cout << "hovno4" << endl;
fseek(dataFile, 0, SEEK_SET);
if ( fread(buffer, sizeof(buffer), 1, dataFile) != 0 ) {
cout << "!!!Data file is larger than expected!!!" << endl;
if ( fOutFaultFile.is_open() ) {
fOutFaultFile << "!!!Data file " << filename << " is larger than expected!!!" << endl;
}
return -2;
}
cout << endl
<< "------------------------" << endl
<< " Buffer analysis " << endl
<< "------------------------" << endl;
//////////////////////////////////////////////////////////////////////////////////////////////////////
// event structure: //
// n-times 2 bytes cells with parameters according to file protocol //
// 2 bytes cell with number of FLOAT parameters //
// pairs of 2 bytes cells with identifier (crate, address, subaddress) and value //
// 2 bytes cell with end of the event sign (first byte, F0) and trigger identifier (second byte) //
//
// structure of 2 bytes cell: |lower|higher|
//////////////////////////////////////////////////////////////////////////////////////////////////////
Int_t bufferPosition = 0;
Int_t corrEventsNumber = 0;
Int_t incorrEventsNumber = 0;
Int_t parametresNumber = 0;
Int_t identifier; //2 bytes (0FFF mask)
//1111(empty positions)111(crateMask)11111(blockPositionMask)1111(addresMask)
const Int_t addressMask = 15;
const Int_t blockPositionMask = 31;
const Int_t crateMask = 7;
Int_t address = 0;
Int_t blockPosition = 0;
Int_t crate = 0;
UInt_t value = 0;
while (bufferPosition < BUFFSIZE) {
parametresNumber = static_cast<Int_t>(buffer[bufferPosition + 2*PERMPARNUMBER]);
//verification of the correct EndOfEvent identifier position
while ( static_cast<Int_t>(buffer[bufferPosition + 2*PERMPARNUMBER + 2/*2 bytes with parametres Number*/ + 4*parametresNumber + 1/*cell with trigger ID*/]) != ENDOFEVENT && bufferPosition < BUFFSIZE) {
bufferPosition++; //search for the beginning of the correct event
if ( static_cast<Int_t>(buffer[bufferPosition + 2*PERMPARNUMBER + 2 + 4*parametresNumber + 1]) == ENDOFEVENT ) {
incorrEventsNumber++; //increased after end of the event identifier
cout << "!!!damaged event!!!" << endl;
if ( fOutFaultFile.is_open() ) {
fOutFaultFile << "!!!damaged event!!!" << endl;
}
}
}//while
//////////////////////////////////
//// correct event processing ////
//////////////////////////////////
if ( static_cast<Int_t>(buffer[bufferPosition + 2*PERMPARNUMBER + 2/*2 bytes with parametres Number*/ + 4*parametresNumber/*vysvetlit*/ + 1/*cell with trigger ID*/]) == ENDOFEVENT ) {
//rawEvent reseting
fRawEvent->Reset();
//trigger
fRawEvent->trigger = static_cast<Int_t>(buffer[bufferPosition + 2*PERMPARNUMBER + 2/*2 bytes with parametres Number*/ + 4*parametresNumber/*vysvetlit*/ /*+ 1*/ /*cell with trigger ID*/]);
//dodelat MWPC`s
//floating parametres
for (Int_t i = 0; i < parametresNumber; i++) {
identifier = static_cast<Int_t>( buffer[bufferPosition + 2*PERMPARNUMBER + 2 + 4*i] + (buffer[bufferPosition + 2*PERMPARNUMBER + 3 + 4*i] << 8) );
address = 0;
address = identifier & addressMask;
blockPosition = 0;
blockPosition = (identifier >> 4) & blockPositionMask;
crate = 0;
crate = (identifier >> 9) & crateMask;
value = 0;
value = static_cast<UInt_t>( buffer[bufferPosition + 2*PERMPARNUMBER + 4 + 4*i] + (buffer[bufferPosition + 2*PERMPARNUMBER + 5 + 4*i] << 8) );
if (crate == 3) {
fRawEvent->C3[blockPosition][address] = static_cast<Int_t>(value&fMaskC3[blockPosition]);
}
}
bufferPosition = bufferPosition + 2*PERMPARNUMBER + 1 + 4*parametresNumber + 2 + 1;
corrEventsNumber++;
treeRaw->Fill();
}//if
}//while
cout << endl
<< "In file " << filename << " was found " << corrEventsNumber << " correct events" << endl
<< "In file " << filename << " was found " << incorrEventsNumber << " incorrect events" << endl;
treeRaw->AutoSave();
if (rootRawFile->GetName() == "ConvertFileOutput.root") { //promyslet, vhodny nazev
rootRawFile->Close();
}
//close data file
if (fclose(dataFile) == EOF) { //parameter file closing
cout << "File " << filename << " closing error\n";
if ( fOutFaultFile.is_open() ) {
fOutFaultFile << "File " << filename << " closing error\n";
}
return -1;
}
cout << endl
<< "Data file " << filename << " was closed" << endl << endl;
return 0;
}
Int_t AculConvert::ReadFileProtocol(const char* protocolname)
{
for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
fMaskC3[i] = 0;
}
//file protocol opening
ifstream protocolFile;
protocolFile.open(protocolname);
if (!protocolFile.is_open()) {
cout << "File " << protocolname << " was not found\n";
if ( fOutFaultFile.is_open() ) {
fOutFaultFile << "File " << protocolname << " was not found\n";
}
return -1;
}
const Int_t lineLength = 200;
Char_t line[lineLength];
UInt_t crate, address, subaddress;
unsigned int mask;
//reading parametres
while (!protocolFile.eof()) {
protocolFile.getline(line, lineLength);
if ( line[0] != '#' ) {
sscanf(line, "%*d %*s %*d %d %d %d %*c %x", &crate, &address, &subaddress, &mask);
if (crate == 3) {
fMaskC3[address] = static_cast<UInt_t>(mask); //writing parameters
}
}
}
protocolFile.close();
return 0;
}
////////////////////////////////////////////////////////////////////////
//// ////
//// Class with functions converting the ACCULINNA raw format ////
//// to ROOT raw format (class AculRaw). Functions converting ////
//// particular file, particular folder and series of folders ////
//// are available. ////
//// ////
////////////////////////////////////////////////////////////////////////
#pragma once
//#include "DllExport.h"
#include "AculRaw.h"
#include <TObject.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#define ENDOFEVENT 0xf0 //identifier of the end of one event
#define PERMPARNUMBER 10 //number of the permanent parameters in one event according to file protocol
using std::cout;
using std::endl;
using std::setw;
class /*DllExport*/ AculConvert
{
private:
AculRaw *fRawEvent;
ofstream fOutFaultFile; //file for error output
UInt_t fMaskC3[BLOCKSNUMBER];
//functions
Int_t ReadFileProtocol(const char* protocolname); //read file protocol and save mask values
public:
AculConvert();
virtual ~AculConvert();
ClassDef(AculConvert, 1);
//functions
Int_t ConvertRun(const char* runname, const Int_t nofirstfolder, const Int_t nolastfolder/*, const Char_t* runaddress = ""*/); //convert more folders
Int_t ConvertFolder(const char* foldername); //convert whole folder
Int_t ConvertFile(const char* datafile, TFile* rootRawFile = NULL); //convert one file, zmenit vnitrni vypisy
};
################################################################################
# AculData input with some variables
################################################################################
ACULDATALIBS := -lCore -lCint -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMathCore -lGpad -lGraf
# Add inputs and outputs from these tool invocations to the build variables
ACULDATA_HEADERS += \
$(ACULDATA)/AculCalibration.h \
$(ACULDATA)/AculConvert.h \
$(ACULDATA)/AculRaw.h \
$(ACULDATA)/ConfigDictionary.h \
$(ACULDATA)/linkdef.h
ACULDATACPP_SRCS += \
$(ACULDATA)/AculCalibration.cpp \
$(ACULDATA)/AculConvert.cpp \
$(ACULDATA)/AculRaw.cpp \
$(ACULDATA)/ConfigDictionary.cpp \
$(ACULDATA)/AculDataCint.cpp
ACULDATAOBJS += \
$(ACULDATA)/AculCalibration.o \
$(ACULDATA)/AculConvert.o \
$(ACULDATA)/AculRaw.o \
$(ACULDATA)/ConfigDictionary.o \
$(ACULDATA)/AculDataCint.o
ACULDATACPP_DEPS += \
$(ACULDATA)/AculCalibration.d \
$(ACULDATA)/AculConvert.d \
$(ACULDATA)/AculRaw.d \
$(ACULDATA)/ConfigDictionary.d
\ No newline at end of file
#include "AculRaw.h"
ClassImp(AculRaw);
AculRaw::AculRaw() // it is possible to replace whole text in the constructor by use of the Reset() function
{
trigger = 0;
for (Int_t i = 0; i < 4; i++) {
mwpcReg[i] = 0;
}
for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
C3[i][j] = 0;
}
}
};
AculRaw::~AculRaw()
{
}
Int_t AculRaw::Reset()
{
trigger = 0;
for (Int_t i = 0; i < 4; i++) {
mwpcReg[i] = 0;
}
for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
C3[i][j] = 0;
}
}
return 0;
}
////////////////////////////////////////////////////////////////
//// ////
//// Class with structure of the TTree for immediate ////
//// translation of the ACCULINNA raw data format to ////
//// ROOT format. This class correspond to the DAQ ////
//// structure. ////
//// ////
////////////////////////////////////////////////////////////////
#pragma once
//#include "DllExport.h"
#include <TObject.h>
#define BLOCKSNUMBER 24
#define ADDRESSNUMBER 16
class AculRaw : public TObject
{
public:
Int_t mwpcReg[4]; //information about MWPCx plane in the register format
Int_t C3[BLOCKSNUMBER][ADDRESSNUMBER]; //information from crate C3
Int_t trigger;
AculRaw();
virtual ~AculRaw();
ClassDef(AculRaw, 1);
//functions
Int_t Reset(); //reseting of the read values
};
#include "ConfigDictionary.h"
ClassImp(ConfigDictionary);
//////////////////////////////////////////////////////////////////////////////
// BEGIN_HTML
// <p><font size="4"><b>Config(uration)Dictionary class</b></font></p>
// <br>
// <i>Author: Bartlomiej Hnatio 2012-08-06</i>
// <br><br>
// This is very useful class to convert strings containing pairs "key"="value"
// into fast dictionaries. This strings are often read from external files.
// From dictionary created in such way one can easily
// extract values in few supported formats (integer, double, boolean, string)
// It can also be used in opposite direction: when dictionary is created with
// given values and keys, one can create config string, which is most
// convenient to write in external files.
// <br><br><br>
// <u>Most simple example of usage:</u>
// <br><br>Suppose you have two variables fD and fI:
// <pre>-------------------
// Double_t fD;
// Int_t fI;
//
//fD = 3.14;
//fI = 2012;
//-----------------------</pre>
// To save its parameters into file you can create ConfigDictionary class
// instance and use <b>SetDouble</b> and <b>SetInt</b> functions to
// insert parameters values with arbitrarly defined keys (let them be
// really long and descriptive in this example):
// <pre>---------------------
//ConfigDictionary CD;
//CD.SetDouble("Most important variable D",fD);
//CD.SetInt("Equally important integer I",fI);
//---------------------</pre>
// Now configuration string saved in <b>ConfigDictionary</b> CD
// can be obtained with <b>ToString</b> method and should look like:
//<pre>
//"Most important variable D"="3.14" "Equally important integer I"="2012"
//</pre>
// <br> It can be easily saved to a file using simplest <i>fstream</i>
// methods. And the advantage is that as a key one can use any string,
// so configuration file created in such way is very self-explanatory.
// <br>
// <br>
// Now lets suppose the opposite action - loading config from file.
// Imagine, that you have 1000 objects of A class, which config was saved
// to file - one object per line:
// <pre>-----------------
//"Most important variable D"="3.14" "Equally important integer I"="2012"
//"Equally important integer I"="1011" "Most important variable D"="8.15"
//"Most important variable D"="13.16" "Equally important integer I"="10"
//(...)
//----------------</pre>
// Please notice that order in which pairs of keys and values are placed
// in file doesn't make any difference to the ConfigDictionary class.
// To recreate objects, you just read each line containing single
// config string and then create with it ConfigDictionary object
// using special constructor with string as argument,
// or using <b>FromString</b> method. After this, you can
// extract parameters (using <b>GetDouble</b> and <b>GetInt</b> methods)
// with using same keys as were used for saving, and ConfigDictionary
// will return their values:
// <pre>----------------------
//string line;//This line you read from file
//ConfigDictionary CD(line);
//fD = CD.GetDouble("Most important variable D");
//fI = CD.GetInt("Equally important integer I");
//--------------------</pre>
//<br>
// And last but not least: what happens, if key requested by user doesn't
// exist in dictionary? (This can be caused by many reasons, mostly errors on
// user side, but not always).
// When you try to extract non-existent key using <b>Get</b> functions,
// an exception is risen. In normal program it should end execution and print
// some information about where program stopped working.
// But lets say that you don't want that, i.e. program may use default
// configs instead of those from files, or not all keys were that important.
// You can surround all uses of <b>Get</b> methods with try/catch clause.
// So when exception is risen, you will catch it, and decide, if you want
// the program to stop running, or anything else. Simple example is
// shown below:
//<pre>-------------------
//ConfigDictionary CD(some_string);
//try{//try to read important variables:
// double d = CD.GetDouble("crucial var");
// bool b = CD.GetBool("most important b");
// int i = CD.GetInt("unique ID");
// //...and so on
//}catch(std::string & e){//catch any kind of exception
// Error("Some crucial variable wasn't read, ending program!");
// return SOME_REALLY_BAD_ERROR_CODE;
//}
//try{//now less important, or optional:
// string name = CD.GetString("least important variable");
// double p = CD.GetDouble("optional parameter");
// //...and so on
//}catch(std::string & f){
// Info("Some optional variables wasn't read!");
//}
//--------------------</pre>
// END_HTML
//////////////////////////////////////////////////////////////////////////////
//_____________________________________________________________________________
ConfigDictionary::ConfigDictionary(){
//just empty map...
}
//_____________________________________________________________________________
ConfigDictionary::ConfigDictionary(std::string params){
//Just creates dictionary using FromString method
FromString(params);
}
//_____________________________________________________________________________
std::string ConfigDictionary::ToString(){
//Builds string that can be easily saved to file.
//Same format that can be read from files or from QtGui
//This should work same way for all uses
std::map<std::string,std::string>::iterator it;//Iterator to map elements
std::stringstream ss;//stream helpful with adding strings
//iterate whole dictionary:
for (it=configMap.begin();it!=configMap.end();++it){
//insert pairs "key"="value":
ss<<"\""<<it->first<<"\""<<"="<<"\""<<it->second<<"\""<<" ";
}
return ss.str();
}
//_____________________________________________________________________________
void ConfigDictionary::FromString(std::string params){
//params - TString containing list of key=value pairs
//
//Changes string formatted:
//"key1"="value1" "key2"="value with spaces 2" ...
//into map with keys and values
//Useful in lots of I/O functions,
//when we want to have a nice, readable format of data
std::stringstream loading(params);
std::string k,v;
while(!loading.fail()){
getline(loading,k,'\"');//removes everything to first "
getline(loading,k,'\"');//All chars between first two "" are the key
getline(loading,v,'\"');//removes all until third "
getline(loading,v,'\"');//All between another pair of "" is the value
if (!loading.fail())
configMap[k]=v;
}
}
//_____________________________________________________________________________
std::string ConfigDictionary::GetString(std::string key)throw(std::string){
//Extracts string from given key
//(if it exist, otherwise raises exception)
if (configMap.find(key) == configMap.end()){
Error("ConfigDictionary::GetString",
"Couldn't find the key: %s!",key.c_str());
throw(key);
}
return configMap[key];
}
//_____________________________________________________________________________
int ConfigDictionary::GetInt(std::string key)throw(std::string){
//Extracts integer from given key
//(if it exist, otherwise raises exception)
if (configMap.find(key) == configMap.end()){
Error("ConfigDictionary::GetInt",
"Couldn't find the key: %s!",key.c_str());
throw(key);
}
int returned=0;
//Convert string to int:
std::stringstream ss(configMap[key]);
ss>>returned;
return returned;
}
//_____________________________________________________________________________
double ConfigDictionary::GetDouble(std::string key)throw(std::string){
//Extracts integer from given key
//(if it exist, otherwise raises exception)
if (configMap.find(key) == configMap.end()){
Error("ConfigDictionary::GetDouble",
"Couldn't find the key: %s!",key.c_str());
throw(key);
}
double returned=0.0;
//Convert string to double:
std::stringstream ss(configMap[key]);
ss>>returned;
return returned;
}
//_____________________________________________________________________________
bool ConfigDictionary::GetBool(std::string key)throw(std::string){
//Extracts boolean from given key
//(if it exist, otherwise raises exception)
if (configMap.find(key) == configMap.end()){
Error("ConfigDictionary::GetBool",
"Couldn't find the key: %s!",key.c_str());
throw(key);
}
//Convert string to bool:
if (configMap[key].compare("true") == 0)
return true;
else return false;
}
//_____________________________________________________________________________
void ConfigDictionary::SetString(std::string key,std::string value){
//Sets value to key, no comments needed here...
configMap[key] = value;
}
//_____________________________________________________________________________
void ConfigDictionary::SetDouble(std::string key,double value){
//Sets value to key, converts double to string first
std::stringstream ss;
ss<<value;
configMap[key] = ss.str();
}
//_____________________________________________________________________________
void ConfigDictionary::SetInt(std::string key,int value){
//Sets value to key, converts int to string first
std::stringstream ss;
ss<<value;
configMap[key] = ss.str();
}
//_____________________________________________________________________________
void ConfigDictionary::SetBool(std::string key,bool value){
//Sets value to key, converts bool to string first
if (value == true) configMap[key] = "true";
else configMap[key] = "false";
}
#ifndef CONFIGDICTIONARY_H
#define CONFIGDICTIONARY_H
#include "TObject.h"
#include "ReturnCodes.h"
#include "TStopwatch.h"
#include "TROOT.h"
#include <string>
#include <map>
#include <sstream>
#include "TLorentzVector.h"
#include <iostream>
#include "TMath.h"
#include "TString.h"
#include "TTree.h"
class ConfigDictionary{
public:
typedef std::map<std::string,std::string>::iterator CDIter;
ConfigDictionary();
ConfigDictionary(std::string);
virtual ~ConfigDictionary(){};//empty virtual destructor
ClassDef(ConfigDictionary,1);
std::string ToString();
void FromString(std::string);
//These throw errors if couldn't find key:
std::string GetString(std::string)throw(std::string);
int GetInt(std::string)throw(std::string);
double GetDouble(std::string)throw(std::string);
bool GetBool(std::string)throw(std::string);
//These will always set 'something' into map:
void SetString(std::string,std::string);
void SetDouble(std::string,double);
void SetInt(std::string,int);
void SetBool(std::string,bool);
CDIter Begin(){return configMap.begin();};
CDIter End(){return configMap.end();};
private:
std::map<std::string,std::string> configMap;
};
#endif
#ifndef RETURN_CODES_H
#define RETURN_CODES_H
//Some useful return codes
const static int SUCCESS = 0;
const static int IOEXCEPTION = -1;
const static int NOTFOUND = -2;
const static int NULLPOINTER = -3;
const static int UNKNOWN = -4;
const static int FAILURE = -5;
const static int CDEXCEPTION = -6;
const static int EMPTYCONTAINER = -7;
const static int NOTDEFINED = -8;
const static int EXCEPTION = -11;
#endif
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class AculRaw;
#pragma link C++ class AculConvert;
#pragma link C++ class AculCalibration;
#pragma link C++ class ConfigDictionary;
#endif
################################################################################
# Be input with some variables
################################################################################
BELIBS := -lCore -lCint -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMathCore -lGpad -lGraf -lPhysics -lGeom #-lAculData -lDetectors -lTELoss
# Add inputs and outputs from these tool invocations to the build variables
BE_HEADERS += \
$(BE)/BeEvent.h \
$(BE)/BePureEvent.h \
$(BE)/BeReaction.h \
$(BE)/BeWork.h \
$(BE)/BinaryReaction.h \
$(BE)/linkdef.h
BECPP_SRCS += \
$(BE)/BeCint.cpp \
$(BE)/BeEvent.cpp \
$(BE)/BePureEvent.cpp \
$(BE)/BeReaction.cpp \
$(BE)/BeWork.cpp \
$(BE)/BinaryReaction.cpp
BEOBJS += \
$(BE)/BeCint.o \
$(BE)/BeEvent.o \
$(BE)/BePureEvent.o \
$(BE)/BeReaction.o \
$(BE)/BeWork.o \
$(BE)/BinaryReaction.o
BECPP_DEPS += \
$(BE)/BeCint.d \
$(BE)/BeEvent.d \
$(BE)/BePureEvent.d \
$(BE)/BeReaction.d \
$(BE)/BeWork.d \
$(BE)/BinaryReaction.d
This source diff could not be displayed because it is too large. You can view the blob instead.
/*
* BeEvent.h
*
* Created on: 9.12.2009
* Author: Vratislav
*/
#pragma once
#include <TObject.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TRandom3.h>
#include <TMath.h>
#include <TVector3.h>
#include <TRotation.h>
#include <TLorentzVector.h>
#include <TLorentzRotation.h>
#include <TCutG.h>
#include "../AculData/AculRaw.h"
#include "../AculData/AculCalibration.h"
#include "../Detectors/AnnularDetector.h"
#include "../TELoss/TELoss.h"
#include "../AculData/ConfigDictionary.h"
#include<iostream>
using std::cout;
using std::endl;
//TMath
using TMath::Power;
using TMath::Sqrt;
#define NOESEC 16 //number of sectors carrying full dE-E information
#define NORINGS 32 //number of all rings in the first detector of each telescope
#define NOHITS 16 //maximal number of hits
#define NOTEL 2 //number of telescopes
//materials used for simulations
//Si
#define SI_RHO 2.321 //in g.cm^-3
#define SI_NELEMENTS 1 //number of elements in material
#define SI_1_A 28.086 //A of the first element
#define SI_1_Z 14. //Z of the first element
#define SI_1_W 1. //weight of the first element
//CsI
#define CSI_RHO 4.51 //in g.cm^-3
#define CSI_NELEMENTS 2 //number of elements in material
#define CSI_1_A 132.91 //A of the first element
#define CSI_1_Z 55. //Z of the first element
#define CSI_1_W 0.5 //weight of the first element
#define CSI_2_A 126.9 //A of the second element
#define CSI_2_Z 53. //Z of the second element
#define CSI_2_W 0.5 //weight of the second element
//target properties (19.8.2012)
#define T_NELEMENTS 1 //number of elements in material
#define T_1_A 1.008 //A of the first element
#define T_1_Z 1. //Z of the first element
#define T_1_W 1. //weight of the first element
#define T_NORMAL_RHO 0.00008895 //density in g/cm^3 at 273 K and 10^5 Pa
#define T_TEMPERATURE 30. //temperature in K
#define T_PRESSURE 4. //pressure in bars
#define T_THICKNESS 6000. //target thickness in mcm
#define T_WIN_THICK 6. //window thickness in mcm
#define T_WIN_RHO 8.
#define T_WIN_NELEMENTS 3 //number of elements in material
#define T_WIN_1_A 51.996 //A of the first element
#define T_WIN_1_Z 24. //Z of the first element
#define T_WIN_1_W .08 //weight of the first element
#define T_WIN_2_A 55.847 //A of the first element
#define T_WIN_2_Z 26. //Z of the first element
#define T_WIN_2_W .74 //weight of the first element
#define T_WIN_3_A 58.69 //A of the first element
#define T_WIN_3_Z 28. //Z of the first element
#define T_WIN_3_W .18 //weight of the first element
//particles
#define kPROTON 11 //ZA
#define kALPHA 24 //ZA
//masses
#define kPMASS 938.272 //in MeV/c^2
#define kNMASS 939.565 //in MeV/c^2
#define kAMASS 3728.400 //in MeV/c^2
#define kLiMASS 5603.050 //in MeV/c^2
#define kBeMASS 5607.338 //in MeV/c^2
class BeEvent {
public: //melo by byt private
Int_t fTrigger; ///also used for state identificator fTrigger == 1J for state J, i.e 10 for 0+, 12 for 0+, 11 for 1-
Int_t fn[NOESEC];
Double_t fCalCorr[NOESEC];
//detectors
AnnularDetector fT11;
AnnularDetector fT12;
AnnularDetector fT13;
AnnularDetector fT21;
AnnularDetector fT22;
AnnularDetector fT23;
//deltaE in first detector, 32 sectors converted to 16 sectors
Double_t fdE1[NOESEC];
Double_t fdE2[NOESEC];
//energy and multiplicity information from CsI detectors
Double_t fE13aSec[NOESEC];
Double_t fE13pSec[NOESEC];
Double_t fE23aSec[NOESEC];
Double_t fE23pSec[NOESEC];
Int_t fMult13a;
Int_t fMult23a;
Int_t fMult13p;
Int_t fMult23p;
//time information from Tx2
Double_t fTau[NOTEL][NOESEC]; //time in ns
//hits
Int_t fNOHits[NOTEL]; //general number of hits with dE-E, phi and theta information
Int_t fNOSHits[NOTEL];
Int_t fNORHits[NOTEL]; //hits in first detector of telescope detected only by rings
Int_t fNOaHits[NOTEL]; //alpha particle hits before choice by time
Int_t fNOpHits[NOTEL]; //proton hits before choice by time
Int_t fNOtHits[NOTEL]; //hits searched by time
Int_t fSN[NOTEL][NOHITS];
Int_t fTauSN[NOTEL][NOHITS]; //time hit sector number
Int_t fRN[NOTEL][NOHITS];
Double_t fE1Sec[NOTEL][NOHITS]; //hit energy in 1st layer (sectors) of particular telescope
Double_t fE1Ring[NOTEL][NOHITS]; //hit energy in 1st layer (rings) of particular telescope
Double_t fE2Sec[NOTEL][NOHITS]; //hit energy in 2nd layer of particular telescope
Double_t fE3Sec[NOTEL][NOHITS]; //hit energy in 3rd layer of particular telescope
Double_t fTheta[NOTEL][NOHITS]; //in rad for lab system
Double_t fPhi[NOTEL][NOHITS]; //in rad for lab system
Double_t fhTau[NOTEL][NOHITS];
Double_t fT[NOTEL][NOHITS]; //in MeV for lab system
Double_t fP[NOTEL][NOHITS]; //in MeV for ?? system
Double_t fdE1calc[NOTEL][NOHITS]; //calculated energy in the first detector, it is used for particle identification
// Double_t fE1[NOHITS]; //zatim bez pouziti
Int_t fID[NOTEL][NOHITS]; //particle identificator = ZA; proton==11, alpha==24
Int_t fPM[2*NOHITS]; //particle markers
// //6Be information
Int_t fNOA; //number of detected alphas related to triple coincidence Be decay
Int_t fNOP; //number of detected protons related to triple coincidence Be decay
//laboratory system
//6Be
TLorentzVector f6BeLab;
Double_t f6BeThetaLab; //Theta_lab in rad
Double_t f6BePhiLab; //Phi_lab in rad
Double_t f6BePcLab; //abs(pc) in MeV
Double_t f6BeIM; //invariant mass in MeV
//neutron
TLorentzVector fNLab;
Double_t fNThetaLab; //Theta_lab in rad
Double_t fNPhiLab; //Phi_lab in rad
Double_t fNPcLab; //invariant mass in MeV
//proton1
TLorentzVector fP1Lab;
Double_t fP1ThetaLab; //Theta_lab in rad
Double_t fP1PhiLab; //Phi_lab in rad
Double_t fP1PcLab; //invariant mass in MeV
//proton2
TLorentzVector fP2Lab;
Double_t fP2ThetaLab; //Theta_lab in rad
Double_t fP2PhiLab; //Phi_lab in rad
Double_t fP2PcLab; //invariant mass in MeV
//alpha
TLorentzVector fALab;
Double_t fAThetaLab; //Theta_lab in rad
Double_t fAPhiLab; //Phi_lab in rad
Double_t fAPcLab; //invariant mass in MeV
//protons relative kinetic energy
Double_t fPPTrel;
//CM system 6Li-p
//6Be
TLorentzVector f6BeCM1;
Double_t f6BeThetaCM1; //Theta_cm in rad
Double_t f6BePhiCM1; //Phi_lab in rad
Double_t f6BePcCM1; //abs(pc) in MeV
//CM system 6Be
//6Be
TLorentzVector f6BeCM;
Double_t f6BeThetaCM; //Theta_cm in rad
Double_t f6BePhiCM; //Phi_lab in rad
Double_t f6BePcCM; //abs(pc) in MeV
//alpha
TLorentzVector fACM;
Double_t fAThetaCM; //Theta_cm in rad
Double_t fAPhiCM; //Phi_lab in rad
Double_t fAPcCM; //abs(pc) in MeV
//proton1
TLorentzVector fP1CM;
Double_t fP1ThetaCM; //Theta_cm in rad
Double_t fP1PhiCM; //Phi_lab in rad
Double_t fP1PcCM; //abs(pc) in MeV
//proton2
TLorentzVector fP2CM;
Double_t fP2ThetaCM; //Theta_cm in rad
Double_t fP2PhiCM; //Phi_lab in rad
Double_t fP2PcCM; //abs(pc) in MeV
//neutron
TLorentzVector fNCM;
Double_t fNThetaCM; //Theta_cm in rad
Double_t fNPhiCM; //Phi_lab in rad
Double_t fNPcCM; //abs(pc) in MeV
//general information
Double_t fQLiP; //Q of the reaction 6Li+p-->6Be+n
//Jacobi coordinates //fixme check them
//"T" system
Double_t fTpp; //"T" system
Double_t fTapp; //"T" system
Double_t fCosThetaTk; //"T" system
//"Y" system
Double_t fTap; //"Y" system
Double_t fTpap; //"Y" system
Double_t fCosThetaYk; //"Y" system
//variables used in Identification functions
//theta angles
Double_t fTheta1; //!
Double_t fQ1; // //todo apparently no information saved
Double_t fT6Li2; //todo apparently no information saved
Double_t fTheta2; //!
//computed energies
Double_t fE11aSec[NOESEC]; //computed energy in T11, fE11cSec = f(E12, x)
Double_t fE11pSec[NOESEC]; //compSetQLiP()uted energy in T11, fE11cSec = f(E12, x)
Double_t fE21aSec[NOESEC]; //computed energy in T21, fE21cSec = f(E22, x)
Double_t fE21pSec[NOESEC]; //computed energy in T21, fE21cSec = f(E22, x)
Double_t fE12aSec[NOESEC]; //computed energy in T12, fE12cSec = f(E13, x)
Double_t fE12pSec[NOESEC]; //computed energy in T12, fE12cSec = f(E13, x)
Double_t fE22aSec[NOESEC]; //computed energy in T22, fE21cSec = f(E23, x)
Double_t fE22pSec[NOESEC]; //computed energy in T22, fE21cSec = f(E23, x)
//used in function DeltaEId. and DeltaECal. functions, used just for verification
//these variables need not be be written into file
//first telescope
// Double_t fX11; //! //draha v T11, pouze docasne
// Double_t fX12; //! //draha v T12, pouze docasne
Double_t fXd112; //! //draha v mrtve vrstve mezi T11 a T12
Double_t fXd123; //! //draha v mrtve vrstve mezi T12 a T13
Int_t fRN1; //!
Int_t fSN1[NOHITS]; //!
//second telescope
// Double_t fX21; //! //draha v T11, pouze docasne
// Double_t fX22; //! //draha v T12, pouze docasne
Double_t fXd212; //! //draha v mrtve vrstve mezi T11 a T12
Double_t fXd223; //! //draha v mrtve vrstve mezi T12 a T13
Int_t fRN2; //!
Int_t fSN2; //!
//used for verification of deltaE calibration method for first telescope:
Double_t fE11v_0s[NOESEC]; //!
Double_t fE11v_0r; //!
Double_t fE11vs[NOESEC]; //!
Double_t fE11vr; //!
Double_t fE12v_0[NOESEC]; //!
Double_t fE12v[NOESEC]; //!
//used for verification of deltaE calibration method for first telescope:
Double_t fE21v_0s[NOESEC]; //!
Double_t fE21v_0r; //!
Double_t fE21vs[NOESEC]; //!
Double_t fE21vr; //!
Double_t fE22v_0[NOESEC]; //!
Double_t fE22v[NOESEC]; //!
//used for deltaE calibration of CsI
//alhpas
Double_t fE12_0a[NOESEC]; //!
Double_t fE12a[NOESEC]; //!
Double_t fE13_0a[NOESEC]; //!
//protons
Double_t fE12_0p[NOESEC]; //!
Double_t fE12p[NOESEC]; //!
Double_t fE13_0p[NOESEC]; //!
//alhpas
Double_t fE22_0a[NOESEC]; //!
Double_t fE22a[NOESEC]; //!
Double_t fE23_0a[NOESEC]; //!
//protons
Double_t fE22_0p[NOESEC]; //!
Double_t fE22p[NOESEC]; //!
Double_t fE23_0p[NOESEC]; //!
static TRandom3 ranE; //!
static TRandom3 ranChoice; //!
TELoss mSiAlphaV; //! //silicon
TELoss mSiAlpha; //! //silicon
TELoss mSiP; //! //silicon
TELoss mCsIAlpha; //! //CsI
TELoss mCsIP; //! //CsI
TELoss mH_P; //! //protons in hydrogen
TELoss mH_Alpha; //! //alpha in hydrogen
TELoss mWin_P; //! //protons in target window
TELoss mWin_Alpha; //! //alpha in target window
//particle indentification cuts
TFile *fCuts; //!
TCutG *cutAlpha1; //! 1st telescope
TCutG *cutProton1; //!
TCutG *cutAlpha2; //! 2nd telescope
TCutG *cutProton2; //!
//calibration parameters
AculCalibration *calSi; //! calibration parameters for Si detectors
AculCalibration *calCsIa; //! alpha parameters for CsI detectors
AculCalibration *calCsIp; //! proton parameters for CsI detectors
TString sSiCal; //! Si calibration parameters
TString sCsIprotonsCal; //! CsI calibration parameters for protons
TString sCsIalphaCal; //! CsI calibration parameters for alphas
//parameters
TString fConfigFile; //!
//geometry used for data reconstruction
Double_t fT1ExpPos; //! in mm
Double_t fT2ExpPos; //! in mm
//beam characteristics used in DA
Double_t fTBeamRec; //! in AMeV used for reconstruction //where is it? in the centre? on the target window?
// Double_t fTBeamRes; //! in MeV, sigma for gaus
Double_t fBeamX; //! in cm, x position of beam at target
Double_t fBeamY; //! in cm, y position of beam at target
Double_t fBeamX_sigma; //! in cm; x position sigma of beam at target
Double_t fBeamY_sigma; //! in cm; x position sigma of beam at target
//detector thicknesses in mcm
Double_t fX11_FD; //!
Double_t fX11; //!
Double_t fX11_BD; //!
Double_t fX12_FD; //!
Double_t fX12; //!
Double_t fX12_BD; //!
//#define XD13F 14. //fixme original value probably in Si equivalent
Double_t fX13_FD; //!
Double_t fX21_FD; //!
Double_t fX21; //!
Double_t fX21_BD; //!
Double_t fX22_FD; //!
Double_t fX22; //!
Double_t fX22_BD; //!
//#define XD23F 14. //fixme original value probably in Si equivalent
Double_t fX23_FD; //!
Double_t fCorrT1_0_7_CsI_A; //! correction for calibration of alphas in T1 CsI detectors
Double_t fCorrT1_8_15_CsI_A; //! correction for calibration of alphas in T1 CsI detectors
Double_t fCorrT1_0_7_CsI_P; //! correction for calibration of protons in T1 CsI detectors
Double_t fCorrT1_8_15_CsI_P; //! correction for calibration of protons in T1 CsI detectors
//telescope 2
Double_t fCorrT2_0_7_CsI_A; //! correction for calibration of alphas in T2 CsI detectors //probably same as MS (17. 12. 2012)
Double_t fCorrT2_8_15_CsI_A; //! correction for calibration of alphas in T2 CsI detectors
Double_t fCorrT2_0_7_CsI_P; //! correction for calibration of protons in T2 CsI detectors
Double_t fCorrT2_8_15_CsI_P; //! correction for calibration of protons in T2 CsI detectors
//thresholds
Double_t fRSTolerance; //! in MeV, tolerance between energy in sectors and rings for one hit
Double_t fProtonThr; //! in MeV, upper threshold for protons
TString fCutFile; //! file containing cuts
TString fCalibFilePath; //! path to files with calibration parameters
public:
BeEvent();
BeEvent(const char* configfile/*, const char* cutfile,*/ /*const TString calibfilepath*/); //todo finish implementation of this constructor
//change beam kinetic energy constant to a parameter
//probably the data member of BeEvent and then
//constructor will be in form of BeEvent(Double_t tbeam)
// BeEvent(BeEvent &);
virtual ~BeEvent();
ClassDef(BeEvent, 1);
void ReadConfigFile();
void FillEvent(AculRaw *rawevent);
void FillEventOld(AculRaw *rawevent, AculCalibration *calSi, AculCalibration *calA, AculCalibration *calP);
//functions for work with particular detectors
//opravit cislovani ringu v prvnim detektoru
void SetChannels(AculRaw *rawevent); //fill all variables with information in channels
void SetSiDetEnergies(AculRaw *rawevent); //fill all variables with energy information for T11, T12, T21, T22
void SetSiDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal); //fill all variables with energy information for T11, T12, T21, T22
void SetSiDetTimes(AculRaw *rawevent);
void SetSiDetTimesOld(AculRaw *rawevent, AculCalibration *cal);
void SetCsIalphaDetEnergies(AculRaw *rawevent);
void SetCsIalphaDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal);
void SetCsIprotonDetEnergies(AculRaw *rawevent); //zatim nepouzivat
void SetCsIprotonDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal); //zatim nepouzivat
void SetDeltaE();
void SetTheta1(Int_t ring);
void SetTheta2(Int_t ring);
//functions for work with hits
void SetSecHitEnergiesSi(Int_t tel); //for first telescope notel==0, for second notel==1
void SetSecHitEnergiesSiOld(Int_t tel); //for first telescope notel==0, for second notel==1
void SetRingHitEnergies(Int_t tel); //for first telescope notel==0, for second notel==1
void SetHitsEnergies(Int_t tel);
void SetEnergyHits(Int_t tel);
void SetThetas(Int_t tel);
void SetThetasOld(Int_t tel);
void SetThetasNew(Int_t tel); //fixme redo this function, it should take into account beam position
void SetPhis(Int_t tel);
void BeParticleIdentificationNew(Int_t tel);
void BeParticleIdentification(Int_t tel);
void ParticleMarkers();
void SetTimeHits(Int_t tel);
void SetTimeHitsNew(Int_t tel);
void SetHitsOld();
void SetHits();
void SetT(Int_t tel); //particle kinetic energy
void SetTNew(Int_t tel); //particle kinetic energy
Double_t CalcKinE(const Int_t telescope, const Int_t hit,
const Double_t cosTheta, const Double_t fDeadCsI,
const Double_t bDeadSi2, const Double_t fDeadSi2,
const Double_t bDeadSi1, const Double_t fDeadSi1);
void SetP(Int_t tel); //particle size of impuls
//general functions
void InitDetectors();
Bool_t WriteCondition();
void Reset();
Bool_t WriteConditionPure();
//physics
void SetProductsLab(); //set kinematical information about 6Be and n in laboratory system
void SetProductsLabNew(); //set kinematical information about 6Be and n in laboratory system
void SetQLiP();
void SetProductsCM1(); //set kinematical information about 6Be and n in CM system 6Li-p
void SetProductsCM(); //set kinematical information about 6Be and n in CM system of 6Be
//functions used for calibration
void Set6LiQ();
// void DeltaEIdentification1old(); //todo this function can be erased
// void DeltaEIdentification1();
// void DeltaEIdentification2();
void DeltaECalibration1(); //first telescope //todo probably can be erased
void DeltaECalibration2(); //second telescope //todo probably can be erased
private:
Double_t CsIEnergy(Int_t _address, Int_t _subaddress, AculRaw *_rawevent, AculCalibration *_cal, Double_t _x0, Double_t _ethr);
Double_t CsIcorrectedEnergy(Int_t _address, Int_t _subaddress, AculRaw *_rawevent, AculCalibration *_cal, Double_t _x0, Double_t _ethr);
Double_t CompareLosses(Double_t de_ex, Double_t e, Double_t x, Int_t particleID); //de_ex experimental dE, e experimental E in second layer, x range in the dE layer, ID==24 for alpha, ID==11 for proton
void HitIdentification(const Int_t tel, const Int_t hit, const Double_t eCsI, const Int_t particle, const Double_t amin, const Double_t amax, const Double_t bmin, const Double_t bmax, const Double_t cmin, const Double_t cmax);
Int_t NOParticlesTauE(Int_t tel, Int_t particleID); //number of particular particles which are detected as hit in proper time window
Int_t GetNOSecPhi(Int_t tel, Int_t esec); //get number of physical sector in first telescope of each detector
//CsI calibrations
void SetCalCorrs();
public:
//relativistic kinematics
static Double_t PC2(Double_t T, Double_t mc2); //neni potrebna, je obsazena v TLorentzVector
static Double_t T(Double_t pc2, Double_t mc2); //v podstate taktez
static Double_t LawOfCosines(Double_t pb, Double_t pc, Double_t alpha);
static Double_t ReactionQ(Double_t Ta1, Double_t Ta2, Double_t mac2, Double_t mbc2, Double_t theta);
static Double_t ReactionQ(Double_t Ea, Double_t Eb, Double_t Ec, Double_t Ed);
//auxiliary
static Double_t ReducedMass(Double_t m1, Double_t m2);
private:
void ReadCuts(/*const char* cutfile*/);
void CreateSiELosses();
void CreateCsIELosses();
void CreateTargetELosses();
};
/*
* BePureEvent.cpp
*
* Created on: 6.6.2011
* Author: vratik
*/
#include "BePureEvent.h"
ClassImp(BePureEvent);
BePureEvent::BePureEvent() {
}
BePureEvent::~BePureEvent() {
}
void BePureEvent::FillEvent(BeEvent *expevent) {
// fn[NOESEC]; //doubtful
// fCalCorr[NOESEC]; //doubtful
for (Int_t i = 0; i < NOTEL; i++) {
//time information from Tx2
for (Int_t k = 0; k < NOESEC; k++) {
fTau[i][k] = expevent->fTau[i][k];
}//for k NOESEC
//hits
fNOHits[i] = expevent->fNOHits[i];
fNOaHits[i] = expevent->fNOaHits[i];
fNOpHits[i] = expevent->fNOpHits[i];
fNOtHits[i] = expevent->fNOtHits[i];
for (Int_t j = 0; j < NOSELHITS; j++) {
fSN[i][j] = expevent->fSN[i][j];
fTauSN[i][j] = expevent->fTauSN[i][j];
fRN[i][j] = expevent->fRN[i][j];
fE1Sec[i][j] = expevent->fE1Sec[i][j];
fE1Ring[i][j] = expevent->fE1Ring[i][j];
fE2Sec[i][j] = expevent->fE2Sec[i][j];
fE3Sec[i][j] = expevent->fE3Sec[i][j];
fTheta[i][j] = expevent->fTheta[i][j];
fPhi[i][j] = expevent->fPhi[i][j];
fhTau[i][j] = expevent->fhTau[i][j];
fT[i][j] = expevent->fT[i][j];
fP[i][j] = expevent->fP[i][j];
fID[i][j] = expevent->fID[i][j];
}//for j NOHITS
}//for i NOTEL
for (Int_t i = 0; i < 2*NOSELHITS; i++) {
fPM[i] = expevent->fPM[i];
}
fNOA = expevent->fNOA;
fNOP = expevent->fNOP;
//laboratory system
//6Be
fBeLab = expevent->f6BeLab;
fBeThetaLab = expevent->f6BeThetaLab;
fBePhiLab = expevent->f6BePhiLab;
fBePcLab = expevent->f6BePcLab;
fBeIM = expevent->f6BeIM;
//neutron
fNLab = expevent->fNLab;
fNThetaLab = expevent->fNThetaLab;
fNPhiLab = expevent->fNPhiLab;
fNPcLab = expevent->fNPcLab;
//proton1
fP1Lab = expevent->fP1Lab;
fP1ThetaLab = expevent->fP1ThetaLab;
fP1PhiLab = expevent->fP1PhiLab;
fP1PcLab = expevent->fP1PcLab;
//proton2
fP2Lab = expevent->fP2Lab;
fP2ThetaLab = expevent->fP2ThetaLab;
fP2PhiLab = expevent->fP2PhiLab;
fP2PcLab = expevent->fP2PcLab;
//alpha
fALab = expevent->fALab;
fAThetaLab = expevent->fAThetaLab;
fAPhiLab = expevent->fAPhiLab;
fAPcLab = expevent->fAPcLab;
//protons relative kinetic energy
fPPTrel = expevent->fPPTrel;
//CM system 6Li-p
//6Be
fBeCM1 = expevent->f6BeCM1;
fBeThetaCM1 = expevent->f6BeThetaCM1;
fBePhiCM1 = expevent->f6BePhiCM1;
fBePcCM1 = expevent->f6BePcCM1;
//CM system 6Be
//6Be
fBeCM = expevent->f6BeCM;
fBeThetaCM = expevent->f6BeThetaCM;
fBePhiCM = expevent->f6BePhiCM;
fBePcCM = expevent->f6BePcCM;
//neutron
fNCM = expevent->fNCM;
fNThetaCM = expevent->fNThetaCM;
fNPhiCM = expevent->fNPhiCM;
fNPcCM = expevent->fNPcCM;
//proton1
fP1CM = expevent->fP1CM;
fP1ThetaCM = expevent->fP1ThetaCM;
fP1PhiCM = expevent->fP1PhiCM;
fP1PcCM = expevent->fP1PcCM;
//proton2
fP2CM = expevent->fP2CM;
fP2ThetaCM = expevent->fP2ThetaCM;
fP2PhiCM = expevent->fP2PhiCM;
fP2PcCM = expevent->fP2PcCM;
//alpha
fACM = expevent->fACM;
fAThetaCM = expevent->fAThetaCM;
fAPhiCM = expevent->fAPhiCM;
fAPcCM = expevent->fAPcCM;
//general information
fQLiP = expevent->fQLiP;
fTpp = expevent->fTpp;
fTapp = expevent->fTapp;
fCosThetaTk = expevent->fCosThetaTk;
fTap = expevent->fTap;
fTpap = expevent->fTpap;
fCosThetaYk = expevent->fCosThetaYk;
return;
}
void BePureEvent::Reset() {
//temporarily useless
return;
}
/*
* BePureEvent.h
*
* Created on: 6.6.2011
* Author: vratik
*/
#ifndef BEPUREEVENT_H_
#define BEPUREEVENT_H_
#pragma once
#include <TObject.h>
#include <TLorentzVector.h>
#include "BeEvent.h"
#define NOESEC 16 //number of sectors carrying full dE-E information
#define NOSELHITS 6 //maximal number of hits //navrhuji 8 nebo 3 (zapisujeme pouze udalosti s Be)
#define NOTEL 2 //number of telescopes
class BePureEvent {
public:
// Int_t fTrigger;
// Int_t fn[NOESEC]; //doubtful
// Double_t fCalCorr[NOESEC]; //doubtful
//time information from Tx2
Double_t fTau[NOTEL][NOESEC]; //doubtful //time in ns
//hits
Int_t fNOHits[NOTEL]; //general number of hits searched by energy with dE-E, phi and theta information
Int_t fNOaHits[NOTEL]; //hits searched by energy
Int_t fNOpHits[NOTEL]; //hits searched by energy
Int_t fNOtHits[NOTEL]; //hits searched by time
Int_t fSN[NOTEL][NOSELHITS]; //doubtful
Int_t fTauSN[NOTEL][NOSELHITS]; //doubtful
Int_t fRN[NOTEL][NOSELHITS]; //doubtful
Double_t fE1Sec[NOTEL][NOSELHITS]; //hit energy in 1st layer (sectors) of particular telescope
Double_t fE1Ring[NOTEL][NOSELHITS]; //hit energy in 1st layer (rings) of particular telescope
Double_t fE2Sec[NOTEL][NOSELHITS]; //hit energy in 2nd layer of particular telescope
Double_t fE3Sec[NOTEL][NOSELHITS]; //hit energy in 3rd layer of particular telescope
Double_t fTheta[NOTEL][NOSELHITS];
Double_t fPhi[NOTEL][NOSELHITS];
Double_t fhTau[NOTEL][NOSELHITS]; //hit time
Double_t fT[NOTEL][NOSELHITS]; //hit energy
Double_t fP[NOTEL][NOSELHITS]; //hit impulse
Int_t fID[NOTEL][NOSELHITS]; //particle identificator = ZA; proton==11, alpha==24
Int_t fPM[2*NOSELHITS]; //particle markers
// //6Be information
Int_t fNOA; //number of detected alphas related to triple coincidence Be decay
Int_t fNOP; //number of detected protons related to triple coincidence Be decay
//laboratory system
//6Be
TLorentzVector fBeLab;
Double_t fBeThetaLab; //Theta_lab in rad
Double_t fBePhiLab; //Phi_lab in rad
Double_t fBePcLab; //abs(pc) in MeV
Double_t fBeIM; //invariant mass in MeV
//neutron
TLorentzVector fNLab;
Double_t fNThetaLab; //Theta_lab in rad
Double_t fNPhiLab; //Phi_lab in rad
Double_t fNPcLab; //invariant mass in MeV
//proton1
TLorentzVector fP1Lab;
Double_t fP1ThetaLab; //Theta_lab in rad
Double_t fP1PhiLab; //Phi_lab in rad
Double_t fP1PcLab; //invariant mass in MeV
//proton2
TLorentzVector fP2Lab;
Double_t fP2ThetaLab; //Theta_lab in rad
Double_t fP2PhiLab; //Phi_lab in rad
Double_t fP2PcLab; //invariant mass in MeV
//alpha
TLorentzVector fALab;
Double_t fAThetaLab; //Theta_lab in rad
Double_t fAPhiLab; //Phi_lab in rad
Double_t fAPcLab; //invariant mass in MeV
//protons relative kinetic energy
Double_t fPPTrel;
//CM system 6Li-p
//6Be
TLorentzVector fBeCM1;
Double_t fBeThetaCM1; //Theta_cm in rad
Double_t fBePhiCM1; //Phi_lab in rad
Double_t fBePcCM1; //abs(pc) in MeV
//CM system 6Be
//6Be
TLorentzVector fBeCM;
Double_t fBeThetaCM; //Theta_cm in rad
Double_t fBePhiCM; //Phi_lab in rad
Double_t fBePcCM; //abs(pc) in MeV
//neutron
TLorentzVector fNCM;
Double_t fNThetaCM; //Theta_cm in rad
Double_t fNPhiCM; //Phi_lab in rad
Double_t fNPcCM; //abs(pc) in MeV
//proton1
TLorentzVector fP1CM;
Double_t fP1ThetaCM; //Theta_cm in rad
Double_t fP1PhiCM; //Phi_lab in rad
Double_t fP1PcCM; //abs(pc) in MeV
//proton2
TLorentzVector fP2CM;
Double_t fP2ThetaCM; //Theta_cm in rad
Double_t fP2PhiCM; //Phi_lab in rad
Double_t fP2PcCM; //abs(pc) in MeV
//alpha
TLorentzVector fACM;
Double_t fAThetaCM; //Theta_cm in rad
Double_t fAPhiCM; //Phi_lab in rad
Double_t fAPcCM; //abs(pc) in MeV
//general information
Double_t fQLiP; //Q of the reaction 6Li+p-->6Be+n
Double_t fTpp; //"T" system
Double_t fTapp; //"T" system
Double_t fCosThetaTk; //"T" system
Double_t fTap; //"Y" system
Double_t fTpap; //"Y" system
Double_t fCosThetaYk; //"Y" system
BePureEvent();
virtual ~BePureEvent();
ClassDef(BePureEvent, 1);
void FillEvent(BeEvent *expevent);
void Reset();
};
#endif /* BEPUREEVENT_H_ */
/*
* BeReaction.cpp
*
* Created on: 24.5.2010
* Author: Vratislav
*/
#include "BeReaction.h"
ClassImp(BeReaction);
Int_t ranseed = 1;
//TRandom3 BeReaction::ranTheta(1277372118);
//TRandom3 BeReaction::ranTheta(0);
//TRandom3 BeReaction::ranMass(1277372118);
TRandom3 BeReaction::ranMass(0);
TRandom3 BeReaction::ranTheta(0);
Double_t BeReaction::LipM[] = {0};
Double_t BeReaction::BeDecayM[] = {0};
Double_t BeReaction::BeStateEnergy[] = {0};
Double_t BeReaction::DipDecayM[] = {0};
TF1 BeReaction::theta("thetaCM", "TMath::Sin(x)", 0., TMath::Pi());
TF1 BeReaction::thetaUniform("thetaCM", "1", 0., TMath::Pi());
TF1 BeReaction::dipEnergy("dpexcitation", "TMath::Sqrt( x*([0] - x) )", 0., 1.);
BeReaction::BeReaction()
{
ReadParameters();
// PrintParameters();
}
BeReaction::~BeReaction() {
}
void BeReaction::ReadParameters(const char *parameterfile)
{
ifstream parfile;
parfile.open(parameterfile); //parameter file opening
if (!parfile.is_open()) {
Warning("BeReaction::ReadParameters", "File %s opening error\n", parameterfile);
return;
}//if
Double_t value;
TString line;
while (parfile.good()) {
parfile >> value;
line.ReadLine(parfile);
if (line.Contains("//")) line.Resize(line.Index("//"));
line.ToLower();
if (line.Contains("6li")) {
LipM[0] = value;
continue;
}
if (line.Contains("proton")) {
LipM[1] = value;
BeDecayM[2] = 2* value;
DipDecayM[2] = value;
DipDecayM[3] = value;
continue;
}
if (line.Contains("neutron")) {
LipM[3] = value;
continue;
}
if (line.Contains("4he") || line.Contains("alpha")) {
BeDecayM[3] = value;
continue;
}
if (line.Contains("6be_g.s._e")) {
BeStateEnergy[0] = value;
continue;
}
if (line.Contains("6be_g.s._gamma")) {
BeStateEnergy[1] = value;
continue;
}
if (line.Contains("6be_e.s._e")) {
BeStateEnergy[2] = value;
continue;
}
if (line.Contains("6be_e.s._gamma")) {
BeStateEnergy[3] = value;
continue;
}
if (line.Contains("6be_(g.s.)_prob")) {
BeStateEnergy[4] = value;
continue;
}
if (line.Contains("6be_3body_thr")) {
BeStateEnergy[5] = value;
continue;
}
}//while
//set 6Be mass
LipM[2] = BeDecayM[3] /*4He*/ + 2*DipDecayM[3] /*proton*/ - BeStateEnergy[5] /*threshold*/;
BeDecayM[0] = BeDecayM[3] + 2*DipDecayM[3] - BeStateEnergy[5];
parfile.close();
if (parfile.is_open()) { //parameter file closing
Warning("BeReaction::ReadParameters", "File %s closing error\n", parameterfile);
return;
}
Info("BeReaction::ReadParameters", "Mass and other parameters were set");
return;
}
void BeReaction::FillProcess(Double_t _LiT, Double_t _LiThetaIn, Double_t _LiPhi)
{
//simple MC generator using phase volume
Double_t operatingMarray[4] = {0}; //auxiliary array for mass setting
// SetLipMasses(operatingMarray); //to simulate in zero approximation
// SetLipMasses_uniformBeMass(operatingMarray); //to determine energy resolution
SetLipMasses_discreteBeMass(operatingMarray); //to determine energy resolution
fLip.SetName("fLip");
// fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, ThetaCMdistr(0., TMath::Pi()), _LiPhi);
// fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, ThetaCMdistrUniform(0., TMath::Pi()), _LiPhi);
fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, ThetaCMdistrDiscrete(), _LiPhi);
SetBeDecayMasses(operatingMarray);
fBeDecay.SetName("fBeDecay");
fBeDecay.FillReaction(fLip.GetThetaA(), fLip.GetTa(), operatingMarray, ThetaCMdistr(0., TMath::Pi()), fLip.GetPhiA() );
SetDipDecayMasses(operatingMarray);
fDipDecay.SetName("fDipDecay");
fcorrect = fDipDecay.FillReaction( fBeDecay.GetThetaA(), fBeDecay.GetTa(), operatingMarray, ThetaCMdistr(0., TMath::Pi()), fBeDecay.GetPhiA() );
if (fcorrect != 1) {
cout << operatingMarray[0] - 2*DipDecayM[2] << endl;
}
}
void BeReaction::FillProcess(Double_t _LiT, Double_t _LiThetaIn,
Double_t _LiPhi, Double_t *_p_alpha, Double_t *_p_p1, Double_t *_p_p2,
Double_t _thetaCMmin, Double_t _thetaCMmax)
{
Double_t thetaCM = ThetaCMdistr(_thetaCMmin, _thetaCMmax);
FillProcess(_LiT, _LiThetaIn, _LiPhi, _p_alpha, _p_p1, _p_p2, thetaCM);
}
void BeReaction::FillProcess(Double_t _LiT, Double_t _LiThetaIn,
Double_t _LiPhi, Double_t *_p_alpha, Double_t *_p_p1, Double_t *_p_p2,
Double_t _thetaCM)
{
//function variables: LiT, LiThetaIn, LiPhiIn, BeThetaCM, impulses array
//_LiT:
//_thetaCMmin: in rad
//_thetaCMmax: in rad
//physical generator based on theoretical calculations
//set masses for 6Li + p --> 6Be + n: mass of 6Be taken from generator
TLorentzVector pa(_p_alpha[0], _p_alpha[1], _p_alpha[2]); //alpha in the a-p-p CM
pa.SetE( E( Power(pa.P(), 2), BeDecayM[3] ) );
TLorentzVector pp1(_p_p1[0], _p_p1[1], _p_p1[2]); //proton in the a-p-p CM
pp1.SetE( E( Power(pp1.P(), 2), LipM[1] ) );
TLorentzVector pp2(_p_p2[0], _p_p2[1], _p_p2[2]); //proton in the a-p-p CM
pp2.SetE( E( Power(pp2.P(), 2), LipM[1] ) );
TLorentzVector pbe; //6Be in the a-p-p CM
pbe = pa + pp1 + pp2;
Double_t operatingMarray[4] = {0}; //auxilliary array for mass setting
operatingMarray[0] = LipM[0];
operatingMarray[1] = LipM[1];
operatingMarray[2] = pbe.E();
operatingMarray[3] = LipM[3];
//reaction 6Li + p --> 6Be + n: BeThetaCM taken from generator
fLip.SetName("fLip");
// fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, _thetaCM*TMath::DegToRad(), _LiPhi);
fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, _thetaCM/*ThetaCMdistr(_thetaCMmin, _thetaCMmax)*/, _LiPhi); //thetaCM from generator is not used
//transformation from a-p-p CM system into lab system
TRotation r1;
r1.SetZAxis(-GetNeutronP(), GetBeP());
TLorentzRotation rot(r1);
//rotation to CM system where axis are parallel with lab system (CMaux)
pa.Transform(rot);
pp1.Transform(rot);
pp2.Transform(rot);
TVector3 beta = GetBe().BoostVector();
//CMaux --> lab
pa.Boost(beta);
pp1.Boost(beta);
pp2.Boost(beta);
//set the variables fBeDecay.fTb(.ThetaA, .PhiA) (alpha),
// fDipDecay.fTa(.ThetaA, .PhiA) (one of the protons),
// fDipDecay.fTb(.ThetaB, .PhiB) (one of the protons)
//set variables for the alpha particle
fBeDecay.SetM(3, BeDecayM[3]);
fBeDecay.SetPhiB(pa.Phi());
fBeDecay.SetThetaB(pa.Theta());
fBeDecay.SetTb(pa.E()-BeDecayM[3]);
//set variables for the first proton
fDipDecay.SetM(2, DipDecayM[2]);
fDipDecay.SetPhiA(pp1.Phi());
fDipDecay.SetThetaA(pp1.Theta());
fDipDecay.SetTa(pp1.E()-DipDecayM[2]);
//set variables for the second proton
fDipDecay.SetM(3, DipDecayM[3]);
fDipDecay.SetPhiB(pp2.Phi());
fDipDecay.SetThetaB(pp2.Theta());
fDipDecay.SetTb(pp2.E()-DipDecayM[3]);
// Info("BeReaction::FillProcess", "%f\t%f\t%f", pp2.E(), pp2.Phi(), pp2.Theta());
return;
}
void BeReaction::PrintParameters()
{
cout <<endl;
for (Int_t i = 0; i < 4; i++) {
cout << "LipM[" << i << "] is " << LipM[i] << endl;
}
cout <<endl;
for (Int_t i = 0; i < 4; i++) {
cout << "BeDecayM[" << i << "] is " << BeDecayM[i] << endl;
}
cout <<endl;
for (Int_t i = 0; i < 6; i++) {
cout << "BeStateEnergy[" << i << "] is " << BeStateEnergy[i] << endl;
}
cout <<endl;
for (Int_t i = 0; i < 4; i++) {
cout << "DipDecayM[" << i << "] is " << DipDecayM[i] << endl;
}
return;
}
Double_t BeReaction::ThetaCMdistr(Double_t tmin, Double_t tmax) const {
return theta.GetRandom(tmin, tmax);
}
Double_t BeReaction::ThetaCMdistrUniform(Double_t tmin, Double_t tmax) const {
return thetaUniform.GetRandom(tmin, tmax);
}
Double_t BeReaction::ThetaCMdistrDiscrete() const {
const Int_t nopoints = 6;
const Double_t angles[nopoints] = {20*TMath::DegToRad(),
45*TMath::DegToRad(),
70*TMath::DegToRad(),
95*TMath::DegToRad(),
120*TMath::DegToRad(),
150*TMath::DegToRad()};
Int_t angle = ranMass.Integer(nopoints);
return angles[angle];
}
void BeReaction::SetLipMasses(Double_t _lpm[]) {
for (Int_t i = 0; i < 4; i++) _lpm[i] = 0.;
//constant masses setting
_lpm[0] = LipM[0]; //mass of 6Li
_lpm[1] = LipM[1]; //mass of proton
_lpm[3] = LipM[3]; //mass of neutron
//variable 6Be mass setting
if ( BeStateEnergy[4] >= ranMass.Uniform(0., 1.) ) {
while (_lpm[2] <= BeDecayM[3] + 2*DipDecayM[3]) {
_lpm[2] = LipM[2] + ranMass.Gaus(BeStateEnergy[0], BeStateEnergy[1]);
// printf("%f\t>\t%f\n", BeStateEnergy[0], BeStateEnergy[1]);
// printf("%f\t>\t%f\n", _lpm[2], fBeDecay.GetMb()+2*fDipDecay.GetMa());
// _lpm[2] = BeDecayM[3] + 2*DipDecayM[4] + ranMass.Gaus(BeStateEnergy[0], BeStateEnergy[1]);
}//while
}//if
else {
while (_lpm[2] <= BeDecayM[3] + 2*DipDecayM[3]) {
_lpm[2] = LipM[2] + ranMass.Gaus(BeStateEnergy[2], BeStateEnergy[3]);
// printf("%f\t>\t%f\n", _lpm[2], fBeDecay.GetMb()+2*fDipDecay.GetMa());
// _lpm[2] = BeDecayM[3] + 2*DipDecayM[4] + ranMass.Gaus(BeStateEnergy[2], BeStateEnergy[3]);
}//while
}//if
// printf("%f\t%f\t%f\t%f\n", LipM[0], LipM[1], LipM[2], LipM[3]);
return;
}
void BeReaction::SetLipMasses_uniformBeMass(Double_t _lpm[]) {
//constant masses setting
_lpm[0] = LipM[0]; //mass of 6Li
_lpm[1] = LipM[1]; //mass of proton
_lpm[3] = LipM[3]; //mass of neutron
//uniform distrubution 6Be mass setting in range (Emin;Emax)
const Double_t Emin = 0.;
const Double_t Emax = 20.;
// while (_lpm[2] <= fBeDecay.GetMb()+2*fDipDecay.GetMa()) {
_lpm[2] = BeDecayM[3] + 2*DipDecayM[3] + ranMass.Uniform(Emin, Emax);
// _lpm[2] = LipM[2] + BeStateEnergy[5] + ranMass.Uniform(Emin, Emax);
// }
return;
}
void BeReaction::SetLipMasses_discreteBeMass(Double_t _lpm[]) {
//constant masses setting
_lpm[0] = LipM[0]; //mass of 6Li
_lpm[1] = LipM[1]; //mass of proton
_lpm[3] = LipM[3]; //mass of neutron
//discrete distrubution 6Be mass 1.4, 3, 6, 9, 13
const Int_t nopoints = 5;
const Double_t masses[nopoints] = {1.4, 3, 6, 9, 13};
Int_t mass = ranMass.Integer(nopoints);
_lpm[2] = BeDecayM[3] + 2*DipDecayM[3] + masses[mass];
return;
}
void BeReaction::SetBeDecayMasses(Double_t _bdm[]) {
//constant masses setting
_bdm[1] = BeDecayM[1]; //zero mass
_bdm[3] = BeDecayM[3]; //mass of 4He
//variable 6Be mass setting
_bdm[0] = fLip.GetMa();
//variable diproton mass setting
const Double_t freeEnergy = fLip.GetMa() - (BeDecayM[0] + BeStateEnergy[5]);
//TF1 dipEnergy("dpexcitation", "TMath::Sqrt( x*([0] - x) )", 0., freeEnergy);
dipEnergy.SetRange(0., freeEnergy);
dipEnergy.SetParameter(0, freeEnergy);
const Double_t random = dipEnergy.GetRandom();
_bdm[2] = BeDecayM[2] + random; //problem
// if (BeDecayM[2] > _bdm[2]) {
// printf("//////////////////\t%f\t%f\t%f\n", BeDecayM[2], random, freeEnergy);
// printf("//////////////////\t%f\t%f\t%f\t%f\n", fLip.GetMa(), BeDecayM[0], fLip.GetMa() -BeDecayM[0], BeStateEnergy[5]);
// }
}
void BeReaction::SetDipDecayMasses(Double_t _2pm[])
{
if (fBeDecay.GetMa()/2. < 938.272) {
printf("\n\n%f\t%f\n\n", fBeDecay.GetMa()/2., fBeDecay.GetMb());
}
_2pm[0] = fBeDecay.GetMa();
DipDecayM[0] = _2pm[0];
_2pm[1] = 0;
_2pm[2] = DipDecayM[2];
_2pm[3] = DipDecayM[3];
}
TVector3 BeReaction::GetAlphaP()
{
TVector3 vAlpha(1., 1., 1.);
Double_t pc = TMath::Sqrt(TMath::Power(GetAlphaT(), 2) + 2*GetAlphaT()*fBeDecay.GetMb());
vAlpha.SetMagThetaPhi(pc, GetAlphaTheta(), GetAlphaPhi());
return vAlpha;
}
TLorentzVector BeReaction::GetAlpha()
{
// TLorentzVector *vAlpha = new TLorentzVector(GetAlphaP(), GetAlphaT()+fBeDecay.GetMb());
// cout << "getalfa Zacatek" << endl;
// cout << GetAlphaT()+fBeDecay.GetMb() << endl;
// cout << fBeDecay.GetMb() << endl;
// cout << GetAlphaT() << endl;
// cout << "getalfa konec" << endl;
return TLorentzVector(GetAlphaP(), GetAlphaT()+fBeDecay.GetMb());
}
TVector3 BeReaction::GetP1P()
{
TVector3 vP1(1., 1., 1.);
Double_t pc = TMath::Sqrt(TMath::Power(GetP1T(), 2) + 2*GetP1T()*fDipDecay.GetMa());
vP1.SetMagThetaPhi(pc, GetP1Theta(), GetP1Phi());
return vP1;
}
TLorentzVector BeReaction::GetP1()
{
// TLorentzVector *vP1 = new TLorentzVector(*GetP1P(), GetP1T()+fDipDecay.GetMa());
return TLorentzVector(GetP1P(), GetP1T()+fDipDecay.GetMa());
}
TVector3 BeReaction::GetP2P()
{
TVector3 vP2(1., 1., 1.);
Double_t pc = TMath::Sqrt(TMath::Power(GetP2T(), 2) + 2*GetP2T()*fDipDecay.GetMb());
vP2.SetMagThetaPhi(pc, GetP2Theta(), GetP2Phi());
return vP2;
}
TLorentzVector BeReaction::GetP2()
{
// TLorentzVector *vP2 = new TLorentzVector(*GetP2P(), GetP2T()+fDipDecay.GetMb());
// Info("BeReaction::GetP2", "\t%f\t%f\t%f", GetP2P().Mag(), GetP2T(), GetP2Theta());
return TLorentzVector(GetP2P(), GetP2T()+fDipDecay.GetMb());
}
//TVector3* BeReaction::GetNeutronP()
TVector3 BeReaction::GetNeutronP()
{
// TVector3 *vN = new TVector3(1., 1., 1.);
TVector3 vN(1., 1., 1.);
Double_t pc = TMath::Sqrt(TMath::Power(GetNeutronT(), 2) + 2*GetNeutronT()*fLip.GetMb());
// vN->SetMagThetaPhi(pc, GetNeutronTheta(), GetNeutronPhi());
vN.SetMagThetaPhi(pc, GetNeutronTheta(), GetNeutronPhi());
return vN;
}
TLorentzVector BeReaction::GetNeutron()
{
// TLorentzVector *vN = new TLorentzVector(*GetNeutronP(), GetNeutronT()+fLip.GetMb());
// TLorentzVector vN = new TLorentzVector(GetNeutronP(), GetNeutronT()+fLip.GetMb());
// return vN;
return TLorentzVector(GetNeutronP(), GetNeutronT()+fLip.GetMb());
}
TVector3 BeReaction::GetBeP()
{
TVector3 vBe(1., 1., 1.);
Double_t pc = TMath::Sqrt(TMath::Power(GetBeT(), 2) + 2*GetBeT()*fLip.GetMa());
vBe.SetMagThetaPhi(pc, GetBeTheta(), GetBePhi());
return vBe;
}
TLorentzVector BeReaction::GetBe()
{
return TLorentzVector(GetBeP(), GetBeT()+fLip.GetMa());
}
Double_t BeReaction::T(Double_t pc2, Double_t mc2)
{
return TMath::Power(mc2*mc2 + pc2, 0.5) - mc2;
}
Double_t BeReaction::E(Double_t pc2, Double_t mc2)
{
return TMath::Power(mc2*mc2 + pc2, 0.5);
}
void BeReaction::Reset()
{
fLip.Reset();
fBeDecay.Reset();
fDipDecay.Reset();
fcorrect = 0;
return;
}
/*
* BeReaction.h
*
* Created on: 24.5.2010
* Author: Vratislav
*/
//#ifndef BEREACTION_H_
//#define BEREACTION_H_
#pragma once
//#include "DllExport.h"
#include <TObject.h>
#include <TROOT.h>
#include <TRandom3.h>
#include <TF1.h>
#include <TString.h>
#include <TMath.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <TLorentzRotation.h>
#include <iostream> //k nicemu
#include <fstream> //k nicemu
#include <iomanip> //k nicemu
#include <sstream> //k nicemu
//#include <string>
#include "BinaryReaction.h"
//#define PARAMETERFILE "sim.par"
using std::cout;
using std::endl;
//using std::string;
//using std::getline;
class BeReaction {
//private:
public:
//1 + 2 --> a + b
BinaryReaction fLip; //6Li + p --> 6Be + n
BinaryReaction fBeDecay; //6Be --> 2p + 4He
BinaryReaction fDipDecay; //2p --> p + p
Int_t fcorrect;
//neni dobre, ze jsou static
static Double_t LipM[4]; //! {E_6Li; E_proton; E_6Be; E_neutron}
static Double_t BeDecayM[4]; //! {E_6Be_g.s.; 0; E_diproton; E_4He}
static Double_t BeStateEnergy[6]; //! {E_g.s. = 0; Gamma_g.s.; E_e.s.; Gamma_e.s.; population ratio; 4He2p threshold}
static Double_t DipDecayM[4]; //! {E_diproton; 0; proton; proton}
// static TRandom3 ranTheta;
static TRandom3 ranMass; //!
static TRandom3 ranTheta; //!
static TF1 theta; //!
static TF1 thetaUniform; //!
static TF1 dipEnergy; //! phase volume
public:
BeReaction();
virtual ~BeReaction();
ClassDef(BeReaction, 1);
void FillProcess(Double_t _LiT, Double_t _LiThetaIn, Double_t _LiPhi);
//using generator
void FillProcess(Double_t _LiT, Double_t _LiThetaIn, Double_t _LiPhi,
Double_t *_p_alpha, Double_t *_p_p1, Double_t *_p_p2,
Double_t _thetaCM);
void FillProcess(Double_t _LiT, Double_t _LiThetaIn, Double_t _LiPhi,
Double_t *_p_alpha, Double_t *_p_p1, Double_t *_p_p2,
Double_t _thetaCMmin, Double_t _thetaCMmax);
// void ReadParameters();
void ReadParameters(const char *parameterfile = "sim.par");
void PrintParameters();
void Reset();
Double_t ThetaCMdistr(Double_t tmin, Double_t tmax) const;
Double_t ThetaCMdistrUniform(Double_t tmin, Double_t tmax) const;
Double_t ThetaCMdistrDiscrete() const;
void SetLipMasses(Double_t _lpm[]);
void SetLipMasses_uniformBeMass(Double_t _lpm[]);
void SetLipMasses_discreteBeMass(Double_t _lpm[]);
void SetBeDecayMasses(Double_t _bdm[]);
void SetDipDecayMasses(Double_t _2pm[]);
//getters in laboratory system
//alpha
Double_t GetAlphaT() const { return fBeDecay.GetTb(); };
Double_t GetAlphaTheta() const { return fBeDecay.GetThetaB(); };
Double_t GetAlphaPhi() const { return fBeDecay.GetPhiB(); };
// TVector3* GetAlphaP();
TVector3 GetAlphaP();
TLorentzVector GetAlpha();
//first proton (higher T)
Double_t GetP1T() const { return (fDipDecay.GetTa() >= fDipDecay.GetTb()) ? fDipDecay.GetTa() : fDipDecay.GetTb(); };
Double_t GetP1Theta() const { return (fDipDecay.GetTa() >= fDipDecay.GetTb()) ? fDipDecay.GetThetaA() : fDipDecay.GetThetaB(); };
Double_t GetP1Phi() const { return (fDipDecay.GetTa() >= fDipDecay.GetTb()) ? fDipDecay.GetPhiA() : fDipDecay.GetPhiB(); };
TVector3 GetP1P();
TLorentzVector GetP1();
//second proton (lower T)
Double_t GetP2T() const { return (fDipDecay.GetTa() <= fDipDecay.GetTb()) ? fDipDecay.GetTa() : fDipDecay.GetTb(); };
Double_t GetP2Theta() const { return (fDipDecay.GetTa() <= fDipDecay.GetTb()) ? fDipDecay.GetThetaA() : fDipDecay.GetThetaB(); };
Double_t GetP2Phi() const { return (fDipDecay.GetTa() <= fDipDecay.GetTb()) ? fDipDecay.GetPhiA() : fDipDecay.GetPhiB(); };
TVector3 GetP2P();
TLorentzVector GetP2();
//neutron
Double_t GetNeutronT() const { return fLip.GetTb(); };
Double_t GetNeutronTheta() const { return fLip.GetThetaB(); };
Double_t GetNeutronPhi() const { return fLip.GetPhiB(); };
// TVector3* GetNeutronP();
TVector3 GetNeutronP();
TLorentzVector GetNeutron();
//6Be
Double_t GetBeT() const { return fLip.GetTa(); };
Double_t GetBeTheta() const { return fLip.GetThetaA(); };
Double_t GetBePhi() const { return fLip.GetPhiA(); };
TVector3 GetBeP();
TLorentzVector GetBe();
//getters
BinaryReaction* GetLip() { return &fLip; };
BinaryReaction* GetBeDecay() { return &fBeDecay; };
BinaryReaction* GetDipDecay() { return &fDipDecay; };
//diagnostics
UInt_t GetMSeed() const { return ranMass.GetSeed(); };
// UInt_t GetTSeed() const { return ranTheta.GetSeed(); };
UInt_t GetLipPhiSeed() const { return fLip.GetPSeed(); };
UInt_t GetBePhiSeed() const { return fBeDecay.GetPSeed(); };
UInt_t GetDipPhiSeed() const { return fDipDecay.GetPSeed(); };
//physics
Double_t T(Double_t pc2, Double_t mc2); //v podstate taktez
Double_t E(Double_t pc2, Double_t mc2); //v podstate taktez
};
//#endif /* BEREACTION_H_ */
/*
* BeWork.cpp
*
* Created on: 2.2.2010
* Author: Vratislav
*/
//unsigned int max_number_threads = 2;
#include "BeWork.h"
//#include "BeThread.h"
ClassImp(BeWork);
#ifndef __CINT__
//boost::mutex gFillMutex;
//boost::mutex gFillMutex;
#endif /* __CINT __ */
BeWork::BeWork()
{
Info("BeWork::BeWork", "CsI resolution is %f for alphas and %f for protons", 100*fCsIResA, 100*fCsIResP);
Info("BeWork::BeWork", "Si resolution is %f keV", fSiRes*1000);
Info("BeWork::BeWork", "Creating TELoss objects ...");
CreateTELosses();
}
BeWork::BeWork(const char* configfile)
{
fConfigFile = configfile;
ReadConfigFile();
SetWorkDir();
Info("BeWork::BeWork", "Work directory: \"%s\"", fWorkDir.Data());
Info("BeWork::BeWork", "CsI resolution is %f for alphas and %f for protons", 100*fCsIResA, 100*fCsIResP);
Info("BeWork::BeWork", "Si resolution is %f keV", fSiRes*1000);
Info("BeWork::BeWork", "Creating TELoss objects ...");
CreateTELosses();
// fCalibFile = "/data2/be/parameterfiles";
Info("BeWork", "End of second constructor\n\n");
}
BeWork::BeWork(BeWork &bework) {
//copy constructor
fSiAlpha = new TELoss( *(bework.fSiAlpha) );
fSiP = new TELoss( *(bework.fSiP) );
fCsIAlpha = new TELoss( *(bework.fCsIAlpha) );
fCsIP = new TELoss( *(bework.fCsIP) );
fTargetAlpha = new TELoss( *(bework.fTargetAlpha) );
fTargetP = new TELoss( *(bework.fTargetP) );
fTargetLi = new TELoss( *(bework.fTargetLi) );
fTargetWinAlpha = new TELoss( *(bework.fTargetWinAlpha) );
fTargetWinP = new TELoss( *(bework.fTargetWinP) );
fTargetWinLi = new TELoss( *(bework.fTargetWinLi) );
//configuration and parameter files
fConfigFile = bework.fConfigFile;
fWorkDir = bework.fWorkDir;
fRawFilePath = bework.fRawFilePath;
fSiRes = bework.fSiRes;
fCsIResP = bework.fCsIResP;
fCsIResA = bework.fCsIResA;
fTBeamMC = bework.fTBeamMC;
fTBeamResMC = bework.fTBeamResMC;
fBeamX_MC = bework.fBeamX_MC;
fBeamY_MC = bework.fBeamY_MC;
fBeamX_sigma_MC = bework.fBeamX_sigma_MC;
fBeamY_sigma_MC = bework.fBeamY_sigma_MC;
//
fT1SimPosition = bework.fT1SimPosition;
fT2SimPosition = bework.fT2SimPosition;
fX11_FD = bework.fX11_FD;
fX11 = bework.fX11;
fX11_BD = bework.fX11_BD;
fX12_FD = bework.fX12_FD;
fX12 = bework.fX12;
fX12_BD = bework.fX12_BD;
fX13_FD = bework.fX13_FD;
fX21_FD = bework.fX21_FD;
fX21 = bework.fX21;
fX21_BD = bework.fX21_BD;
fX22_FD = bework.fX22_FD;
fX22 = bework.fX22;
fX22_BD = bework.fX22_BD;
fX23_FD = bework.fX23_FD;
Info("BeWork::BeWork", "Copy constructor finished");
}
BeWork::~BeWork()
{
Info("BeWork::~BeWork", "Destructor called");
delete fSiAlpha;
delete fSiP;
delete fCsIAlpha;
delete fCsIP;
delete fTargetAlpha;
delete fTargetP;
delete fTargetLi;
delete fTargetWinAlpha;
delete fTargetWinP;
delete fTargetWinLi;
Info("BeWork::~BeWork", "Destructor finished");
}
void BeWork::ReadConfigFile() {
ifstream cfile;
//generator file opening
if (!fConfigFile.IsNull()) {
cfile.open(fConfigFile.Data());
if (!cfile.is_open()) {
Error("BeWork::ReadConfigFile", "Configuration file was not opened!!!");
return;
}//if
}//if
TString line;
TString configuration;
while (cfile.good()) {
//read line from file
line.ReadLine(cfile);
if (line.Contains("//")) line.Resize(line.Index("//"));
line.Append(" ");
configuration.Append(line);
}
// cout << configuration.Data() << endl;
//lower and upper case are important
ConfigDictionary CDpath(configuration.Data());
fRawFilePath = CDpath.GetString("rawFilePath");
// fCutFile = CDpath.GetString("cutFile");
//numbers only
configuration.ToLower();
ConfigDictionary CD(configuration.Data());
fBeOnly = CD.GetBool("beonly");
fsRatioMax = CD.GetDouble("sratiomax");
fsRatioMin = CD.GetDouble("sratiomin");
// printf("%d\t%d", fBeOnly, CD.GetBool("beonly"));
fSiRes = CD.GetDouble("sires");
fCsIResP = CD.GetDouble("csiresp");
fCsIBestResP = CD.GetDouble("csirespmin");
fCsIResA = CD.GetDouble("csiresa")/1.5;
fCsIBestResA = CD.GetDouble("csiresamin");
fTBeamMC = CD.GetDouble("tbeam");
fTBeamResMC = CD.GetDouble("tbeamres");
fBeamX_MC = CD.GetDouble("beamx");
fBeamY_MC = CD.GetDouble("beamy");
fBeamX_sigma_MC = CD.GetDouble("beamx_sigma");
fBeamY_sigma_MC = CD.GetDouble("beamy_sigma");
fT1SimPosition = CD.GetDouble("t1simposition");
fT2SimPosition = CD.GetDouble("t2simposition");
//detector thicknesses
fX11_FD = CD.GetDouble("t11_frontdead");
fX11 = CD.GetDouble("t11");
fX11_BD = CD.GetDouble("t11_backdead");
fX12_FD = CD.GetDouble("t12_frontdead");
fX12 = CD.GetDouble("t12");
fX12_BD = CD.GetDouble("t12_backdead");
//#define XD13F 14. //fixme original value probably in Si equivalent
fX13_FD = CD.GetDouble("t13_frontdead");
fX21_FD = CD.GetDouble("t21_frontdead");
fX21 = CD.GetDouble("t21");
fX21_BD = CD.GetDouble("t21_backdead");
fX22_FD = CD.GetDouble("t22_frontdead");
fX22 = CD.GetDouble("t22");
fX22_BD = CD.GetDouble("t22_backdead");
//#define XD23F 14. //fixme original value probably in Si equivalent
fX23_FD = CD.GetDouble("t23_frontdead");
//config file closing
cfile.close();
if (cfile.is_open()) {
Warning("BeWork::ReadConfigFile", "File %s closing error\n", fConfigFile.Data());
}//if
return;
}
void BeWork::SetWorkDir() {
fWorkDir = fConfigFile;
fWorkDir.Remove(fConfigFile.Last('/'), fConfigFile.Length());
}
const char* BeWork::GetWorkDir() {
return fWorkDir.Data();
}
Int_t BeWork::FillExpFile(const char* inputrawfile, const char* outputfile,
Long64_t noevents, Option_t *opt)
{
//convert raw data in channel units into all needed kinematical units, new file containing tree "be" filled by "beevent" is created
//
//inputrawfile: .root file containing raw data
//noevents: number of events to be filled, if noevents is set to 0, whole file will be processed
//calibfilepath: path to files with calibration data "beSi.cal", "beCsIa.cal" and "beCsIp.cal"
//outputfile: .root file with calibrated data
//option: config: add config file name into the outputfile name
TString ofile;
// if (opt
ofile.Form("%s/%s", fWorkDir.Data(), outputfile);
TString ifile;
ifile.Form("%s/%s", fRawFilePath.Data(), inputrawfile);
//read calibration parameters and thresholds
TFile *fr = new TFile(ifile.Data(), "READ");
if (fr->IsOpen() == 0) {
Warning("BeWork::FillExpFile", "File %s was not opened and won't be processed", ifile.Data());
return 0;
}
//RAW data tree opening
TTree *tr = (TTree*)fr->Get("RAW");
if (!tr) {
Warning("BeWork::FillExpFile", "Tree \"RAW\" was not found in file %s .", ifile.Data());
return 0;
}
AculRaw *eventr = new AculRaw();
Long64_t noEntries = tr->GetEntries(/*getCondition*/);
if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
tr->SetBranchAddress("channels", &eventr);
//output file containing tree creating
TFile *fw = new TFile(ofile.Data(), "RECREATE"); //dodelat overeni vytvoreni
BeEvent *eventw = new BeEvent(fConfigFile.Data());
TTree *tw = new TTree("be", "strom s energetickou a jinou informaci"); //predelat jmena, ...
tw->Bronch("BePhys", "BeEvent", &eventw, 1024000, 99); //large buffer needed for faster processing in the course of analysis
Info("BeWork::FillExpFile", "In file %s %d events will be processed", ifile.Data(), (int)noEntries);
if (fBeOnly) Info("BeWork::FillExpFile", "Only events containing Be will be saved");
else Info("BeWork::FillExpFile", "All events will be saved");
Info("BeWork::FillExpFile", "Output file %s will be filled", ofile.Data());
fw->cd();
for (Int_t i = 0; i <= (Int_t)noEntries; i++) {
printProgBar(i, noEntries);
fr->cd();
tr->GetEntry(i);
eventw->FillEvent(eventr);
fw->cd();
if (fBeOnly) {
if (eventw->WriteCondition()) {
tw->Fill();
} //zkontrolovat
} else
tw->Fill();
}//for
printf("\n\n");
fw->cd();
tw->AutoSave();
fw->cd();
tw->Write();
fw->Close();
return (Int_t)noEntries;
}
Int_t BeWork::FillBeFile(const char* inputfile, Long64_t noevents, const char* rtree, const char* rbranch,
const char* outputfile, const char* wtree, const char* wbranch, Option_t *opt) {
//inputfile: .root file containing calibrated data
//noevents: number of events to be filled, if noevents is set to 0, whole file will be processed
//rtree: name of the tree containing information with structure BeEvent class to be read
//rbranch: name of the branch to be read
//outputfile: .root file with data containing Be event only, which will be filled into TTree named "beonly"
//wtree: name of the branch to be wrote
//wbranch: name of the tree containing information with structure BePureEvent class to be wrote
//option:
//output file name creating
if (!strlen(outputfile)) {
Error("BeWork::FillBeFile", "Name of outputfile must be set");
return 0;
}
TString output(outputfile);
TFile fr(inputfile, "READ");
if (fr.IsOpen() == 0) {
Warning("BeWork::FillBeFile", "File %s was not opened and won't be processed", inputfile);
return 0;
}
//RAW data tree opening
TTree *tr = (TTree*)fr.Get(rtree);
if (!tr) {
Warning("BeWork::FillBeFile", "Tree \"be\" was not found in file %s .", inputfile);
return 0;
}
BeEvent *eventr = new BeEvent(fConfigFile.Data());
Long64_t noEntries = tr->GetEntries(/*getCondition*/);
if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
tr->SetBranchAddress(rbranch, &eventr);
//tree with simulated beam coordinates opening
TTree *trbeam = 0;
Double_t x = 0;
Double_t y = 0;
Double_t z = 0;
Double_t thetaMC = 0;
Double_t E_IM = 0;
TLorentzVector *beamLab = new TLorentzVector();
TLorentzVector *alphaLab = new TLorentzVector();
TLorentzVector *p1Lab = new TLorentzVector();
TLorentzVector *p2Lab = new TLorentzVector();
TLorentzVector *alphaCM = new TLorentzVector();
TLorentzVector *p1CM = new TLorentzVector();
TLorentzVector *p2CM = new TLorentzVector();
TLorentzVector *beCM = new TLorentzVector();
Double_t sTpp;
Double_t sTapp;
Double_t sCosThetaTk;
TString treereadname(rtree);
if (treereadname.Contains("sim")) {
Info("BeWork::FillBeFile", "Tree containing information about simulated beam opening");
trbeam = (TTree*)fr.Get("sbeam");
if (!trbeam) {
Warning("BeWork::FillBeFile", "Tree \"sbeam\" was not found in file %s .", inputfile);
return 0;
}
trbeam->SetBranchAddress("bx", &x);
trbeam->SetBranchAddress("by", &y);
trbeam->SetBranchAddress("bz", &z);
trbeam->SetBranchAddress("thetaMC", &thetaMC);
trbeam->SetBranchAddress("E_IM", &E_IM);
trbeam->SetBranchAddress("beamLab.", &beamLab);
trbeam->SetBranchAddress("sALab.", &alphaLab);
trbeam->SetBranchAddress("sP1Lab.", &p1Lab);
trbeam->SetBranchAddress("sP2Lab.", &p2Lab);
trbeam->SetBranchAddress("sACM.", &alphaCM);
trbeam->SetBranchAddress("sP1CM.", &p1CM);
trbeam->SetBranchAddress("sP2CM.", &p2CM);
trbeam->SetBranchAddress("sBeCM.", &beCM);
trbeam->SetBranchAddress("sTpp", &sTpp);
trbeam->SetBranchAddress("sTapp", &sTapp);
trbeam->SetBranchAddress("sCosThetaTk", &sCosThetaTk);
}//if
//output file containing tree creating
TFile fw(output.Data(), "RECREATE"); //dodelat overeni vytvoreni
BePureEvent *eventw = new BePureEvent();
TTree *tw = new TTree(wtree, "tree containing Be events only"); //predelat jmena, ...
tw->Bronch(wbranch, "BePureEvent", &eventw, 1024000, 99); //large buffer needed for faster processing in the course of analysis
//tree containing beam information
TTree *twbeam = 0;
if (treereadname.Contains("sim")) {
Info("BeWork::FillBeFile", "Tree containing information about simulated beam creating");
twbeam = new TTree("sbeam", "tree with beam position and ThetaCM read from generator (not used for simulation)");
twbeam->Branch("bx", &x, "bx/D");
twbeam->Branch("by", &y, "by/D");
twbeam->Branch("bz", &z, "bz/D");
twbeam->Branch("thetaMC", &thetaMC, "thetaMC/D");
twbeam->Branch("E_IM", &E_IM, "E_IM/D");
twbeam->Bronch("beamLab", "TLorentzVector", &beamLab, 1024000, 99);
twbeam->Bronch("sALab.", "TLorentzVector", &alphaLab, 1024000, 99);
twbeam->Bronch("sP1Lab.", "TLorentzVector", &p1Lab, 1024000, 99);
twbeam->Bronch("sP2Lab.", "TLorentzVector", &p2Lab, 1024000, 99);
twbeam->Bronch("sACM.", "TLorentzVector", &alphaCM, 1024000, 99);
twbeam->Bronch("sP1CM.", "TLorentzVector", &p1CM, 1024000, 99);
twbeam->Bronch("sP2CM.", "TLorentzVector", &p2CM, 1024000, 99);
twbeam->Bronch("sBeCM.", "TLorentzVector", &beCM, 1024000, 99);
twbeam->Branch("sTpp", &sTpp, "sTpp/D");
twbeam->Branch("sTapp", &sTapp, "sTapp/D");
twbeam->Branch("sCosThetaTk", &sCosThetaTk, "sCosThetaTk/D");
}
Info("BeWork::FillBeFile", "In file %s %d events will be processed", inputfile, (int)noEntries);
Info("BeWork::FillBeFile", "Output file %s will be filled", output.Data());
for (Int_t i = 0; i < (Int_t)noEntries; i++) {
printProgBar(i, (Int_t)noEntries);
fr.cd();
tr->GetEntry(i);
if (trbeam) trbeam->GetEntry(i);
if ( eventr->WriteConditionPure() ) {
eventw->FillEvent(eventr);
fw.cd();
tw->Fill();
if (twbeam) twbeam->Fill();
}//if
}//for
Info("BeWork::FillBeFile", "%d events containing Be were found", (int)tw->GetEntries());
printf("\n\n");
fw.cd();
tw->AutoSave();
Info("BeWork::FillBeFile", "Tree \"%s\" with Be events saving", tw->GetName());
tw->Write();
if (twbeam) {
Info("BeWork::FillBeFile", "Tree \"%s\" with beam information saving", twbeam->GetName());
twbeam->Write();
}
delete eventw;
fw.Close();
fr.cd();
delete eventr;
fr.Close();
return (Int_t)noEntries;
}
Int_t BeWork::FillBeExpFile(const char *inputfile, const char* outputfile, Long64_t noevents) {
//fill file containing tree with BePureEvent class structure
//inputfile:
//noevents:
//outputfile: if empty, character sequence "Exp" in inputfile will be changed to "Be"
TString ifile;
ifile.Form("%s/%s", fWorkDir.Data(), inputfile);
TString ofile;
ofile.Form("%s/%s", fWorkDir.Data(), outputfile);
// TString output(inputfile);
// if (strlen(outputfile)) {
// output = outputfile;
// }
// else {
// output.ReplaceAll("Exp", "Be");
// }
const char* rtreename = "be";
const char* rbranchname = "BePhys";
const char* wtreename = "beonly";
const char* wbranchname = "BeEvents";
return FillBeFile(ifile.Data(), noevents, rtreename, rbranchname, ofile.Data(), wtreename, wbranchname, "");
}
Int_t BeWork::FillBeSimFile(const char *inputfile, const char* outputfile, Long64_t noevents) {
//inputfile:
//noevents:
//outputfile: if empty, character sequence "Exp" in inputfile will be changed to "Be"
TString ifile = fWorkDir + '/' + inputfile;
TString ofile = fWorkDir + '/' + outputfile;
const char* rtreename = "simbe";
const char* rbranchname = "BeSimPhys";
const char* wtreename = "simbeonly";
const char* wbranchname = "BeSimEvents";
return FillBeFile(ifile.Data(), noevents, rtreename, rbranchname, ofile.Data(), wtreename, wbranchname, "");
}
void BeWork::MixSimBeFiles(const Int_t infiles, const char* treename, ...) {
// obsolete function which probably won't be developed more
//
//dodelat podminku "soubor" && vstupy != 0
// if (!i) return 0;
va_list params;
va_start(params, treename);
TString file[infiles];
Int_t events[infiles];
for (Int_t i = 0; i < infiles; i++) {
file[i] = va_arg(params, char*);
events[i] = va_arg(params, Int_t);
}
va_end(params);
//open output mix file
BePureEvent *revent = new BePureEvent();
TFile fw("mixed.root", "RECREATE");
if (!fw.IsOpen()) {
Error("BeWork::MixSimBeFiles", "File %s was not opened", "mixed.root");
return;
}
TTree *tw = new TTree(treename, "tree of mixed simulated events");
//!!!! make possible to mix also BeEvent trees
tw->Bronch("BeSimMix", "BePureEvent", &revent, 1024000, 99);
for (Int_t i = 0; i < infiles; i++) {
// cout << file[i] << endl;
//file opening
TFile fr(file[i].Data());
if (!fr.IsOpen()) {
Warning("BeWork::MixSimBeFiles", "File \"%s\" was not open", file[i].Data());
continue;
}
TTree *tr = (TTree*)fr.Get(treename);
// tr->GetListOfBranches()->Print();
tr->SetBranchAddress("BeSimEvents", &revent);
Long64_t noEntries = tr->GetEntries(/*getCondition*/);
if (events[i] != 0 && noEntries > events[i]) { noEntries = events[i]; }
// tr->GetEntries();
Info("BeWork::MixSimBeFiles", "%d events from file %s will be processed", (Int_t)noEntries, file[i].Data());
for (Long64_t j = 0; j < noEntries; j++) {
BeWork::printProgBar((Int_t)j, (Int_t)noEntries);
tr->GetEntry(j);
tw->Fill();
}//for j
fr.Close();
if (fr.IsOpen()) {
Warning("BeWork::MixSimBeFiles", "File \"%s\" closing error", file[i].Data());
return;
}
}//for i
fw.cd();
tw->Write();
fw.Close();
return;
}
TTree* BeWork::OpenTree(const char* inputfile, const char* treename,
const Color_t color, const char* friendtreename)
{
// open file with saved TTree named "treename"
// and return pointer to this tree
TFile *fr = new TFile(inputfile, "READ");
if (fr->IsOpen() == 0) {
Warning("FillFile", "File %s was not opened and won't be processed", inputfile);
return 0;
}
TTree *tr = (TTree*)fr->Get(treename);
if (!tr) {
Warning("OpenTree", "Tree \"%s\" was not found in file %s", treename, inputfile);
return 0;
}
//add friend
TString fr_name = friendtreename;
if ( !fr_name.IsNull() ) {
tr->AddFriend(fr_name.Data(), fr);
}
tr->SetMarkerStyle(20);
tr->SetLineColor(color);
tr->SetMarkerColor(color);
gROOT->cd();
return tr;
}
TChain* BeWork::OpenChain(const char* inputfile, int first, int last,
const char* treename, const Color_t color, const char* friendtreename, Option_t* option)
{
// open series of files with identical names numbered
// from first to last with saved TTree named
// "treename" and return pointer to chain constisting
// all of this trees
TString opt = option;
opt.ToLower();
TChain *ch = new TChain(treename);
Char_t fileName[100]; //potreba delsi pole ?232?
// if (opt.Contains("gp")) {
for (Int_t i = first; i <= last; i++) {
sprintf(fileName, "%s%03d.root", inputfile, i);
if (opt.Contains("verbose")) {
cout << fileName << endl;
}
ch->Add(fileName);
}
//adding friend chain
TString fr_name = friendtreename;
if ( !fr_name.IsNull() ) {
TChain *fr_ch = new TChain(fr_name.Data());
TString filename;
for (Int_t i = first; i <= last; i++) {
filename.Form("%s%03d.root", inputfile, i);
// sprintf(fileName, "%s%03d.root", inputfile, i);
// cout << fileName << endl;
if (opt.Contains("verbose")) {
cout << filename << endl;
}
// ch->Add(fileName);
fr_ch->Add(filename.Data());
}
// ch->AddFriend(fr_name);
ch->AddFriend(fr_ch);
}//if
ch->SetMarkerStyle(20);
ch->SetLineColor(color);
ch->SetMarkerColor(color);
gROOT->cd();
return ch;
}
TGeoManager* BeWork::BuildGeometry()
{
//geometry units are cm
//dodelat nejakou kontrolu, jestli uz nahodou takova geometrie neexistuje
TGeoManager *geometry = new TGeoManager("SETUP", "6Be structure experimental setup");
//materials
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
// static TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085, 14, 2.332);
//media
TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
// static TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);
TGeoVolume *top = geometry->MakeBox("Top", medVacuum, 100., 100., 100.);
geometry->SetTopVolume(top);
top->SetLineColor(kMagenta);
//target
TGeoVolume *target = MakeTarget();
target->SetName("target");
top->AddNode(target, 1);
//first telescope
const Double_t z11 = fT1SimPosition; //in cm
const Double_t z12 = z11 + 1.; //10.1; //in cm //original value
// const Double_t z12 = z11 + 0.65; //10.1; //in cm
// const Double_t z13 = z11 + 3.5; //12.6; //in cm //original value
const Double_t z13 = z11 + 2.9; //12.6; //in cm
TGeoVolume *annDSD1 = MakeAnnularDetector(64, 32, fX11, fX11_FD, fX11_BD); //backdl (from previous detector) + frontdl = 1.77
annDSD1->SetName("T11");
top->AddNode(annDSD1, 1, new TGeoTranslation(0, 0, z11));
TGeoVolume *annSSD1 = MakeAnnularDetector(64, 1, fX12, fX12_FD, fX12_BD); //backdl (from previous detector) + frontdl = 1.7 + 2.36 = 4.06 == XD12
annSSD1->SetName("T12");
top->AddNode(annSSD1, 1, new TGeoTranslation(0, 0, z12));
TGeoVolume *annCsI1 = MakeCsIDetector(); //backdl (from previous detector) + frontdl = 2. + ?? = ??? ==
// TGeoVolume *annCsI1 = MakeCsIDetectorMS("T13");
annCsI1->SetName("T13");
top->AddNode(annCsI1, 1, new TGeoTranslation(0., 0., z13));
//*/
//second telescope
const Double_t z21 = fT2SimPosition; //in cm
const Double_t z22 = z21 + 1.; //31.1; //in cm //original value
// const Double_t z22 = z21 + 0.65; //31.1; //in cm
// const Double_t z23 = z21 + 3.5; //33.6; //in cm //original value
const Double_t z23 = z21 + 2.9; //33.6; //in cm
TGeoVolume *annDSD2 = MakeAnnularDetector(64, 32, fX21, fX21_FD, fX21_BD); //!!!!!! thickness of T21 was changed !!!!!!
annDSD2->SetName("T21");
top->AddNode(annDSD2, 1, new TGeoTranslation(0, 0, z21));
TGeoVolume *annSSD2 = MakeAnnularDetector(64, 1, fX22, fX22_FD, fX22_BD);
annSSD2->SetName("T22");
top->AddNode(annSSD2, 1, new TGeoTranslation(0, 0, z22));
TGeoVolume *annCsI2 = MakeCsIDetector();
// TGeoVolume *annCsI2 = MakeCsIDetectorMS("T23");
annCsI2->SetName("T23");
top->AddNode(annCsI2, 1, new TGeoTranslation(0., 0., z23));
geometry->CloseGeometry();
Info("BeWork::BuildGeometry", "Distances in virtual setup: zT11 = %4.1f cm zT21 = %4.1f cm", fT1SimPosition, fT2SimPosition);
return geometry;
}
TGeoVolume* BeWork::MakeAnnularDetector(Int_t nosecs, Int_t norings,
Double_t thickness, Double_t fdeadlayer, Double_t bdeadlayer)
{
//nosecs: number of sectors
//norings: number of rings
//thickness: detector thickness in microns
//fdeadlayer: frontside deadlayer thickness in microns
//bdeadlayer: backside deadlayer thickness in microns
//annular detector
const Double_t siThickness = thickness*MicToCm();
const Double_t frontdl = fdeadlayer*MicToCm();
const Double_t backdl = bdeadlayer*MicToCm();
const Double_t sSiRmin = 1.6; //r_min of sensitive region in cm
const Double_t sSiRmax = 4.2; //r_max of sensitive region in cm
const Double_t siRmin = 1.4; //r_min of Si plate in cm
// const Double_t siRmin = 1.6; //r_min of Si plate in cm, version of MS
Double_t deskThickness = 0.38;
if (siThickness > deskThickness) { deskThickness = siThickness; }
// const Double_t dist = 0.01; //in cm; distance between two neighbourgh rings or sectors //!!!!!!!!!! original value
const Double_t dist = 0.00001; //in cm; distance between two neighbourgh rings or sectors //!!!!!!!!!! almost zero
// const Double_t dist = 1.; //distance between two neighbourgh rings or sectors
//materials
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
/*static*/ TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085, 14, 2.332);
//media
TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
/*static*/ TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);
// cout << endl << endl;
// printf("\t%3.3f\t%3.3f\n", deskThickness, siThickness);
// cout << endl << endl;
//detector box
// TGeoBBox *detbox = new TGeoBBox(sSiRmax + 1., sSiRmax + 1., deskThickness);
TGeoBBox *detbox = new TGeoBBox(sSiRmax + 1., sSiRmax + 1., siThickness/2.+0.1);
TGeoVolume *annSSD = new TGeoVolume("annSSD", detbox, medVacuum);
TGeoTube *senslayer = new TGeoTube(siRmin, sSiRmax, siThickness/2.);
TGeoVolume *siDisc = new TGeoVolume("discSi", senslayer, medSi);
siDisc->SetLineColor(3);
annSSD->AddNode(siDisc, 1);
/////////////
////rings
/////////////
TGeoTube *ringtube = new TGeoTube(sSiRmin, sSiRmax, backdl/2.);
TGeoVolume *rlayer = new TGeoVolume("rlayer", ringtube, medVacuum);
Double_t rmin = 0, rmax = 0;
for (Int_t i = 0; i < norings; i++) {
rmin = sSiRmin + i*(sSiRmax-sSiRmin)/norings;
rmax = sSiRmin + (i+1)*(sSiRmax-sSiRmin)/norings - dist;
TGeoTube *shapering = new TGeoTube(rmin, rmax, backdl/2.);
TGeoVolume *ring = new TGeoVolume("ring", shapering, medSi);
ring->SetLineColor(2);
rlayer->AddNode(ring, i);
}
annSSD->AddNode(rlayer, 1, new TGeoTranslation(0, 0, siThickness/2.+backdl/2.));
////////////////
//// sectors
////////////////
const Double_t tubsphi = 180./nosecs;
// cout << tubsphi << endl;
TGeoTube *sectube = new TGeoTube(sSiRmin, sSiRmax, frontdl/2.);
TGeoVolume *slayer = new TGeoVolume("slyer", sectube, medVacuum);
//sensitive strip
TGeoTubeSeg *tubesec = new TGeoTubeSeg(sSiRmin, sSiRmax, frontdl/2., -tubsphi, tubsphi);
tubesec->SetName("tubesec");
TGeoBBox *space = new TGeoBBox(sSiRmax, dist/2., frontdl);
space->SetName("space");
TGeoUnion *sub1 = new TGeoUnion(space, space, new TGeoRotation("r1", tubsphi, 0, 0), new TGeoRotation("r2", -tubsphi, 0, 0));
TGeoCompositeShape *cs1 = new TGeoCompositeShape("cs1", sub1);
TGeoSubtraction *sub2 = new TGeoSubtraction(tubesec, cs1, 0, 0);
TGeoCompositeShape *ssec = new TGeoCompositeShape("ssec", sub2);
TGeoVolume *senssec = new TGeoVolume("senssec", ssec, medSi);
senssec->SetLineColor(2);
for (Int_t i = 0; i < nosecs; i++) {
slayer->AddNode(senssec, i, new TGeoRotation("secrot", 2*i*tubsphi+tubsphi, 0, 0));
}
TGeoRotation platerot("platerot", 90, 180, 0);
TGeoTranslation platetrans("platetrans", 0, 0, -(siThickness/2.+frontdl/2.));
annSSD->AddNode(slayer, 1, new TGeoCombiTrans(platetrans, platerot) );
return annSSD;
}
TGeoVolume* BeWork::MakeCsIDetector()
{
//all sizes are in cm
//one CsI sector parameters
const Double_t segX1 = 1.8; //size of outer edge
const Double_t segX2 = 0.6; //size of inner edge
const Double_t segY = 1.9; //thickness
const Double_t segZ = 3.0; //length from the center to the edge
const Double_t rIn = 3.75; //diameter of the holel; radius = 1.875 cm
const Double_t phisec = 180./16.;
// const Double_t fdl = 12.;
//materials
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085, 14, 2.332);
//media
TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);
// cout << endl << endl;
// printf("\t%3.3f\t%3.3f\t%3.3f\n", rIn + segZ + 1., segY/2.+0.2, rIn/2. + segZ/2.);
// cout << endl << endl;
TGeoBBox *detbox = new TGeoBBox(rIn + segZ + 1., rIn + segZ + 1., segY/2.+0.2);
TGeoVolume *annDet = new TGeoVolume("annDet", detbox, medVacuum);
TGeoTrd1 *csI = new TGeoTrd1(segX1/2., segX2/2., segY/2., segZ/2.);
TGeoVolume *csSec = new TGeoVolume("csSec", csI, medSi);
csSec->SetLineColor(kGray+1);
Double_t phi = 0;
for (Int_t i = 0; i < 16; i++) {
phi = 2*i*phisec+phisec;
TGeoRotation platerot("platerot", -phi, 90., 0.);
TGeoTranslation platetrans("platetrans", (rIn/2. + segZ/2.)*TMath::Sin(phi*TMath::DegToRad()), (rIn/2. + segZ/2.)*TMath::Cos(phi*TMath::DegToRad()), 0);
annDet->AddNode(csSec, i, new TGeoCombiTrans(platetrans, platerot));
}
return annDet;
}
TGeoVolume* BeWork::MakeCsIDetectorMS(TString name)
{
//all sizes are in cm
const Int_t nosecs = 16;
const Double_t siThickness = 1.9;
// const Double_t frontdl = 4.;
// const Double_t backdl = 4.;
const Double_t sSiRmin = 1.6; //r_min of sensitive region in cm
const Double_t sSiRmax = 4.2; //r_max of sensitive region in cm
// const Double_t siRmin = 1.4; //r_min of Si plate in cm
Double_t deskThickness = 0.38;
if (siThickness > deskThickness) { deskThickness = siThickness; }
// const Double_t dist = 0.01; //in cm; distance between two neighbourgh rings or sectors //!!!!!!!!!! original value
const Double_t dist = 0.0001; //in cm; distance between two neighbourgh rings or sectors //!!!!!!!!!! almost zero
// const Double_t dist = 1.; //distance between two neighbourgh rings or sectors
//materials
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085, 14, 2.332);
//media
TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);
TGeoBBox *detbox = new TGeoBBox(sSiRmax + 1., sSiRmax + 1., siThickness/2.+0.1);
TGeoVolume *annDet = new TGeoVolume("annDet", detbox, medVacuum);
////////////////
//// sectors
////////////////
const Double_t tubsphi = 180./nosecs;
// cout << tubsphi << endl;
TGeoTube *sectube = new TGeoTube(sSiRmin, sSiRmax, siThickness/2.);
TGeoVolume *slayer = new TGeoVolume(name.Data(), sectube, medVacuum);
//sensitive strip
TGeoTubeSeg *tubesec = new TGeoTubeSeg(sSiRmin, sSiRmax, siThickness/2., -tubsphi, tubsphi);
tubesec->SetName("tubesec");
TGeoBBox *space = new TGeoBBox(sSiRmax, dist/2., siThickness);
space->SetName("space");
TGeoUnion *sub1 = new TGeoUnion(space, space, new TGeoRotation("r1", tubsphi, 0, 0), new TGeoRotation("r2", -tubsphi, 0, 0));
TGeoCompositeShape *cs1 = new TGeoCompositeShape("cs1", sub1);
TGeoSubtraction *sub2 = new TGeoSubtraction(tubesec, cs1, 0, 0);
TGeoCompositeShape *ssec = new TGeoCompositeShape("ssec", sub2);
TGeoVolume *csSec = new TGeoVolume("csSec", ssec, medSi);
csSec->SetLineColor(kGray+1);
for (Int_t i = 0; i < nosecs; i++) {
slayer->AddNode(csSec, i, new TGeoRotation("secrot", 2*i*tubsphi+tubsphi, 0, 0));
}
TGeoRotation platerot("platerot", 90, 180, 0);
TGeoTranslation platetrans("platetrans", 0., 0., 0.);
annDet->AddNode(slayer, 1, new TGeoCombiTrans(platetrans, platerot) );
// TGeoTrd1 *csI = new TGeoTrd1(segX1/2., segX2/2., segY/2., segZ/2.);
// TGeoVolume *csSec = new TGeoVolume("csSec", csI, medSi);
// csSec->SetLineColor(kGray+1);
// Double_t phi = 0;
// for (Int_t i = 0; i < 16; i++) {
// phi = 2*i*phisec+phisec;
// TGeoRotation platerot("platerot", -phi, 90., 0.);
// TGeoTranslation platetrans("platetrans", (rIn/2. + segZ/2.)*TMath::Sin(phi*TMath::DegToRad()), (rIn/2. + segZ/2.)*TMath::Cos(phi*TMath::DegToRad()), 0);
// annDet->AddNode(csSec, i, new TGeoCombiTrans(platetrans, platerot));
// }
return annDet;
}
TGeoVolume* BeWork::MakeTarget()
{
//one CsI sector parameters
const Double_t rmin = 0.;
const Double_t rmax = 2.;
const Double_t dz = (T_THICKNESS+T_WIN_THICK*2)/10000; //from mcm to cm
const Double_t d_mat = T_THICKNESS/10000; //from mcm to cm
//materials
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
// TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085, 14, 2.332);
//media
TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
// TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);
//stainless steel box for target material
TGeoTube *tartube = new TGeoTube(rmin, rmax, dz/2.);
TGeoVolume *target = new TGeoVolume("target", tartube, medVacuum);
target->SetLineColor(kGray);
//target material
TGeoTube *tar_mat_tube = new TGeoTube(rmin, rmax, d_mat/2.);
TGeoVolume *tar_mat = new TGeoVolume("tarmaterial", tar_mat_tube, medVacuum);
tar_mat->SetLineColor(kBlue);
target->AddNode(tar_mat, 1);
return target;
}
void BeWork::FillSimFile(const char* outputfile, const Int_t noevents, const char* generator, Option_t *opt,
Double_t beThetaMin, Double_t beThetaMax,
UInt_t gseed/*, const char* cutfile*/)
{
//TTree "simbe" will be filled and saved into outputfile, tree "simbe" contains two branches, "BeSimPhys" with virtually detected particles coming
// out of investigated reaction and the branch "BeSimKin" with whole kinematic information
//outputfile:
//noevents: number of simulated events, if 0 then all events from generator
//generator: name of file with generated impulses
//gseed: seed for beam position random generator
//opt: root file option
//beamX: x coordinate position of the center of beam in cm
//beamY: y coordinate position of the center of beam in cm
//beamXfwhm: x width of beam in cm
//beamYfwhm: y width of beam in cm
//beThetaMin: minimum thetaBe in CM in
//beThetaMax: maximum thetaBe in CM in
//beamenergy: beam energy in MeV
//beamtheta: beam inclination in rad
//beamphi: beam inclination in rad
TString ofile = fWorkDir + '/' + outputfile;
// cout << ofile.Data() << endl;
TGeoManager *geom = BuildGeometry();
BeEvent *simevent = new BeEvent(fConfigFile); //fixme initialize fCutsFile
//output file containing tree of two branches
TFile *fw = new TFile(ofile.Data(), opt);
if (!fw) {
Error("BeWork::FillSimFile", "File %s was not created", ofile.Data());
return;
}
TTree *tw = new TTree("simbe", "tree of simulated events");
tw->Bronch("BeSimPhys", "BeEvent", &simevent, 1024000, 99); //velky buffer !!!!!!!!!!!!!!!velky problem s vetvi, prilis komprimovana informace
//beam position information in Friend Tree
TTree *tbeam = new TTree("sbeam", "tree with beam position");
//uzavreno v cyklu
Int_t noEntries = 0;
TString gen = generator;
if (gen.Length() != 0) {
noEntries = CountLines(generator);
if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
Info("BeWork::FillSimFile", "File \"%s\" will be used as generator", generator);
Info("BeWork::FillSimFile", "In file %s %d events will be processed", fw->GetName(), noEntries);
}//if
else {
Info("BeWork::FillSimFile", "In file %s %d events will be processed", fw->GetName(), noevents);
}
Info("BeWork::FillSimFile", "Output file %s will be filled", fw->GetName());
MonteCarloState(noevents, generator, tw, simevent, geom,
tbeam, beThetaMin, beThetaMax, gseed);
// dodelat dalsi dva generatory ------ jake?
fw->cd();
tw->Write();
tbeam->Write();
fw->Close();
delete geom;
delete simevent;
return;
}
void BeWork::MonteCarloState(const Int_t noevents, const char* generator, TTree *writetree, BeEvent *besimevent, TGeoManager *geometry,
TTree *beampos, Double_t reactionAngleMin, Double_t reactionAngleMax, UInt_t gseed)
{
//if generator =NULL, default internal generator will be used
//promenne: simevent, reaction
// energy losses
//beamenergy: in MeV
//beamtheta: in rad
//beamphi: in rad
//reactionAngleMin: in rad
//reactionAngleMax: in rad
Info("BeWork::MonteCarloState", "Function was called");
TString g(generator);
ifstream gen;
//generator file opening
if (!g.IsNull()) {
gen.open(g.Data());
if (!gen.is_open()) {
Error("BeWork::MonteCarloState", "Physical generator was not opened!");
return;
}//if
}//if
//variables for use with generator
Double_t E, p_alpha[3], p_p1[3], p_p2[3], thetaCM; //E is mass or exc. energy of Be, thetaCM is angle of (p,n) reaction
//uzavreno v cyklu
Int_t noEntries = 0;
if (gen.is_open()) {
noEntries = CountLines(generator);
if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
Info("BeWork::MonteCarloState", "%d events will be processed", noEntries);
}//if
else {
noEntries = noevents;
Info("BeWork::MonteCarloState", "%d events will be processed", noEntries);
}
// Int_t noEntries = CountLines(generator);
// if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
// Info("BeWork::FillSimFile", "In file %s %d events will be processed", fw->GetName(), noEntries);
// Info("BeWork::FillSimFile", "Output file %s will be filled", fw->GetName());
//beam size and position: change to function parameter
TRandom3 ranPosition(gseed);
TRandom3 ranEnergy(gseed+1);
TRandom3 ranThetaCM(gseed+3);
Double_t x = fBeamX_MC;
Double_t y = fBeamY_MC;
Double_t z = 0;
Double_t beamT = 0.;
beampos->Branch("bx", &x, "bx/D");
beampos->Branch("by", &y, "by/D");
beampos->Branch("bz", &z, "bz/D");
beampos->Branch("thetaMC", &thetaCM, "thetaMC/D");
beampos->Branch("E_IM", &E, "E_IM/D");
Double_t stateRatio = 0.;
beampos->Branch("sRatio", &stateRatio, "sRatio/D");
TLorentzVector *beamLab = new TLorentzVector();
TLorentzVector *alphaLab = new TLorentzVector();
TLorentzVector *p1Lab = new TLorentzVector();
TLorentzVector *p2Lab = new TLorentzVector();
beampos->Bronch("beamLab.", "TLorentzVector", &beamLab, 1024000, 99);
beampos->Bronch("sALab.", "TLorentzVector", &alphaLab, 1024000, 99);
beampos->Bronch("sP1Lab.", "TLorentzVector", &p1Lab, 1024000, 99);
beampos->Bronch("sP2Lab.", "TLorentzVector", &p2Lab, 1024000, 99);
TLorentzVector *alphaCM = new TLorentzVector();
TLorentzVector *p1CM = new TLorentzVector();
TLorentzVector *p2CM = new TLorentzVector();
TLorentzVector *beCM = new TLorentzVector();
beampos->Bronch("sACM.", "TLorentzVector", &alphaCM, 1024000, 99);
beampos->Bronch("sP1CM.", "TLorentzVector", &p1CM, 1024000, 99);
beampos->Bronch("sP2CM.", "TLorentzVector", &p2CM, 1024000, 99);
beampos->Bronch("sBeCM.", "TLorentzVector", &beCM, 1024000, 99);
//Jacobi "T" system
const Double_t MTx = BeEvent::ReducedMass(kPMASS, kPMASS);
const Double_t MTy = BeEvent::ReducedMass(kPMASS + kPMASS, kAMASS);
TVector3 kTx;
TVector3 kTy;
Double_t sTpp;
Double_t sTapp;
Double_t sCosThetaTk;
beampos->Branch("sTpp", &sTpp, "sTpp/D");
beampos->Branch("sTapp", &sTapp, "sTapp/D");
beampos->Branch("sCosThetaTk", &sCosThetaTk, "sCosThetaTk/D");
//Jacobi "Y" system
const Double_t MYx = BeEvent::ReducedMass(kPMASS, kAMASS);
const Double_t MYy = BeEvent::ReducedMass(kPMASS + kAMASS, kPMASS);
TVector3 kYx;
TVector3 kYy;
Double_t sTap;
Double_t sTpap;
Double_t sCosThetaYk;
beampos->Branch("sTap", &sTap, "sTap/D");
beampos->Branch("sTpap", &sTpap, "sTpap/D");
beampos->Branch("sCosThetaYk", &sCosThetaYk, "sCosThetaYk/D");
BeReaction reaction;
//detectors resolution
TRandom3 *detres = new TRandom3();
TRandom3 *proton_mix = new TRandom3();
// const Double_t beamResolution =
for (Int_t i = 0; i < noEntries; i++) {
// cout << i << endl;
printProgBar(i, noEntries);
//particle tracking
besimevent->Reset();
reaction.Reset();
//reaction position in the thickness of the target
z = ranPosition.Uniform(-T_THICKNESS/10000./2., T_THICKNESS/10000./2.); //in cm
beamT = ranEnergy.Gaus(fTBeamMC*6, fTBeamResMC);
beamT = fTargetLi->GetE(beamT, z*CmToMic());
beamLab->SetPxPyPzE(0., 0., TMath::Sqrt(BeEvent::PC2(beamT, kLiMASS)),
beamT + kLiMASS);
// printf("%3.2f\t%3.2f\n", beamLab->E(), beamT + kLiMASS);
// printf("%3.2f\t%3.2f\n", beamLab->Pz(), TMath::Sqrt(BeEvent::PC2(beamT, kLiMASS)));
// cout << beamLab->Px() << endl
// << beamLab->Py() << endl;
// cout << beamT << endl;
// Info("BeWork::MonteCarloState", "beam z is %f mic, beam energy is %f MeV", z*CmToMic(), benergy);
// Info("BeWork::MonteCarloState",
// "beam z is %f mic, beam energy is %f MeV\n", -z * CmToMic(),
// fTargetLi->GetE0(beamenergy, -z * CmToMic()) );
//phase volume
if (g.IsNull()) {
// reaction.FillProcess(beamenergy, 0., 0.); //pomala fce
reaction.FillProcess(beamT, 0., 0.); //pomala fce
thetaCM = reaction.fLip.GetThetaCM();
E = reaction.fLip.GetMa();
}
//generator
else {
if (gen.good()) {
gen >> E >> p_alpha[0] >> p_alpha[1] >> p_alpha[2] >>
p_p1[0] >> p_p1[1] >> p_p1[2] >>
p_p2[0] >> p_p2[1] >> p_p2[2] >> stateRatio;
if (E > 20) continue; //to high energies cannot be populated with our low energy beam
//fixme beam energy is different for different z
// reaction.FillProcess(beamenergy, 0., 0., p_alpha, p_p1, p_p2, reactionAngleMin, reactionAngleMax);
if (stateRatio < fsRatioMin || stateRatio > fsRatioMax) continue; //to omit unusable events from generator
alphaCM->SetPxPyPzE(p_alpha[0], p_alpha[1], p_alpha[2],
BeEvent::T(
Power(p_alpha[0], 2) + Power(p_alpha[1], 2)
+ Power(p_alpha[2], 2), kAMASS) + kAMASS);
p1CM->SetPxPyPzE(p_p1[0], p_p1[1], p_p1[2],
BeEvent::T(
Power(p_p1[0], 2) + Power(p_p1[1], 2)
+ Power(p_p1[2], 2), kPMASS) + kPMASS);
p2CM->SetPxPyPzE(p_p2[0], p_p2[1], p_p2[2],
BeEvent::T(
Power(p_p2[0], 2) + Power(p_p2[1], 2)
+ Power(p_p2[2], 2), kPMASS) + kPMASS);
*beCM = *p1CM + *p2CM + *alphaCM;
//Jacobi "T" system
kTx = (p1CM->Vect() - p2CM->Vect()) * 0.5;
kTy =
(4. * (p1CM->Vect() + p2CM->Vect())
- 2. * alphaCM->Vect()) * (1. / 6.);
sTpp = kTx.Mag2() / (2 * MTx);
sTapp = kTy.Mag2() / (2 * MTy);
sCosThetaTk = kTx.Dot(kTy) / (kTx.Mag() * kTy.Mag());
//Jacobi "Y" system
kYx = (4. * p1CM->Vect() - alphaCM->Vect()) * (1./5.);
kYy = ((p1CM->Vect() + alphaCM->Vect())
- 5. * p2CM->Vect()) * (1. / 6.);
sTap = kYx.Mag2() / (2 * MYx);
sTpap = kYy.Mag2() / (2 * MYy);
sCosThetaYk = kYx.Dot(kYy) / (kYx.Mag() * kYy.Mag());
thetaCM = ranThetaCM.Uniform(reactionAngleMin, reactionAngleMax);
// reaction.FillProcess(beamT, 0., 0., p_alpha, p_p1, p_p2,
// reactionAngleMin, reactionAngleMax);
reaction.FillProcess(beamT, 0., 0., p_alpha, p_p1, p_p2,
thetaCM);
if (!gen.good()) {
Warning("BeWork::FillSimKinFile", "End of file %s was reached at %d event ", generator, i);
break;
}//if
}//if
}//else
TLorentzVector vAlpha(reaction.GetAlpha());
TLorentzVector vP1;
TLorentzVector vP2;
if (proton_mix->Uniform() < 0.5) {
vP1 = reaction.GetP1();
vP2 = reaction.GetP2();
}
else {
vP1 = reaction.GetP2();
vP2 = reaction.GetP1();
}
// alphaT = vAlpha.E() - kAMASS;
// p1T = vP1.E() - kPMASS;
// p2T = vP2.E() - kPMASS;
//beam position
//beam size
//Gauss
x = ranPosition.Gaus(fBeamX_MC, fBeamX_sigma_MC);
y = ranPosition.Gaus(fBeamY_MC, fBeamY_sigma_MC);
//rectangle (square)
// x = ranPosition.Uniform(-beamXfwhm, beamXfwhm); //Gaus(beamX, beamXfwhm/2.35);
// y = ranPosition.Uniform(-beamYfwhm, beamYfwhm);
//reaction position in the thickness of the target
// z = ranPosition.Uniform(-T_THICKNESS/10000./2., T_THICKNESS/10000./2.); //in cm
// Info("BeWork::MonteCarloState", "beam z is %f", z*CmToMic());
//elipse (circle) is not possible
//should be treated on level of TCut
// Info("BeWork::MonteCarloState", "Input energy of alpha is %f", alphaT);
*alphaLab = vAlpha;
*p1Lab = vP1;
*p2Lab = vP2;
ParticleTracking(&vAlpha, 24,
geometry, besimevent, detres, x, y, z); //funguje pomerne rychle
ParticleTracking(&vP1, 11,
geometry, besimevent, detres, x, y, z);
ParticleTracking(&vP2, 11,
geometry, besimevent, detres, x, y, z);
//work with energies to fill the all BeEvent variables
besimevent->SetDeltaE();
besimevent->SetHits();
// fw->cd();
// writetree->Fill(); // pomala funkce
// beampos->Fill();
if (fBeOnly) {
if (besimevent->WriteCondition()) {
writetree->Fill();
beampos->Fill();
} //zkontrolovat
} else {
writetree->Fill(); // pomala funkce
beampos->Fill();
}
}//for all events
//closing of generator file
if (!g.IsNull()) {
gen.close();
if (gen.is_open()) {
Warning("BeWork::FillSimKinFile", "File %s closing error\n", generator);
}//if
}//if
//delete all TLorentzVectors:
//fixme
delete detres;
delete proton_mix;
delete beamLab;
delete alphaLab;
delete p1Lab;
delete p2Lab;
delete alphaCM;
delete p1CM;
delete p2CM;
return;
}
void BeWork::ParticleTracking(const TLorentzVector *particle, const Int_t pID,
TGeoManager *geom, BeEvent *event, TRandom3 *detres,
Double_t beamX, Double_t beamY, Double_t beamZ)
{
//particle: treated particle
//silosses: object for energy losses calculations in Si detectors
//csilosses: object for energy losses calculations in CsI(Tl) detectors
//geom: geometry to be used
//event: object of event class to be filled
//detres:
//beamX: in cm
//beamY: in cm
Double_t kinE0 = 0;
TELoss *silosses;
TELoss *csilosses;
TELoss *targetlosses;
TELoss *target_win_losses;
if (pID == kPROTON) { //proton
kinE0 = particle->E() - kPMASS;
silosses = fSiP;
csilosses = fCsIP;
targetlosses = fTargetP;
target_win_losses = fTargetWinP;
}
else {
if (pID == kALPHA) { //alpha
kinE0 = particle->E() - kAMASS;
silosses = fSiAlpha;
csilosses = fCsIAlpha;
targetlosses = fTargetAlpha;
target_win_losses = fTargetWinAlpha;
}
else {
Info("BeWork::ParticleTracking",
"I can only track either protons or alphas.");
return;
}
}//else
Double_t kinE = 0;
const Double_t cos_theta = TMath::Cos(particle->Theta());
TVector3 pardir(particle->Vect().Unit()); //particle direction
geom->InitTrack(beamX, beamY, beamZ, pardir.Px(), pardir.Py(), pardir.Pz());
TString anode;
Double_t dist; //in cm
TString detector; //mother node
Int_t ring, sec;
const Int_t kNoStrip = -1;
Double_t deltaE = 0;
const Double_t fdlCsI = fX13_FD; //in mcm
// const Double_t fdlCsI = 1.; //in mcm
while (!geom->IsOutside()) {
if (kinE0 == 0) {
break; }
geom->FindNextBoundary();
ring = kNoStrip;
sec = kNoStrip;
while (geom->GetMother()) { //inside Si detector //XXX problem: overit
//funguje, pouze pokud se castice nezastavi, jinak se zasekne
geom->FindNextBoundary();
anode = geom->FindNode()->GetName();
if (anode.Contains("senssec") ) { //as dead layer
//energy losses in front dead layer of Si detector
sec = geom->FindNode()->GetNumber();
dist = geom->GetStep();
// cout << dist << endl;
kinE0 = silosses->GetE(kinE0, CmToMic()*dist);
}//if
if (anode.Contains("discSi") ) {
//energy losses in sensitive layer of Si detector
dist = geom->GetStep();
kinE = silosses->GetE(kinE0, CmToMic()*dist);
deltaE = kinE0 - kinE;
kinE0 = kinE;
detector = geom->GetMother()->GetName(); //used to identify detector
//for write the energy
}//if
if (anode.Contains("ring") ) { //as dead layer
//energy losses in back dead layer of Si detector
ring = geom->FindNode()->GetNumber();
dist = geom->GetStep();
kinE0 = silosses->GetE(kinE0, CmToMic()*dist);
}//if
if (anode.Contains("csSec") ) { //energy losses in CsI detector
// cout << "CsI" << endl;
sec = geom->FindNode()->GetNumber();
ring = 0;
//dead layer
kinE0 = csilosses->GetE(kinE0, fdlCsI/cos_theta);
//sensitive layer
dist = geom->GetStep();
kinE = csilosses->GetE(kinE0, CmToMic()*dist);
deltaE = kinE0 - kinE;
kinE0 = kinE;
detector = geom->GetMother()->GetName(); //used to identify detector
//for write the energy
// printf("CsI: %d\t%3.2f\t%3.2f\t%s\n", sec, dist, deltaE, detector.Data());
}//if
if (anode.Contains("tarmaterial") ) { //energy losses in target
dist = geom->GetStep();
kinE0 = targetlosses->GetE(kinE0, CmToMic()*dist);
}//if
if (anode.Contains("target") ) { //energy losses in target windows
dist = geom->GetStep();
kinE0 = target_win_losses->GetE(kinE0, CmToMic()*dist);
}//if
geom->Step();
}//while inside detector
Double_t E_write;
//write in the simulated energy deposites
if ( event && (sec != kNoStrip) && (ring != kNoStrip) ) {
//first telescope
if (detector.Contains("T11")) {
E_write = deltaE + detres->Gaus(0., fSiRes); //calculated energy loss + detector resolution
if (E_write >= event->calSi->fC[3][0]) event->fT11.fESec[sec/2] = E_write;
E_write = deltaE + detres->Gaus(0., fSiRes); //calculated energy loss + detector resolution
if (E_write >= event->calSi->fC[5][0]) event->fT11.fERing[ring] = E_write;
}//if T11
if (detector.Contains("T12")) {
E_write = deltaE + detres->Gaus(0., fSiRes); //calculated energy loss + detector resolution
if (E_write >= event->calSi->fC[13][sec/4]) event->fT12.fESec[sec/4] = E_write;
event->fTau[0][sec/4] = 90.;
}//if T12
if (detector.Contains("T13")) {
if (deltaE*fCsIResA>fCsIBestResA) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResA/2.35); //pro alfy //calculated energy loss + detector resolution
else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResA/2.35);
if (E_write >= event->calCsIa->fD[15][sec]) event->fE13aSec[sec] = E_write;
if (deltaE*fCsIResP>fCsIBestResP) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResP/2.35); //calculated energy loss + detector resolution
else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResP/2.35);
if (E_write >= event->calCsIp->fD[15][sec]) event->fE13pSec[sec] = E_write;
}//if T13
//second telescope
if (detector.Contains("T21")) {
E_write = deltaE + detres->Gaus(0., fSiRes); //calculated energy loss + detector resolution
if (E_write >= event->calSi->fC[8][0]) event->fT21.fESec[sec/2] = E_write;
E_write = deltaE + detres->Gaus(0., fSiRes); //calculated energy loss + detector resolution
if (ring < 16 && E_write >= event->calSi->fC[10][0]) event->fT21.fERing[ring] = E_write;
if (ring > 15 && E_write >= event->calSi->fC[11][0]) event->fT21.fERing[ring] = E_write;
}//if T21
if (detector.Contains("T22")) {
E_write = deltaE + detres->Gaus(0., fSiRes); //calculated energy loss + detector resolution
if (E_write >= event->calSi->fC[14][sec/4]) event->fT22.fESec[sec/4] = E_write;
event->fTau[1][sec/4] = 90.;
}//if T22
if (detector.Contains("T23")) {
if (deltaE*fCsIResA>fCsIBestResA) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResA/2.35); //pro alfy //calculated energy loss + detector resolution
else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResA/2.35);
if (E_write >= event->calCsIa->fD[16][sec]) event->fE23aSec[sec] = E_write;
if (deltaE*fCsIResP>fCsIBestResP) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResP/2.35); //calculated energy loss + detector resolution
else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResP/2.35);
if (E_write >= event->calCsIp->fD[16][sec]) event->fE23pSec[sec] = E_write;
}//if T23
}//if writein
geom->Step();
}//while tracking through whole detector system
return;
}
void BeWork::FillSimKinFile(Int_t INPUTS, Double_t beamtheta, Double_t beamphi,
const char* outputfile, Option_t *opt) {
//simulation of kinematic variables using BinaryReaction class;
//make TTree object named "KIN" in file
//INPUTS: number of events
//outputfilename: name of root file
//opt: file access option
BeReaction *sReaction = new BeReaction();
Info("BeWork:FillSimKinFile()",
"MC simulation using phase volume generator will be done");
if (INPUTS == 0) {
cout << "\nAssign number of events \n";
return;
// cin >> INPUTS;
}
cout << "\n\nCalculating...\n\n";
TFile *fw = TFile::Open(outputfile, opt);
if (!fw) {
Error("BeWork::FillSimKinFile", "File %s was not created", outputfile);
return;
}
TTree *tree = new TTree("KIN", "two body kinematics");
tree->Bronch("MCsimple", "BeReaction", &sReaction);
for (Int_t i = 1; i <= INPUTS; i++) {
printProgBar(i, INPUTS);
sReaction->FillProcess(fTBeamMC*6, beamtheta, beamphi);
tree->Fill();
}
cout << endl;
tree->Write();
fw->Close();
if (fw->IsOpen()) {
Warning("BeWork::FillSimKinFile", "File %s closing error", outputfile);
}
delete fw;
delete sReaction;
}
void BeWork::FillSimKinFile(const char* generator, Int_t INPUTS,
/*Double_t beamenergy,*/ Double_t beamtheta, Double_t beamphi,
const char* outputfile, Option_t *opt) {
// function variables: path to generator, INPUTS, outputfile option
Info("BeWork::FillSimKinFile",
"MC simulation using physical generator will be done");
// cout << beamenergy << endl;
// cout << beamtheta << endl;
// cout << beamphi << endl;
ifstream gen;
gen.open(generator);
if (!gen.is_open()) {
Error("BeWork::FillSimKinFile", "Physical generator was not opened!");
return;
} //if
//class which will be the structure of TTree branch filled into
BeReaction *sReaction = new BeReaction();
//checking number of inputs
if (INPUTS == 0) {
Info("FillSimKinFile", "O inputs was filled");
return;
// cout << "\nAssign number of events \n";
// cin >> INPUTS;
}
cout << "\n\nCalculating...\n\n";
//outputfile, tree and branch creation
TFile *fw = TFile::Open(outputfile, opt); //name of the output file
if (!fw) {
Error("BeWork::FillSimKinFile", "File %s was not opened", outputfile);
return;
}
TTree *tree = new TTree("KIN", "some suitable title"); //name of the tree, title of the three
tree->Bronch("MCgen", "BeReaction", &sReaction); //name of the branch
//filling the tree
Double_t E, p_alpha[3], p_p1[3], p_p2[3], thetaCM;
for (Int_t i = 1; i <= INPUTS; i++) {
printProgBar(i, INPUTS); //zrychlit fci
if (gen.good()) {
gen >> E >> p_alpha[0] >> p_alpha[1] >> p_alpha[2] >> p_p1[0]
>> p_p1[1] >> p_p1[2] >> p_p2[0] >> p_p2[1] >> p_p2[2]
>> thetaCM;
// sReaction->FillProcess(beamenergy, beamtheta, beamphi, p_alpha, p_p1, p_p2, thetaCM);
sReaction->FillProcess(fTBeamMC*6, beamtheta, beamphi, p_alpha,
p_p1, p_p2, 0., TMath::Pi());
if (!gen.good()) {
Warning("BeWork::FillSimKinFile", "End of file %s was reached",
generator);
tree->Fill();
break;
} //if
} //if
tree->Fill();
} //for
cout << endl;
tree->Write();
fw->Close();
if (fw->IsOpen()) {
Warning("BeWork::FillSimKinFile", "File %s closing error", outputfile);
}
gen.close();
if (gen.is_open()) {
Warning("BeWork::FillSimKinFile", "File %s closing error\n", generator);
} //if
delete fw;
delete sReaction;
return;
}
void BeWork::printProgBar(Int_t percent)
{
// std::string bar;
TString bar;
for(int i = 0; i < 50; i++){
if( i < (percent/2)){
// bar.replace(i,1,"=");
bar.Replace(i,1,"=");
}
else {
if( i == (percent/2)){
// bar.replace(i,1,">");
bar.Replace(i,1,">");
}
else {
// bar.replace(i,1," ");
bar.Replace(i,1," ");
}
}
}
cout << "\r" "[" << bar << "] ";
cout.width( 3 );
// cout << percent << "% " << flush;
if (percent < 100) cout << percent << "% " << flush;
else cout << percent << "% " << endl;
return;
}
void BeWork::printProgBar(Int_t i, Int_t inputs)
{
// printProgBar(0);
// if ( (i%(inputs/100)) != 0 ) return;
Int_t percent = i*100/inputs;
if ( inputs >= 100 && (i%(inputs/100)) == 0 ) {
printProgBar(percent);
}
if (i == inputs-1) { printProgBar(100); }
return;
}
Int_t BeWork::CountLines(const char* file, Int_t maxlinelength) {
char countline[maxlinelength];
TString fname(file);
ifstream rfile;
//generator file opening
if (!fname.IsNull()) {
rfile.open(fname.Data());
if (!rfile.is_open()) {
Error("BeWork::CountLines", "File \"%s\" was not opened!", fname.Data());
return 0;
}//if
}//if
Int_t j = 0;
while (!rfile.eof()) {
rfile.getline(countline, maxlinelength);
j++;
}
rfile.close();
if (rfile.is_open()) {
Warning("BeWork::CountLines", "File %s closing error\n", fname.Data());
}//if
return j;
}
void BeWork::CreateTELosses() {
//energy looses in silicon
//alpha in Si
fSiAlpha = new TELoss();
fSiAlpha->SetEL(SI_NELEMENTS, SI_RHO);
fSiAlpha->AddEL(SI_1_Z, SI_1_A, SI_1_W);
fSiAlpha->SetZP(2., 4.); //alphas
//100000 je urcite v poradku, 80000 staci, 50000 malo
fSiAlpha->SetEtab(100000, 200.); //deltaE
//proton in Si
fSiP = new TELoss();
fSiP->SetEL(SI_NELEMENTS, SI_RHO);
fSiP->AddEL(SI_1_Z, SI_1_A, SI_1_W);
fSiP->SetZP(1., 1.); //protons
fSiP->SetEtab(100000, 200.); //deltaE //100000 je urcite v poradku, 80000 staci, 50000 malo
//energy looses in CsI(Tl)
//alpha in CsI
fCsIAlpha = new TELoss();
fCsIAlpha->SetEL(CSI_NELEMENTS, CSI_RHO);
fCsIAlpha->AddEL(CSI_1_Z, CSI_1_A, CSI_1_W);
fCsIAlpha->AddEL(CSI_2_Z, CSI_2_A, CSI_2_W);
fCsIAlpha->SetZP(2., 4.); //alphas
//100000 je urcite v poradku, 80000 staci, 50000 malo
fCsIAlpha->SetEtab(100000, 200.); //deltaE
//proton in CsI
fCsIP = new TELoss();
fCsIP->SetEL(CSI_NELEMENTS, CSI_RHO);
fCsIP->AddEL(CSI_1_Z, CSI_1_A, CSI_1_W);
fCsIP->AddEL(CSI_2_Z, CSI_2_A, CSI_2_W);
fCsIP->SetZP(1., 1.); //protons
fCsIP->SetEtab(100000, 200.); //deltaE //100000 je urcite v poradku, 80000 staci, 50000 malo
//energy looses in target
const Double_t targetRho = T_NORMAL_RHO*(273/T_TEMPERATURE)*(T_PRESSURE/1);
//alpha in target
fTargetAlpha = new TELoss();
fTargetAlpha->SetEL(T_NELEMENTS, targetRho);
fTargetAlpha->AddEL(T_1_Z, T_1_A, T_1_W);
fTargetAlpha->SetZP(2., 4.); //alphas
//100000 je urcite v poradku, 80000 staci, 50000 malo
fTargetAlpha->SetEtab(100000, 200.); //deltaE
//proton in target
fTargetP = new TELoss();
fTargetP->SetEL(T_NELEMENTS, targetRho);
fTargetP->AddEL(T_1_Z, T_1_A, T_1_W);
fTargetP->SetZP(1., 1.); //protons
fTargetP->SetEtab(100000, 200.); //deltaE //100000 je urcite v poradku, 80000 staci, 50000 malo
//beam in target
fTargetLi = new TELoss();
fTargetLi->SetEL(T_NELEMENTS, targetRho);
fTargetLi->AddEL(T_1_Z, T_1_A, T_1_W);
fTargetLi->SetZP(3., 6.); //protons
fTargetLi->SetEtab(100000, 200.); //deltaE //100000 je urcite v poradku, 80000 staci, 50000 malo
//alpha in target window
fTargetWinAlpha = new TELoss();
fTargetWinAlpha->SetEL(T_WIN_NELEMENTS, T_WIN_RHO);
fTargetWinAlpha->AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
fTargetWinAlpha->AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
fTargetWinAlpha->AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
fTargetWinAlpha->SetZP(2., 4.); //alphas
//100000 je urcite v poradku, 80000 staci, 50000 malo
fTargetWinAlpha->SetEtab(100000, 200.); //deltaE
//proton in target window
fTargetWinP = new TELoss();
fTargetWinP->SetEL(T_NELEMENTS, T_WIN_RHO);
fTargetWinP->AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
fTargetWinP->AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
fTargetWinP->AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
fTargetWinP->SetZP(1., 1.); //protons
fTargetWinP->SetEtab(100000, 200.); //deltaE //100000 je urcite v poradku, 80000 staci, 50000 malo
//beam in target window
fTargetWinLi = new TELoss();
fTargetWinLi->SetEL(T_NELEMENTS, T_WIN_RHO);
fTargetWinLi->AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
fTargetWinLi->AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
fTargetWinLi->AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
fTargetWinLi->SetZP(3., 6.); //protons
fTargetWinLi->SetEtab(100000, 200.); //deltaE //100000 je urcite v poradku, 80000 staci, 50000 malo
Info("BeWork::CreateTELosses", "TELoss objects were initialized.");
}
/*
* BeWork.h
*
* Created on: 2.2.2010
* Author: Vratislav
*/
//#ifndef BEWORK_H_
//#define BEWORK_H_
#pragma once
#include <TObject.h>
#include <TString.h>
#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoBBox.h>
#include <TGeoTube.h>
#include <TGeoTrd1.h>
#include <TGeoBoolNode.h>
#include <TGeoCompositeShape.h>
#include <TLorentzVector.h>
#include <TUnixSystem.h>
//#include <TMemFile.h>
#include "BeEvent.h"
#include "BeReaction.h"
#include "BePureEvent.h"
#include "../AculData/AculCalibration.h"
#include "../AculData/ConfigDictionary.h"
#include "../TELoss/TELoss.h"
#ifndef __CINT__
#endif /* __CINT __ */
#include <stdarg.h>
using std::flush;
class BeWork {
private:
//public:
// TELosses:
TELoss *fSiAlpha; //!
TELoss *fSiP; //!
TELoss *fCsIAlpha; //!
TELoss *fCsIP; //!
TELoss *fTargetAlpha; //!
TELoss *fTargetP; //!
TELoss *fTargetLi; //!
TELoss *fTargetWinAlpha; //!
TELoss *fTargetWinP; //!
TELoss *fTargetWinLi; //!
//configuration and parameter files
TString fConfigFile; //!
TString fWorkDir; //!
TString fRawFilePath; //!
// TString fCutFile; //!
Bool_t fBeOnly;
//parameter used for choice of generator
Double_t fsRatioMin; //!
Double_t fsRatioMax; //!
//parameters
//detectors resolution used in simulation
Double_t fSiRes; //! in MeV, sigma in Si detectors
Double_t fCsIResP; //! in %, FWHM in CsI for protons
Double_t fCsIBestResP; //! in MeV, the best absolute resolution in CsI for protons
Double_t fCsIResA; //! in %, FWHM in CsI for alphas
Double_t fCsIBestResA; //! in MeV, the best absolute resolution in CsI for alphas
//beam characteristics, used in MC
Double_t fTBeamMC; //! in AMeV //where is it? in the centre? on the target window?
Double_t fTBeamResMC; //! in MeV, sigma for gaus
Double_t fBeamX_MC; //! in cm, x position of beam at target
Double_t fBeamY_MC; //! in cm, x position of beam at target
Double_t fBeamX_sigma_MC; //! in cm; x position sigma of beam at target
Double_t fBeamY_sigma_MC; //! in cm; x position sigma of beam at target
Double_t fT1SimPosition; //! position of the first telescope in cm; used in MC
Double_t fT2SimPosition; //! position of the second telescope in cm; used in MC
//detector thicknesses in mcm
Double_t fX11_FD; //!
Double_t fX11; //!
Double_t fX11_BD; //!
Double_t fX12_FD; //!
Double_t fX12; //!
Double_t fX12_BD; //!
//#define XD13F 14. //fixme original value probably in Si equivalent
Double_t fX13_FD; //!
Double_t fX21_FD; //!
Double_t fX21; //!
Double_t fX21_BD; //!
Double_t fX22_FD; //!
Double_t fX22; //!
Double_t fX22_BD; //!
//#define XD23F 14. //fixme original value probably in Si equivalent
Double_t fX23_FD; //!
private:
// static TRandom3 ranPosition;
static void printProgBar(Int_t percent);
static void printProgBar(Int_t i, Int_t inputs);
void CreateTELosses();
void SetWorkDir();
public:
BeWork();
BeWork(const char* configfile);
BeWork(BeWork &);
virtual ~BeWork();
ClassDef(BeWork, 1);
void ReadConfigFile();
const char* GetWorkDir();
//experimental data analysis
// Int_t FillExpRun();
Int_t FillExpFile(const char* inputrawfile, const char* outputfile,
Long64_t noevents = 0, Option_t *opt = ""); //experimental data processing
Int_t FillBeExpFile(const char *inputfile, const char* outputfile, Long64_t noevents = 0);
//simulations
void FillSimFile(const char* outputfile, const Int_t noevents,
const char* generator = "", Option_t *opt = "CREATE",
Double_t beThetaMin = 0., Double_t beThetaMax = TMath::Pi(),
UInt_t gseed = 0/*, const char* cutfile = ""*/); //simulation of measured variables using BinaryReaction class
Int_t FillBeSimFile(const char *inputfile, const char* outputfile, Long64_t noevents = 0);
static void MixSimBeFiles(const Int_t infiles, const char* treename, ...);
TGeoManager* BuildGeometry(); //using cm size unit
static TGeoVolume* MakeAnnularDetector(Int_t nosec, Int_t norings,
Double_t siThickness, Double_t frontdl, Double_t backdl);
static TGeoVolume* MakeCsIDetector();
static TGeoVolume* MakeCsIDetectorMS(TString name);
static TGeoVolume* MakeTarget();
void ParticleTracking(const TLorentzVector *particle, const Int_t pID,
TGeoManager *geom, BeEvent *event = 0, TRandom3 *detres = 0,
Double_t beamX = 0., Double_t beamY = 0., Double_t beamZ = 0.);
void MonteCarloState(const Int_t noevents, const char* generator,
TTree *writetree, BeEvent *besimevent, TGeoManager *geometry,
TTree* beampos, Double_t reactionAngleMin,
Double_t reactionAngleMax, UInt_t gseed);
//could be one function, just using the switch between the sources of simulated events
//fill the file with kinematical information only by simulated events
void FillSimKinFile(Int_t INPUTS = 0,
Double_t beamtheta = 0., Double_t beamphi = 0.,
const char* outputfile = "collision.root",
Option_t *opt = "CREATE"); //simulation of kinematic variables using BinaryReaction class; opt: file access option
void FillSimKinFile(const char* generator, Int_t INPUTS = 0,
Double_t beamtheta = 0., Double_t beamphi = 0.,
const char* outputfile = "collision_g.root",
Option_t *opt = "CREATE"); //simulation of kinematic variables using G. generator; opt: file access option
//general functions
Int_t FillBeFile(const char* inputfile, Long64_t noevents,
const char* rtree, const char* rbranch, const char* outputfile,
const char* wtree, const char* wbranch, Option_t *opt = ""); //file with Be event only data processing
static TTree* OpenTree(const char* inputfile, const char* treename,
const Color_t col = 1, const char* friendtreename = "");
static TChain* OpenChain(const char* inputfile, int first, int last,
const char* treename, const Color_t color = 1, const char* friendtreename = "", Option_t* option = "");
static Double_t CmToMic() { return 10000.; };
static Double_t MicToCm() { return 0.0001; };
static Int_t CountLines(const char* file, Int_t maxlinelength = 1000);
};
//#endif /* BEWORK_H_ */
/*
* BinaryReaction.cpp
*
* Created on: 24.5.2010
* Author: Vratislav
*/
#include "BinaryReaction.h"
BinaryReaction::BinaryReaction() {
}
BinaryReaction::~BinaryReaction() {
}
ClassImp(BinaryReaction);
//TRandom3 *ranpokus = new TRandom3(1);
TRandom3 BinaryReaction::ranPhi(1);
Int_t BinaryReaction::FillReaction(Double_t _fThetaIn, Double_t _fT1, Double_t _fM[], Double_t _fThetaCM, Double_t _fPhiIn)
{
//_fThetaIn in rad
//_fT1 in MeV
//_fM[] in MeV
//_fThetaCM in rad
//_fPhiIn in rad
fThetaIn = _fThetaIn;
fT1 = _fT1;
for (Int_t i = 0; i < 4; i++) {
fM[i] = _fM[i];
}
fThetaCM = _fThetaCM;
fPhiIn = _fPhiIn;
if (FlatCalculate() != 1) { return -1; }
SetPhi();
//SetOrientation();
SetOrientationV();
return 1;
}
Int_t BinaryReaction::FlatCalculate()
{
const Double_t p1 = Sqrt(Power(fT1, 2) + 2*fT1*fM[0]); //in MeV
const Double_t E = fT1 + fM[0] + fM[1]; //in MeV
//beta and gamma for Lor. trans.
const Double_t beta = p1/E;
const Double_t gamma = 1/Sqrt(1 - Power(beta, 2));
//calculating of T3cm and T4cm
//T3cm = ( Tcm^2 + 2*Tcm*m4*c^2 )/(2*( Tcm + m4*c^2 + m3*c^2 )), where Tcm = T3cm + T4cm
//Ecm = Tcm + m3*c^2 + m4*c^2
//Ecm^2 = E^2 - (p1*c)^2, where E = T1 + m1*c^2 + m2*c^2 and (p1*c)^2 = T1^2 + 2T1*m1*c^2
Double_t Ecm = 0.; //energy in CMS, MeV
Double_t Tcm = 0.; //summary kinetic energy in CMS, MeV
Double_t TaCM = 0., TbCM = 0.;
Double_t pcm = 0.; //impuls in CMS, MeV
Ecm = Sqrt( Power(E, 2) - Power(p1, 2) );
// Tcm = Sqrt( Power(fM[0]+fM[1], 2) + 2*fT1*fM[1] ) - fM[2] - fM[3];
// cout << Tcm << endl;
Tcm = Ecm - fM[2] - fM[3];
// cout << Tcm << endl;
static Int_t counter = 0;
// cout << Tcm << endl;
if (Tcm < 0) {
counter++;
//cout << "\nThere is no solution, Tcm is smaller than threshold energy\n";
Warning("FlatCalculate", "%s: No solution, Tcm is smaller than threshold E", GetName());
cout << "Tcm is " << Tcm << endl
<< "Ecm is " << Ecm << endl
<< "fTa is " << fTa << endl
<< "fTb is " << fTb << endl
<< "fM[0] is " << fM[0] << endl
<< "m_dp is " << 2*fM[2] << endl
<< "fT1 is " << fT1 << endl
<< "E is " << E << endl;
printf("\n\n!!!!!!!!!\t%d\t!!!!!!!!!!!\n\n", counter);
return -1;
}
TaCM = ( Power(Tcm, 2) + 2*Tcm*fM[3] )/( 2*(Tcm + fM[3] + fM[2]) );
TbCM = ( Power(Tcm, 2) + 2*Tcm*fM[2] )/( 2*(Tcm + fM[3] + fM[2]) );
//p3cm = -p4cm ==> |p3cm| = |p4cm| = pcm, (pcm*c)^2 = T3cm^2 + 2*T3cm*m3*c^2
pcm = Sqrt( Power(TaCM, 2) + 2*TaCM*fM[2] );
//Lorentz transformation of impuls from CMS to Lab:
Double_t pax, pay, pbx, pby; //in MeV
Double_t EaCM, EbCM; //in MeV
Double_t paq = 0., pbq = 0.; //in MeV
//evaluate of T3
EaCM = TaCM + fM[2];
pay = pcm*Sin(fThetaCM);
pax = gamma*(pcm*Cos(fThetaCM) + beta*EaCM);
paq = Power(pax, 2) + Power(pay, 2);
fTa = Sqrt( Power(fM[2], 2) + paq ) - fM[2];
//evaluate of T4
EbCM = TbCM + fM[3];
pby = pcm*Sin(fThetaCM + TMath::Pi());
pbx = gamma*(pcm*Cos(fThetaCM + TMath::Pi()) + beta*EbCM);
pbq = Power(pbx, 2) + Power(pby, 2);
fTb = Sqrt( Power(fM[3], 2) + pbq ) - fM[3];
//calculating of theta3 and theta4
//cotg(theta3) = 0 ==> theta3 = pi/2, cotg(alfa) = 1/tg(alfa)
//tg(theta3) = ( pcm*sin(thetacm) )/( gamma*(pcm*cos(thetacm) + beta*Ecm/c) );
if ( ( pcm*Sin(fThetaCM) )/( gamma*(pcm*Cos(fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[2], 2))) ) > 10000000 ) { //the upper member of the cotg(theta)
fThetaA = TMath::Pi()/2;
}
else {
fThetaA = ATan( (pcm*Sin(fThetaCM))/(gamma*(pcm*Cos(fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[2], 2)))) );
if (fThetaA < 0) {
fThetaA = fThetaA + TMath::Pi();
}
}
if ( ( pcm*Sin(TMath::Pi() - fThetaCM) )/( gamma*(pcm*Cos(TMath::Pi() - fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[3], 2))) ) > 10000000 ) { //the upper member of the cotg(theta)
fThetaB = TMath::Pi()/2;
}
else {
fThetaB = ATan( (pcm*Sin(TMath::Pi() - fThetaCM))/(gamma*(pcm*Cos(TMath::Pi() - fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[3], 2)))) );
if (fThetaB < 0) {
fThetaB = fThetaB + TMath::Pi();
}
}
return 1;
}
void BinaryReaction::SetPhi()
{
//PhiA and PhiB assignment
//fPhiA = ranpokus->Uniform(0., 2.*TMath::Pi());
/*fPhiA = ranPhi.Uniform(0., 2.*TMath::Pi());
if ( fPhiA < TMath::Pi() ) { fPhiB = fPhiA + TMath::Pi(); }
else { fPhiB = fPhiA - TMath::Pi(); }*/
fPhiA = ranPhi.Uniform(0., 2.*TMath::Pi());
fPhiB = fPhiA - TMath::Pi();
}
void BinaryReaction::SetOrientation()
{
if (fThetaA + fThetaIn <= TMath::Pi()) {
fThetaA = fThetaA + fThetaIn;
}
else {
fThetaA = 2*TMath::Pi() - fThetaA + fThetaIn;
if ( fPhiA < TMath::Pi() ) { fPhiA = fPhiA + TMath::Pi(); }
else { fPhiA = fPhiA - TMath::Pi(); }
}
if (fPhiA + fPhiIn > 2*TMath::Pi()) { fPhiA = fPhiA + fPhiIn - 2*TMath::Pi(); }
else { fPhiA = fPhiA + fPhiIn; }
if (fThetaB - fThetaIn >= 0) {
fThetaB = fThetaB - fThetaIn;
}
else {
fThetaB = -1*(fThetaB - fThetaIn);
if ( fPhiB < TMath::Pi() ) { fPhiB = fPhiB + TMath::Pi(); }
else { fPhiB = fPhiB - TMath::Pi(); }
}
if (fPhiB + fPhiIn > 2*TMath::Pi()) { fPhiB = fPhiB + fPhiIn - 2*TMath::Pi(); }
else { fPhiB = fPhiB + fPhiIn; }
}
void BinaryReaction::SetOrientationV()
{
//z' parallel to rdirection; x' is in theta_rdirection plane; y' perpendicular
TVector3 rdirection(0,0,1);
rdirection.SetTheta(fThetaIn);
rdirection.SetPhi(fPhiIn);
//rdirection.SetMag(1.);
TVector3 a(1);
a.SetTheta(fThetaA);
a.SetPhi(fPhiA);
//a.SetMag(1.);
TVector3 b(1);
b.SetTheta(fThetaB);
b.SetPhi(fPhiB);
//b.SetMag(1.);
a.RotateUz(rdirection);
b.RotateUz(rdirection);
fThetaA = a.Theta();
fPhiA = a.Phi();
fThetaB = b.Theta();
fPhiB = b.Phi();
}
void BinaryReaction::Reset()
{
for (Int_t i = 0; i <= 3; i++) {
fM[i] = 0.;
}
fThetaIn = 0.;
fPhiIn = 0.;
fT1 = 0.;
fThetaCM = 0.;
fTa = 0.;
fTb = 0.;
fThetaA = 0.;
fThetaB = 0.;
fPhiA = 0.;
fPhiB = 0.;
return;
}
Double_t BinaryReaction::GetXa(Double_t z)
{
//univerzalni cast
TVector3 v(1., 1., 1.);
v.SetPhi(fPhiA);
v.SetTheta(fThetaA);
// pokud je Theta == pi/2 ==> nejaky jiny vystup
if (v.CosTheta()*z <= 0) { return 0.; }
v.SetMag(TMath::Abs(z/v.CosTheta()));
return v.X();
}
Double_t BinaryReaction::GetYa(Double_t z)
{
//univerzalni cast
TVector3 v(1., 1., 1.);
v.SetPhi(fPhiA);
v.SetTheta(fThetaA);
// pokud je Theta == pi/2 ==> nejaky jiny vystup
if (v.CosTheta()*z <= 0) { return 0.; }
v.SetMag(TMath::Abs(z/v.CosTheta()));
return v.Y();
}
Double_t BinaryReaction::GetXb(Double_t z)
{
//univerzalni cast
TVector3 v(1., 1., 1.);
v.SetPhi(fPhiB);
v.SetTheta(fThetaB);
// pokud je Theta == pi/2 ==> nejaky jiny vystup
if (v.CosTheta()*z <= 0) { return 0.; }
v.SetMag(TMath::Abs(z/v.CosTheta()));
return v.X();
}
Double_t BinaryReaction::GetYb(Double_t z)
{
//univerzalni cast
TVector3 v(1., 1., 1.);
v.SetPhi(fPhiB);
v.SetTheta(fThetaB);
// pokud je Theta == pi/2 ==> nejaky jiny vystup
if (v.CosTheta()*z <= 0) { return 0.; }
v.SetMag(TMath::Abs(z/v.CosTheta()));
return v.Y();
}
/*
* BinaryReaction.h
*
* Created on: 24.5.2010
* Author: Vratislav
*/
//#ifndef BEBINARYREACTION_H_
//#define BEBINARYREACTION_H_
#pragma once
//#include "DllExport.h"
#include <iostream>
#include <TObject.h>
#include <TNamed.h>
#include <TROOT.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TVector3.h>
using std::cout;
using std::endl;
using TMath::Sqrt;
using TMath::Power;
using TMath::Sin;
using TMath::Cos;
using TMath::ATan;
class BinaryReaction : public TNamed {
private:
Double_t fM[4]; //in MeV
//in laboratory system
Double_t fThetaIn; //in rad
Double_t fPhiIn; //in rad
Double_t fT1; //in MeV
//in CM system
Double_t fThetaCM; //in rad
//in laboratory system
Double_t fTa; //in MeV
Double_t fTb; //in MeV
Double_t fThetaA; //in rad
Double_t fThetaB; //in rad
Double_t fPhiA; //in rad
Double_t fPhiB; //in rad
static TRandom3 ranPhi; //!
public:
BinaryReaction();
virtual ~BinaryReaction();
ClassDef(BinaryReaction, 1);
Double_t GetThetaCM() const { return fThetaCM; };
Double_t GetThetaA() const { return fThetaA; };
Double_t GetPhiA() const { return fPhiA; };
Double_t GetTa() const { return fTa; };
Double_t GetMa() const { return fM[2]; };
Double_t GetThetaB() const { return fThetaB; };
Double_t GetPhiB() const { return fPhiB; };
Double_t GetTb() const { return fTb; };
Double_t GetMb() const { return fM[3]; };
Double_t GetXa(Double_t z); //get X projection in focal plane XY + z
Double_t GetXb(Double_t z); //get X projection in focal plane XY + z
Double_t GetYa(Double_t z); //get Y projection in focal plane XY + z
Double_t GetYb(Double_t z); //get Y projection in focal plane XY + z
Int_t FillReaction(Double_t _fThetaIn, Double_t _fT1, Double_t _fM[], Double_t _fThetaCM, Double_t _fPhi);
//osetrit zadavani parametru, ktere maji fyzikalni smysl
//vsude hlidat signum energie
//void FillInput(Double_t _fThetaIn, Double_t _fT1, Double_t _fM[], Double_t _fThetaCM, Double_t fPhi);
Int_t FlatCalculate(); //calculate theta A and B for phi = 0 in reaction frame //predelat pomoci vektoru, pravdepodobne velmi pomala fce
void SetPhi(); //random choice of phi in reaction frame
void SetOrientation(); //does not work properly
void SetOrientationV(); //transform thetas and phis from reaction frame into laboratory frame
//taking into account the projectile flying direction
void SetPhiA(Double_t phia) { fPhiA = phia; }
void SetPhiB(Double_t phib) { fPhiB = phib; }
void SetThetaA(Double_t thetaa) { fThetaA = thetaa; }
void SetThetaB(Double_t thetab) { fThetaB = thetab; }
void SetTa(Double_t ta) { fTa = ta; }
void SetTb(Double_t tb) { fTb = tb; }
void SetM(Int_t i, Double_t m) { fM[i] = m; }
void Reset();
//diagnostics
UInt_t GetPSeed() const { return ranPhi.GetSeed(); };
};
//#endif /* BEBINARYREACTION_H_ */
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class BeEvent;
#pragma link C++ class BeWork;
#pragma link C++ class BinaryReaction;
#pragma link C++ class BeReaction;
#pragma link C++ class BePureEvent;
#endif
//////////////////////////////////////////////////////////
// //
// AnnularDetector //
// //
//
//Some description of this very useful class,
//its properties and a short example how to
//use it. This text may be partly used in
//PhD thesis.
//
//Description of the detector itself.
//
// //
//////////////////////////////////////////////////////////
#include "AnnularDetector.h"
ClassImp(AnnularDetector);
TRandom3 AnnularDetector::ranTheta(1);
TRandom3 AnnularDetector::ranPhi(1);
AnnularDetector::AnnularDetector()
{
}
AnnularDetector::~AnnularDetector()
{
}
Int_t AnnularDetector::GetChMultiplicity(Bool_t sec)
{
//Calculate the number of strips which detected signal in one event
Int_t mult = 0;
if (sec == 1) {
for (Int_t i = 0; i < 32; i++) {
if (fChSec[i] != 0) { mult++; }
}
}
else {
for (Int_t i = 0; i < 32; i++) {
if (fChRing[i] != 0) { mult++; }
}
}
return mult;
}
Int_t AnnularDetector::GetEMultiplicity(Bool_t sec)
{
//Calculate the number of strips which detected signal over particular energy threshold in one event
Int_t mult = 0;
if (sec == 1) {
for (Int_t i = 0; i < 32; i++) {
if (fESec[i] > 0) { mult++; }
}
}
else {
for (Int_t i = 0; i < 32; i++) {
if (fERing[i] > 0) { mult++; }
}
}
return mult;
}
Int_t AnnularDetector::Reset()
{
//Reset variables with energy and angle information
for (Int_t i = 0; i < 32; i++) {
fChRing[i] = 0;
fChSec[i] = 0;
fERing[i] = 0.;
fESec[i] = 0.;
}
fMultChRing = 0;
fMultChSec = 0;
fMultERing = 0;
fMultESec = 0;
return 0;
}
Int_t AnnularDetector::FindSec()
{
//Find first sector where E != 0
Int_t i = 0;
for (i = 0; i < fSecNumber; i++) {
if (fESec[i] > 0.) { break; }
}
if (i == fSecNumber) {
return (-1);
}
return i;
}
Int_t AnnularDetector::FindRing()
{
//Find first ring where E != 0
if (fDSD == kFALSE) { Warning("FindRing", "May not be used for this detector (SSD)"); }
Int_t i = 0;
for (i = 0; i < fRingsNumber; i++) {
if (fERing[i] > 0.) { break; }
}
if (i == fRingsNumber) {
return -1;
}
return i;
}
Double_t AnnularDetector::GetTheta(Int_t ring)
{
//Get the angle theta in radians which was measured in the ring.
if (ring >= fRingsNumber) {
Warning("GetTheta", "Number of rings in this detector is equal to %d.", fRingsNumber);
return 0.; //zamyslet se
}
if (ring < 0) {
Warning("GetTheta", "You cannot input negative ring number.");
return 0.; //zamyslet se
}
Double_t theta;
theta = TMath::ATan(( fInnerR + ( ring+ranTheta.Uniform(1) )*(fOuterR - fInnerR)/fRingsNumber )/fZ);
return theta;
}
Double_t AnnularDetector::GetPhi(Int_t sec)
{
//Get the angle phi in radians which was measured in the sector sec.
if (sec >= fSecNumber) {
Warning("GetPhi", "Number of sectors in this detector is equal to %d.", fRingsNumber);
return -1.; //zamyslet se
}
if (sec < 0) {
Warning("GetPhi", "You cannot input negative sector number.");
return -1.; //zamyslet se
}
// Double_t theta;
// theta = TMath::ATan(( fInnerR + ( ring+ranTheta.Uniform(1) )*(fOuterR - fInnerR)/fRingsNumber )/fZ);
return (2*TMath::Pi()/fSecNumber)*( sec + ranTheta.Uniform(1) );
}
Double_t AnnularDetector::GetX(Int_t ring, Int_t sec) {
TVector2 flat; //vector in detector plane
Double_t r = fInnerR + ( ring+ranTheta.Uniform(1) )*(fOuterR - fInnerR)/fRingsNumber;
// Double_t phi = 2 * TMath::Pi()
// - (2*TMath::Pi()/fSecNumber)*( sec + ranTheta.Uniform(1) )
// - TMath::Pi() / 2.;
Double_t phi = //2 * TMath::Pi()
- (2 * TMath::Pi() / fSecNumber) * (sec + ranTheta.Uniform(1))
+ TMath::Pi() / 2.;
// printf("sec: %d\tphi: %3.2f\n", sec, phi);
flat.SetMagPhi(r, phi);
// printf("sec: %d\tphi: %3.2f\tphiV: %3.2f\t\t", sec, phi*TMath::RadToDeg(), flat.Phi()*TMath::RadToDeg());
return flat.X();
}
Double_t AnnularDetector::GetY(Int_t ring, Int_t sec) {
TVector2 flat; //vector in detector plane
Double_t r = fInnerR + ( ring+ranTheta.Uniform(1) )*(fOuterR - fInnerR)/fRingsNumber;
Double_t phi = //2 * TMath::Pi()
- (2*TMath::Pi()/fSecNumber)*( sec + ranTheta.Uniform(1) )
+ TMath::Pi() / 2.;
flat.SetMagPhi(r, phi);
return flat.Y();
}
void AnnularDetector::InitDetector(Bool_t dsd, Int_t nosecs, Int_t norings,
Double_t innerr, Double_t outerr, Double_t thick,
Double_t deadthick, Double_t z)
{
//dsd = 1 for DSSD; dsd = 0 for SSSD
//nosecs - number of sectors
//norings - number of rings
//innerr - inner radius in mm
//outerr - outer radius in mm
//thick - thickness of the detector in microns
//deadthick - thickness of the dead layer in microns
//z - distance of the detector from the target
SetDSD(dsd);
SetSecNumber(nosecs);
if (fDSD == 1) { SetRingsNumber(norings); }
SetInnerR(innerr);
SetOuterR(outerr);
SetThickness(thick);
SetDeadLayerThickness(deadthick);
SetZ(z);
return;
}
/*
* AnnularDetector.h
*
* Created on: 8.12.2009
* Author: Vratislav
*/
#pragma once
#include <TObject.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TVector2.h>
//docasne
#include<iostream>
using std::cout;
using std::endl;
class AnnularDetector /*: public TObject*/ {
private:
/////vsechno se bude zadavat ze souboru
Bool_t fDSD; //!
Int_t fRingsNumber; //!
Int_t fSecNumber; //!
Double_t fInnerR; //! //mm, sensitive region
Double_t fOuterR; //! //mm, sensitive region
Double_t fThickness; //! //mcm
Double_t fDeadLayerThickness; //! //mcm,
Double_t fZ; //! //mm
// Double_t fPhi[32]; //
static TRandom3 ranTheta;
static TRandom3 ranPhi;
public: //docasne verejne
Int_t fChSec[32];
Int_t fChRing[32];
Double_t fESec[32];
Double_t fERing[32];
public:
Int_t fMultChSec;
Int_t fMultChRing;
Int_t fMultESec;
Int_t fMultERing;
public:
AnnularDetector();
virtual ~AnnularDetector();
void SetDSD(Bool_t dsd) { fDSD = dsd; } //prakticky vsechny settery
void SetRingsNumber(Int_t norings) { fRingsNumber = norings; } //jsou v rootu virtual
void SetSecNumber(Int_t nosec) { fSecNumber = nosec; }
void SetInnerR(Double_t r1) { fInnerR = r1; }
void SetOuterR(Double_t r2) { fOuterR = r2; }
void SetThickness(Double_t x) { fThickness = x; }
void SetDeadLayerThickness(Double_t dx) { fDeadLayerThickness = dx; }
void SetZ(Double_t z) { fZ = z; }
void InitDetector(Bool_t dsd, Int_t nosecs, Int_t norings,
Double_t innerr, Double_t outerr, Double_t thick,
Double_t deadthick, Double_t z); //dsd, nor, nos, ir(mm), or(mm), x(mcm), xdead(mcm), z(mm)
Int_t GetSecNumber() { return fSecNumber; }
Int_t GetRingsNumber() { return fRingsNumber; }
Double_t GetThickness() { return fThickness; }
Double_t GetDeadLayerThickness() { return fDeadLayerThickness; }
Double_t GetZPosition() { return fZ; }
Int_t GetChMultiplicity(Bool_t sec = 1);
Int_t GetEMultiplicity(Bool_t sec = 1);
Int_t Reset(); //reset information relating to particular event
Int_t FindSec(); //returns first sector where E != 0, returns -1 if E == 0 for all sectors
Int_t FindRing(); //returns first ring where E != 0, returns -1 if E == 0 for all rings
Double_t GetTheta(Int_t ring);
Double_t GetPhi(Int_t sec);
Double_t GetX(Int_t ring, Int_t sec);
Double_t GetY(Int_t ring, Int_t sec);
ClassDef(AnnularDetector, 1);
};
################################################################################
# Detectors input with some variables
################################################################################
DETECTORSLIBS := -lCint -lCore -lMathCore
# Add inputs and outputs from these tool invocations to the build variables
DETECTORS_HEADERS += \
$(DETECTORS)/AnnularDetector.h \
$(DETECTORS)/linkdef.h
DETECTORSCPP_SRCS += \
$(DETECTORS)/AnnularDetector.cpp \
$(DETECTORS)/DetectorsCint.cpp
DETECTORSOBJS += \
$(DETECTORS)/AnnularDetector.o \
$(DETECTORS)/DetectorsCint.o
DETECTORSCPP_DEPS += \
$(DETECTORS)/AnnularDetector.d \
$(DETECTORS)/DetectorsCint.d
################################################################################
# Detectors input with some variables
################################################################################
DETECTORSLIBS := -lCint -lCore -lMathCore
# Add inputs and outputs from these tool invocations to the build variables
DETECTORSCPP_SRCS += \
./AnnularDetector.cpp \
./DetectorsCint.cpp
DETECTORSOBJS += \
./AnnularDetector.o \
./DetectorsCint.o
DETECTORSCPP_DEPS += \
./AnnularDetector.d \
./DetectorsCint.d
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class AnnularDetector;
#endif
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class AnnularDetector;
//#pragma link C++ class AculConvert;
//#pragma link C++ class AculCalibration;
//#pragma link C++ class AculCalibratedEvent;
#endif
!======================================================================
!== Energy loss in material
!==
!== File: ELOSS.f90
!== Date: 08/12/2000
!======================================================================
! ELOSS.f90
!
! FUNCTIONS/SUBROUTINES exported from ELOSS.dll:
! ELOSS - subroutine
!
subroutine ELOSS(NEL,ZEL,AEL,WEL,DEN,ZP,AP,NE,ETAB,RE,ZW,AW)
! Expose subroutine ELOSS to users of this DLL
!
!DEC$ ATTRIBUTES DLLEXPORT::ELOSS
IMPLICIT REAL*8(A-H,O-Z)
INTERFACE
SUBROUTINE rde(E, R, R1, R2, IRE)
REAL*8 E, R, R1, R2
INTEGER IRE
END SUBROUTINE rde
END INTERFACE
! Variables
! Input:
INTEGER NEL, NE
REAL*8 ZEL(NEL),AEL(NEL),WEL(NEL), DEN,ZP,AP, ETAB(NE)
! Output:
REAL*8 RE(NE),ZW,AW
! Local
COMMON adj(5),as(5),z(5),f(5),um(5),oz,ro,azm,aj,zion,aion,&
zmed,amed,jc
! Body of ELOSS
zion = ZP
aion = AP
jc = NEL
ro = DEN
!== two parameters with question
ire = 1
oz = 1.
DO i = 1, NEL
as(i) = AEL(i)
z(i) = ZEL(i)
um(i) = WEL(i)
IF(z(i) .LE. 13.) THEN
adj(i) = 12.*z(i) + 7.
ELSE
adj(i) = 9.76*z(i) + 58.8 / z(i)**0.19
END IF
END DO
uma=0.
zam=0.
umed=0.
zme=0.
aj=0.
DO i = 1, jc
uma=uma+um(i)*as(i)
umed=umed+um(i)
zme=zme+um(i)*z(i)
END DO
amed=uma/umed
zmed=zme/umed
DO i = 1, jc
f(i)=um(i)*as(i)/uma
END DO
DO i = 1, jc
zam=zam+f(i)*z(i)/as(i)
END DO
azm=1./zam
aj=0.
DO i = 1, jc
aj=aj+f(i)*z(i)*log(adj(i))/as(i)
END DO
aj=aj*azm
aj=exp(aj)
ZW = zmed
AW = amed
DO j = 1, NE
!== result [mg/cm^2]
CALL rde(ETAB(j), RE(j), rel1, rel, ire)
END DO
end subroutine ELOSS
SUBROUTINE rde(e,range,rel1,rel,ix)
IMPLICIT REAL*8(A-H,O-Z)
INTERFACE
REAL*8 FUNCTION c(x)
REAL*8 x
END FUNCTION c
REAL*8 FUNCTION c1(x)
REAL*8 x
END FUNCTION c1
END INTERFACE
! calculates range and de/dx for compounds
dimension ala(3),ala1(3)
dimension a(3,3),a1(4,4),a2(4,4),b(2,3),cc(5)
common adj(5),as(5),z(5),f(5),um(5),oz,ro,azm,aj,zion,aion,&
zmed,amed,jc
dimension alaj1(4),altau1(4)
data a/ -.75265, .073736, .040556, 2.5398,-.312,.018664, &
-.24598, .11548, -.0099661/
data a1/-8.0155, 0.36916, -1.4307e-02, 3.4718e-03, &
1.8371, -1.452e-02, -3.0142e-02, 2.3603e-03, &
0.045233, -9.5873e-04, 7.1303e-03, -6.8538e-04, &
-5.9898e-03,-5.2315e-04, -3.3802e-04, 3.9405e-05/
data a2/-8.725, 0.8309, -0.13396, 0.012625, &
1.8797, 0.11139, -0.064808, 0.0054043, &
0.74192, -0.528805, 0.1264232, -0.00934142, &
0.752005, -0.5558937, 0.12843119, -0.009306336/
data b/ 0.422377e-06, 3.858019e-09, 0.0304043e-06, -0.1667989e-09, &
-0.00038106e-06, 0.00157955e-09/
if(e.le.0.002)then
range=0.
rel1=0.
rel=0.
return
endif
! this a flag
ibis=-1
en = e/aion
tau = en*1.008
altau = log(tau)
alaj = log(aj)
alaj1(1)= 1.
altau1(1)= 1.
DO kk = 2,4
alaj1(kk)=alaj**(kk-1)
altau1(kk)=altau**(kk-1)
END DO
t = tau/938.256
tt = 2.*t + t*t
beta= sqrt(tt)/(1.+t)
s1=0.
DO i=1,4
DO j=1,4
s1 = s1 + a2(j,i)*alaj1(j)*altau1(i)
END DO
END DO
ala(1) = azm*exp(s1)
s2=0.
DO i=2,4
DO j=1,4
s2 = s2 + (i-1)*a2(j,i)*alaj1(j)*altau1(i-1)
END DO
END DO
ala1(1)=ala(1)*s2/tau
s1=0.
DO i=1,4
DO j=1,4
s1 = s1 + a1(j,i)*alaj1(j)*altau1(i)
END DO
END DO
ala(3) = azm*exp(s1)
s2=0.
DO i=2,4
DO j=1,4
s2 = s2 + (i-1)*a1(j,i)*alaj1(j)*altau1(i-1)
END DO
END DO
ala1(3)=ala(3)*s2/tau
s1=0.
DO i=1,3
DO j=1,3
s1 = s1 + a(j,i)*alaj1(j)*altau1(i)
END DO
END DO
ala(2)=azm*exp(s1)/1000.
s2=0.
DO i=2,3
DO j=1,3
s2 = s2 + (i-1)*a(j,i)*alaj1(j)*altau1(i-1)
END DO
END DO
ala1(2)=ala(2)*s2/tau
25 continue
coefa=3.
coefb=1.
alaa=(ala(1)*(tanh(coefa*(.98 - en))+1.)/2. &
+ ala(2)*(tanh(coefa*(en - .98))+1.)/2.) &
* (tanh(coefb*(8.0 - en))+1.)/2. &
+ ala(3)*(tanh(coefb*(en - 8.))+1.)/2.
alim1=0.
alim2=0.
if(coefa*(en-.98).lt.85)then
alim1=1.008/cosh(coefa*(.98-en))/cosh(coefa*(.98-en))
endif
if(coefb*(en-8.).lt.85)then
alim2=1.008/cosh(coefb*(8.-en))/cosh(coefb*(8.-en))
endif
alaa1=(ala1(1)*(tanh(coefa*(.98-en))+1.)/2.+ &
ala1(2)*(tanh(coefa*(en-.98))+1.)/2.)* &
(tanh(coefb*(8.-en))+1.)/2.+ &
ala1(3)*(tanh(coefb*(en-8.))+1.)/2.+ &
coefa/2.*(tanh(coefb*(8.-en))+1.)/2.* &
(ala(2)*alim1-ala(1)*alim1)+ &
coefb/2.*(ala(3)*alim2- &
(ala(1)*(tanh(coefa*(.98-en))+1.)/2.+ &
ala(2)*(tanh(coefa*(en-.98))+1.)/2.)*alim2)
hi=137.*beta/zion
bz=(31.8+3.86*(aj**.625))*azm*.000001
bz=bz*(zion**2.666667)*c(hi)
bz1=(4.357+.5288*(aj**.625))*azm*.001
bz1=bz1*(zion**1.666667)*c1(hi)
bep=beta*beta
rel1=zion*zion/(alaa1+bz1*((1.-bep)**1.5)/931.141/beta)
rel1=rel1/1000.
range=(alaa+bz)*aion/1.008/zion/zion
range=range*1000.
! Atention!! this version do not work correctly for ix.ne.1
if(ix.ne.1)return
ank=.153574*ro/azm
z23=zion**.666667
abet=beta*125./z23
zef=zion*(1.-exp(-abet))
zef2=zef*zef
om=1022000.*bep/(1.-bep)
cbet=0.
DO k=1,jc
cc(k)=0.
DO i=1,2
DO j=1,3
cc(k)=cc(k)+b(i,j)*((1./bep-1.)**j)*((adj(k)**(i+1)))
END DO
END DO
cbet=cbet+f(k)*cc(k)/as(k)
END DO
cbet=cbet*azm
del1=ank*zef2*(log(om/oz)-bep)/bep
del1=del1/ro
del2=2.*ank*zef2*(log(oz/aj)-cbet)/bep
del2=del2/ro
rel2=rel1-del1
rel3=rel1-rel2+del2
if(del1)8,8,9
8 rel=rel1
goto 12
9 if(del1+del2-rel2)10,10,11
10 rel=rel3
goto 12
11 if(del1.lt.rel1)then
rel=rel2
else
rel=rel1
endif
12 return
END
REAL*8 FUNCTION c(x)
REAL*8 x
IF(x .LE. 0.2) THEN
c = -0.00006 + (0.05252 + 0.12857*x)*x
ELSE IF(x .LE. 2.) THEN
c = -0.00185 + (0.07355 + (0.07172 - 0.02723*x)*x)*x
ELSE IF(x .LE. 3.) THEN
c = -0.0793 + (0.3323 - (0.1234 - 0.0153*x)*x)*x
ELSE
c = 0.22
END IF
END
REAL*8 FUNCTION c1(x)
REAL*8 x
IF(x .LE. 0.2) THEN
c1 = 0.05252+.25714*x
ELSE IF(x .LE. 2.0) THEN
c1 = 0.07355 + (0.14344 - 0.08169*x)*x
ELSE IF(x .LE. 3.0) THEN
c1 = 0.3323 - (0.2468 - 0.0459*x)*x
ELSE
c1 = 0.
END IF
END
#include "TELoss.h"
ClassImp(TELoss)
extern "C" void eloss_(int *nel, double *zel, double *ael, double *wel, double *den
,double *zp, double *ap, int *ne, double *etab, double *rtab
,double *zw, double *aw);
//extern "C" __declspec( dllimport ) void _stdcall
// ELOSS(int *nel, double *zel, double *ael, double *wel, double *den
// ,double *zp, double *ap, int *ne, double *etab, double *rtab
// ,double *zw, double *aw);
TELoss::TELoss(void)
{
ne = 0;
mel = 0;
} //-----------------------------------------------------------------
TELoss::TELoss(int _mel, double _den)
{
ne = 0;
mel = 0;
SetEL(_mel,_den);
} //-----------------------------------------------------------------
TELoss::TELoss(TELoss &teloss)
{
mel = teloss.mel;
nel = teloss.nel;
// zel = teloss.zel;
// ael = teloss.ael;
// wel = teloss.wel;
zel = TArrayD(teloss.zel);
ael = TArrayD(teloss.ael);
wel = TArrayD(teloss.wel);
den = teloss.den;
zp = teloss.zp;
ap = teloss.ap;
ne = teloss.ne;
// cout << etab.GetSize() << endl;
// etab = teloss.etab;
// rtab = teloss.rtab;
// detab = teloss.detab;
etab = TArrayD(teloss.etab);
rtab = TArrayD(teloss.rtab);
detab = TArrayD(teloss.detab);
zw = teloss.zw;
aw = teloss.aw;
a = teloss.a;
b = teloss.b;
} //-----------------------------------------------------------------
TELoss::~TELoss()
{
} //-----------------------------------------------------------------
void TELoss::SetEL(int _mel, double _den)
{
den = _den;
if(_mel != mel)
{
mel = _mel;
zel.Set(mel);
ael.Set(mel);
wel.Set(mel);
}
nel = 0;
} //-----------------------------------------------------------------
int TELoss::AddEL(double _zel, double _ael, double _wel)
{
//===============================================================
//== Add one element
//===============================================================
if(nel >= mel) return -1;
zel[nel] = _zel;
ael[nel] = _ael;
wel[nel] = _wel;
return ++nel;
} //-----------------------------------------------------------------
void TELoss::SetZP(double _zp, double _ap)
{
zp = _zp;
ap = _ap;
} //-----------------------------------------------------------------
void TELoss::SetEtab(int _ne, double e2)
{
//===============================================================
//== Set list of energies and calculate R(E)
//===============================================================
// cout << "sem" << endl;
int i;
double e1 = 0.;
double es;
double coef = 10. / den; //== convertion mg/cm^2 -> micron
// cout << "sem" << endl;
// cout << _ne << endl;
// cout << ne << endl;
if(_ne != ne)
{
ne = _ne; //asi pocet elementu
// cout << "sem" << endl;
// cout << ne << endl;
etab.Set(ne);
// cout << "sem" << endl;
rtab.Set(ne);
}
// cout << "sem" << endl;
if ( ((e2-e1)/(static_cast<double>(ne))) > 0.01) { Warning("SetEtab", "Table resolution can be too low for some cases."); }
es = (e2 - e1) / (double)(ne-1); //energeticky skok v tabulce
for(i = 0; i < ne; ++i)
{
etab[i] = e1;
e1 += es;
}
eloss_(&nel, zel.GetArray(), ael.GetArray(), wel.GetArray(), &den
,&zp, &ap
,&ne, etab.GetArray(), rtab.GetArray(), &zw, &aw);
for(i = 0; i < ne; ++i) rtab[i] = coef * rtab[i];
//if(RE) delete RE;
//if(ER) delete ER;
//RE = new TSpline3("RE",etab.GetArray(),rtab.GetArray(),ne);
//ER = new TSpline3("ER",rtab.GetArray(),etab.GetArray(),ne);
} //-----------------------------------------------------------------
void TELoss::SetDeltaEtab(double r)
{
//===============================================================
//== Calculate and set list of dE(E)
//===============================================================
int i;
// double es;
//double coef = 10. / den; //== convertion mg/cm^2 -> micron
if(ne == 0) {
Warning("SetDeltaEtab", "Etab was not initialized yet");
return;
}
detab.Set(ne);
// if ( ((e2-e1)/(static_cast<double>(ne))) > 0.01) { Warning("SetEtab", "Table resolution can be too low for some cases."); }
// es = (e2 - e1) / (double)(ne-1); //energeticky skok v tabulce
for(i = 0; i < ne; ++i)
{
if ( GetE(etab[i], r) != 0. ) {
detab[i] = etab[i] - GetE(etab[i], r);
// e1 += es;
}
else { detab[i] = 0.; }
}
// eloss_(&nel, zel.GetArray(), ael.GetArray(), wel.GetArray(), &den
// ,&zp, &ap
// ,&ne, etab.GetArray(), rtab.GetArray(), &zw, &aw);
//
// for(i = 0; i < ne; ++i) rtab[i] = coef * rtab[i];
return;
}
double TELoss::GetR(double e0,double e)
{
//==================================================================
//== Calculate new energy using linear interpolation
//==================================================================
//double r0 = RE->Eval(e0);
//return ER->Eval(r0 - r);
double r0 = linear(ne, etab.GetArray(), rtab.GetArray(), e0);
double r = linear(ne, etab.GetArray(), rtab.GetArray(), e);
// return aitken3(ne, rtab.GetArray(), etab.GetArray(), r0-r);
return r0 - r;
} //--------------------------------------------------------------------
double TELoss::GetRold(double e0,double e)
{
//==================================================================
//== Calculate new energy
//==================================================================
//double r0 = RE->Eval(e0);
//return ER->Eval(r0 - r);
double r0 = aitken3(ne, etab.GetArray(), rtab.GetArray(), e0);
double r = aitken3(ne, etab.GetArray(), rtab.GetArray(), e);
// return aitken3(ne, rtab.GetArray(), etab.GetArray(), r0-r);
return r0 - r;
} //--------------------------------------------------------------------
double TELoss::GetE(double e0,double r)
{
//==================================================================
//== Calculate new energy using linear interpolation
//==================================================================
// Double_t r0 = ylinear(ne, e0);
Double_t r0 = linear(ne, etab.GetArray(), rtab.GetArray(), e0);
if (r0 < r) { return 0; }
return linear(ne, rtab.GetArray(), etab.GetArray(), r0-r);
} //--------------------------------------------------------------------
double TELoss::GetE0dE(double de, double r)
{
//==================================================================
//== Calculate new energy from deltaE using linear interpolation
//==================================================================
SetDeltaEtab(r); //problem s interpolaci
if ( de > TMath::MaxElement(ne, detab.GetArray()) ) { return -1; }
int i0 = 0;
while ( detab[i0] == 0 ) {
i0++;
}
//vypis
// cout << "i0 = \t" << i0 << "\tdetab[i0] = \t" << detab[i0] << "\tetab[i0] = \t" << etab[i0] << endl;
int i1 = ne-1;
int ix;
while(i0 < i1 - 1)
{
ix = (i0 + i1) / 2;
if ( de > detab[ix] ) i1 = ix;
else i0 = ix;
}
// if (r0 < r) { return 0; }
// cout << "r0 =\t\t" << r0 << endl;
// cout << "r0-r =\t\t" << r0-r << endl;
// cout << "E(r0-r) =\t\t" << xlinear(ne, r0-r) << endl;
// cout << "E(r0-r) =\t\t" << linear(ne, rtab.GetArray(), etab.GetArray(), r0-r) << endl;
return linear(ne, detab.GetArray(), etab.GetArray(), i0, de);
} //--------------------------------------------------------------------
double TELoss::GetE0dE(double de)
{
if ( de > TMath::MaxElement(ne, detab.GetArray()) ) { return -1; }
int i0 = 0;
while ( detab[i0] == 0 ) {
i0++;
}
// cout << "detab[i0] je " << detab[i0] << endl;
int i1 = ne-1;
int ix;
while(i0 < i1 - 1)
{
ix = (i0 + i1) / 2;
if ( de > detab[ix] ) i1 = ix;
else i0 = ix;
}
// cout << "detab[i0] je " << detab[i0] << endl
// << "detab[i1] je " << detab[i1] << endl;
// if (r0 < r) { return 0; }
// cout << "r0 =\t\t" << r0 << endl;
// cout << "r0-r =\t\t" << r0-r << endl;
// cout << "E(r0-r) =\t\t" << xlinear(ne, r0-r) << endl;
// cout << "E(r0-r) =\t\t" << linear(ne, rtab.GetArray(), etab.GetArray(), r0-r) << endl;
return linear(ne, detab.GetArray(), etab.GetArray(), i0, de);
//return linear(ne, detab.GetArray(), etab.GetArray(), de);
}
double TELoss::GetEold(double e0,double r)
{
//==================================================================
//== Calculate new energy using aitken interpolation
//==================================================================
//double r0 = RE->Eval(e0);
//return ER->Eval(r0 - r);
double r0 = aitken3(ne, etab.GetArray(), rtab.GetArray(), e0);
return aitken3(ne, rtab.GetArray(), etab.GetArray(), r0-r);
} //--------------------------------------------------------------------
double TELoss::GetE0(double e0,double r)
{
//==================================================================
//== Calculate new energy using linear interpolation
//==================================================================
Double_t r0 = linear(ne, etab.GetArray(), rtab.GetArray(), e0);
return linear(ne, rtab.GetArray(), etab.GetArray(), r0+r);
} //--------------------------------------------------------------------
double TELoss::GetE0old(double e,double r)
{
//==================================================================
//== Calculate new energy
//==================================================================
//double r0 = RE->Eval(e);
//return ER->Eval(r0 + r);
double r0 = aitken3(ne, etab.GetArray(), rtab.GetArray(), e);
return aitken3(ne,rtab.GetArray(),etab.GetArray(),r0+r);
} //-----------------------------------------------------------------
void TELoss::PrintRE()
{
int i;
for(i = 0; i < ne; ++i)
{
printf("%10.2lf %10.3lf\n",etab[i],rtab[i]);
}
} //-----------------------------------------------------------------
void TELoss::PrintdEE()
{
int i;
for(i = 0; i < ne; ++i)
{
printf("%10.3lf %10.3lf\n",etab[i],detab[i]);
}
} //-----------------------------------------------------------------
void TELoss::PrintREtoFile()
{
FILE *fw;
fw = fopen("outputRE.txt", "w");
int i;
for(i = 0; i < ne; ++i)
{
fprintf(fw, "%10.2lf \t%10.3lf\n",etab[i],rtab[i]);
}
fclose(fw);
}
void TELoss::PrintdEEtoFile(const char* outfile)
{
FILE *fw;
fw = fopen(outfile, "w");
int i;
for(i = 0; i < ne; ++i)
{
fprintf(fw, "%10.3lf \t%10.3lf\n",etab[i],detab[i]);
}
fclose(fw);
}
double TELoss::aitken3(int ntab, double *xtab, double *ytab, double x)
{
//====================================================================
//== Interpolation by 4 points
//====================================================================
// int ndeg = 4;
double y[4];
int i0,i;
i0 = bsearch(ntab,xtab,x) - 1;
if(i0 < 0) i0 = 0;
else if(i0 + 3 >= ntab) i0 = ntab - 4;
for(i = 0; i < 4; ++i) y[i] = ytab[i0 +i];
return aitken(3,xtab+i0,y,x);
} //-----------------------------------------------------------------
int TELoss::bsearch(int ntab, double *xtab, double x)
{
int i0 = 0;
int i1 = ntab;
int ix;
while(i0 < i1 - 1)
{
ix = (i0 + i1) / 2;
if(x < xtab[ix]) i1 = ix;
else i0 = ix;
}
return i0;
} //-----------------------------------------------------------------
double TELoss::aitken(int ntab, double *xtab, double *ytab, double x)
{
int i,j;
for(j = 0; j < ntab; ++j)
{
for(i = j+1; i <= ntab; ++i)
{
ytab[i] = ((x - xtab[j])*ytab[i] - (x - xtab[i])*ytab[j])
/(xtab[i] - xtab[j]);
}
}
return ytab[ntab];
} //-----------------------------------------------------------------
Double_t TELoss::linear(int ntab, Double_t *xtab, Double_t *ytab, double x)
{
Double_t y;
Double_t a = 0;
Double_t b = 0;
int i0 = 0;
i0 = bsearch(ntab, xtab, x) - 1;
if (i0 < 0) { i0 = 0; }
else {
if ( i0+2 >= ntab ) { i0 = i0 - 1; }
}
// cout << "xtab[i0] = \t" << xtab[i0] << "\txtab[i0+2] = \t" << xtab[i0+2] << endl;
// cout << "ytab[i0] = \t" << ytab[i0] << "\tytab[i0+2] = \t" << ytab[i0+2] << endl;
a = ( ytab[i0] - ytab[i0+2] )/( xtab[i0] - xtab[i0+2] );
b = ( ytab[i0+2]*xtab[i0] - ytab[i0]*xtab[i0+2] )/( xtab[i0] - xtab[i0+2] );
Double_t lerror = TMath::Abs(( (1/a)*ytab[i0+1] - (b/a) ) - xtab[i0+1]);
if (lerror > 0.01) {
Warning("linear", "Linear interpolation does not give sufficient precision: %f", lerror);
cout << "xtab[i0+1] \t" << xtab[i0+1] << "\tytab[i0+1] \t" << ytab[i0+1] << endl;
}
y = a*x + b;
return y;
}
Double_t TELoss::linear(int ntab, Double_t *xtab, Double_t *ytab, Int_t i0, double x)
{
Double_t y;
Double_t a = 0;
Double_t b = 0;
//int i0 = 0;
//i0 = bsearch(ntab, xtab, x) - 1;
if (i0 < 0) { i0 = 0; }
else {
if ( i0+2 >= ntab ) { i0 = i0 - 1; }
}
// cout << "xtab[i0] = \t" << xtab[i0] << "\txtab[i0+2] = \t" << xtab[i0+2] << endl;
// cout << "ytab[i0] = \t" << ytab[i0] << "\tytab[i0+2] = \t" << ytab[i0+2] << endl;
a = ( ytab[i0] - ytab[i0+2] )/( xtab[i0] - xtab[i0+2] );
b = ( ytab[i0+2]*xtab[i0] - ytab[i0]*xtab[i0+2] )/( xtab[i0] - xtab[i0+2] );
Double_t lerror = TMath::Abs(( (1/a)*ytab[i0+1] - (b/a) ) - xtab[i0+1]);
if (lerror > 0.001) { Warning("linear_ext", "Linear interpolation does not give sufficient precision: %f", lerror); }
// cout << "a je " << a << "\tb je " << b << endl;
y = a*x + b;
return y;
}
#pragma once
#include <TROOT.h>
#include <TObject.h>
#include <TSpline.h>
#include <TArrayD.h>
#include <TMath.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using std::cout;
using std::endl;
class TELoss : public TObject
{
private:
int mel;
int nel;
TArrayD zel, ael, wel;
double den;
double zp,ap;
int ne;
TArrayD etab, rtab, detab; //tables: energy, range, deltaE
double zw,aw;
double a, b; //linear interpolation coeficients, y = a*x + b
int bsearch(int ntab, double *xtab, double x);
double aitken(int ntab, double *xtab, double *ytab, double x);
double aitken3(int ntab, double *xtab, double *ytab, double x);
Double_t linear(int ntab, Double_t *xtab, Double_t *ytab, double x);
Double_t linear(int ntab, Double_t *xtab, Double_t *ytab, Int_t i0, double x);
public:
TELoss(void);
TELoss(int _mel, double _den);
TELoss(TELoss &);
virtual ~TELoss(void);
void SetEL(int _mel, double _den);
void SetDen(double _den) {den = _den;}
int AddEL(double _zel, double _ael, double _wel=1.);
void SetZP(double _zp, double _ap);
void SetEtab(int _ne, double e2); //from 0 MeV to e2 MeV, _ne elements
void SetDeltaEtab(double r); //pro kazdou tlostku bude zvlastni tabulka, nebo taky ne
double GetZ() { return zp; };
double GetA() { return ap; };
double GetE(double e0, double r);
double GetE0dE(double de, double r); //de in MeV, r in microns
double GetE0dE(double de); //de in MeV
double GetE0(double e, double r);
double GetEold(double e0, double r);
double GetE0old(double e, double r);
double GetR(double e0, double e);
double GetRold(double e0, double e);
void PrintRE();
void PrintdEE();
void PrintREtoFile();
void PrintdEEtoFile(const char* outfile = "outputdEE.txt");
ClassDef(TELoss, 1)
};
################################################################################
# TELoss input with some variables
################################################################################
TELOSSLIBS := -lCore -lCint -lMathCore -lMatrix -lHist -lgfortran
# Add inputs and outputs from these tool invocations to the build variables
TELOSSCPP_SRCS += \
$(TELOSS)/TELoss.cpp \
$(TELOSS)/TELossCint.cpp
TELOSSOBJS += \
$(TELOSS)/ELOSS.o \
$(TELOSS)/TELoss.o \
$(TELOSS)/TELossCint.o
TELOSSCPP_DEPS += \
$(TELOSS)/TELoss.d \
$(TELOSS)/TELossCint.d
{
gSystem->Load("/media/commondata/work/Eclipse/TELoss/Debug/libTELoss.so");
// TELoss target;
TELoss mSi;
TELoss mSi3;
// target.SetEL(1, 0.0014625);
// target.AddEL(1., 3., 1.);
// target.SetZP(1., 1.);
// target.SetEtab(1000, 0., 50.);
mSi.SetEL(1, 2.321);
mSi.AddEL(14., 28.086, 1);
// mSi.SetZP(1., 1.); //protony
mSi.SetZP(2., 4.); //alfy
mSi.SetEtab(100000, 200.); //deltaE
mSi.SetDeltaEtab(300);
mSi3.SetEL(1, 2.321);
mSi3.AddEL(14., 28.086, 1);
// mSi3.SetZP(1., 1.); //protony
mSi3.SetZP(2., 3.); //3He
mSi3.SetEtab(100000, 200.); //deltaE
mSi3.SetDeltaEtab(300);
// TCanvas *c1 = new TCanvas("jmeno", "titulek", 617, 0, 1058, 972);
// mSi.PrintdEEtoFile("out.txt");
// TGraph *gr1 = new TGraph("out.txt");
// gr1->Draw("A*");
//
// mSi3.PrintdEEtoFile("out3.txt");
// TGraph *gr2 = new TGraph("out3.txt");
// gr2->SetMarkerColor(2);
// gr2->Draw("*");
// mSi.PrintREtoFile();
// TGraph *gr2 = new TGraph("outputRE.txt");
// gr2->Draw("A*");
// char out[20], gr[10];
// TObjArray glist(0);
// TMultiGraph *mg = new TMultiGraph();
// TGraph *g;
// for (Int_t i = 0; i < 30; i++) {
// sprintf(out, "output%d.txt", i);
// cout << out << endl;
// mSi.SetDeltaEtab(300+i);
// mSi.PrintdEEtoFile(out);
// g = new TGraph(out);
// glist.Add(g);
// cout << out << endl;
// }
// cout << mSi.GetE(7.686, 10) << endl;
// cout << mSi.GetE(7.686, 1) << endl;
// Double_t finalEnergy = -1.;
// Double_t Energy = steelwin.GetE0(finalEnergy, 30.);
// Energy = target.GetE0(Energy, 2000.);
// printf("%f\n", Energy);
}
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class TELoss;
#endif
\ No newline at end of file
{
// gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo","*", "TStreamerInfo", "RIO", "TStreamerInfo()");
gSystem->Load("~/work/makefilesBe/libAculData.so");
gSystem->Load("~/work/makefilesBe/libTELoss.so");
gSystem->Load("~/work/makefilesBe/libDetectors.so");
gSystem->Load("~/work/makefilesBe/libDetectors.so");
gSystem->Load("~/work/makefilesBe/libBe.so");
THtml h;
h.SetInputDir("~/work/makefilesBe");
// h.MakeIndex("");
// h.MakeClass("TELoss");
// h.MakeClass("TELoss");
// h.MakeClass("TELoss");
h.MakeAll();
/*h.SetInputDir("/media/commondata/work/Eclipse/AculData/Debug:/media/commondata/work/Eclipse/AculData:"
"../Detectors/Debug:../Detectors:"
"../TELoss/Debug:../TELoss:"
"../Be/Debug:../Be:"
);*/
// h.MakeIndex("");
/*h.SetInputDir("/media/commondata/work/makefiles/AculData/Debug:/media/commondata/work/Eclipse/AculData:"
"../Detectors/Debug:../Detectors:"
"../TELoss/Debug:../TELoss:"
"../Be/Debug:../Be:"
);*/
// "$ROOTSYS/include/root");
// h.MakeAll();
// h.MakeClass("AculRaw");
// h.MakeClass("AculConvert");
// h.MakeClass("AculCalibration");
//
// h.MakeClass("AnnularDetector");
//
// h.MakeClass("TELoss");
//
// h.MakeClass("TObject");
// h.MakeClass("TClass");
}
{
gSystem->Load("/home/chudoba/work/makefilesBe/libAculData.so");
gSystem->Load("/home/chudoba/work/makefilesBe/libDetectors.so");
gSystem->Load("/home/chudoba/work/makefilesBe/libTELoss.so");
gSystem->Load("/home/chudoba/work/makefilesBe/libBe.so");
}
{
// gSystem->Load("/media/commondata/work/makefiles/libAculData.so");
gSystem->Load("/home/w2/work/makefiles/libTELoss.so");
// gSystem->Load("/media/commondata/work/makefiles/libAculData.so");
}
################################################################################
# Makefile building all dynamic libraries:
# libTELoss.so, libEvent.so, libDetectors.so, libSimulations.so
#
# OBJS - files type *.o (subdir.mk)
# LIBS - list of the loaded libraries (objects.mk)
# RM - command "rm -rf"
# CPP_DEPS - files type *.d (subdir.mk)
# .PHONY -
# .SECONDARY -
#
# Vratislav Chudoba
################################################################################
RM := rm -rf
CC := g++
F90 := gfortran
ROOTLIBS = $(shell root-config --libdir)
ROOTINCS = $(shell root-config --incdir)
PWD = $(shell pwd)
INSTALLFOLDER = $(HOME)/AculLib
ACULDATA = $(PWD)/AculData
ELOSS = $(PWD)/ELoss
TELOSS = $(PWD)/TELoss
DETECTORS = $(PWD)/Detectors
BE = $(PWD)/Be
THREADS = $(PWD)/Threads
APPS = $(PWD)/apps
-include $(ACULDATA)/AculData.mk
-include $(ELOSS)/ELOSS.mk
-include $(TELOSS)/TELoss.mk
-include $(DETECTORS)/Detectors.mk
-include $(BE)/Be.mk
#-include $(THREADS)/Threads.mk
-include $(APPS)/apps.mk
all: libAculData.so \
libTELoss.so \
libDetectors.so \
libBe.so
#ROOT html documentation, it will be done as a program which will be alsa compiled by this makefile, program will be as a last condition after all of the libraries
htmldoc: libTELoss.so
-$(RM) htmldoc
root -l -q html.cxx
#copy the libraries and includes somewhere into $HOME, copy the ./htmldoc folder into some location in $HOME
install:
-mkdir $(INSTALLFOLDER)
-mkdir $(INSTALLFOLDER)/lib
-mkdir $(INSTALLFOLDER)/include
cp lib* $(INSTALLFOLDER)/lib
cp $(ACULDATA)/*.h $(INSTALLFOLDER)/include
cp $(TELOSS)/*.h $(INSTALLFOLDER)/include
cp $(DETECTORS)/*.h $(INSTALLFOLDER)/include
cp -r ./htmldoc $(INSTALLFOLDER)/htmldoc
clean:
-$(RM) $(ACULDATAOBJS) $(ACULDATACPP_DEPS)
-$(RM) $(ACULDATA)/AculDataCint.* libAculData.so
-@echo ' '
-$(RM) $(TELOSSOBJS) $(TELOSSCPP_DEPS)
-$(RM) $(TELOSS)/TELossCint.* libTELoss.so
-@echo ' '
-$(RM) $(DETECTORSOBJS) $(DETECTORSCPP_DEPS)
-$(RM) $(DETECTORS)/DetectorsCint.* libDetectors.so
-@echo ' '
-$(RM) $(BEOBJS) $(BECPP_DEPS)
-$(RM) $(BE)/BeCint.* libBe.so
-@echo ' '
-$(RM) $(THREADSOBJS)
-$(RM) $(THREADS)/ThreadsCint.* libThreads.so
-@echo ' '
-$(RM) htmldoc
-@echo ' '
# Those *Cint* files below need special treating:
$(ACULDATA)/AculDataCint.cpp:
-@echo 'Pre-building AculDataCint.cpp and AculDataCint.h files'
-rootcint -f $(ACULDATA)/AculDataCint.cpp -c -p $(ACULDATA_HEADERS)
-@echo ' '
$(TELOSS)/TELossCint.cpp:
-@echo 'Pre-building TELossCint.cpp and TELossCint.h files'
-rootcint -f $(TELOSS)/TELossCint.cpp -c -p $(TELOSS)/TELoss.h $(TELOSS)/linkdef.h
-@echo ' '
$(DETECTORS)/DetectorsCint.cpp:
-@echo 'Pre-building DetectorsCint.cpp and DetectorsCint.h files'
-rootcint -f $(DETECTORS)/DetectorsCint.cpp -c -p $(DETECTORS_HEADERS)
-@echo ' '
$(BE)/BeCint.cpp:
-@echo 'Pre-building BeCint.cpp and BeCint.h files'
-rootcint -f $(BE)/BeCint.cpp -c -p $(BE_HEADERS)
-@echo ' '
#*.so files
libAculData.so: $(ACULDATAOBJS)
@echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker'
$(CC) -L $(ROOTLIBS) -shared -o"libAculData.so" $(ACULDATAOBJS) $(ACULDATALIBS)
@echo 'Finished building target: $@'
@echo ' '
libTELoss.so: $(TELOSSOBJS)
@echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker'
$(CC) -L . -L $(ROOTLIBS) -shared -o"libTELoss.so" $(TELOSSOBJS) $(TELOSSLIBS)
@echo 'Finished building target: $@'
@echo ' '
libDetectors.so: $(DETECTORSOBJS)
@echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker'
g++ -L $(ROOTLIBS) -shared -o"libDetectors.so" $(DETECTORSOBJS) $(DETECTORSLIBS)
@echo 'Finished building target: $@'
@echo ' '
libBe.so: libAculData.so libDetectors.so libTELoss.so $(BEOBJS)
@echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker'
g++ -L . -L $(ROOTLIBS) -shared -o"libBe.so" $(BEOBJS) $(BELIBS)
@echo 'Finished building target: $@'
@echo ' '
fillExpFile: $(APPSOBJS)
@echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker'
g++ -L . -L $(ROOTLIBS) -o"fillExpFile" $(APPSOBJS) $(APPSLIBS)
@echo 'Finished building target: $@'
@echo ' '
.PHONY: all clean
#.SECONDARY: AculData_pre-build TELoss_pre-build Detectors_pre-build libAculData.so libTELoss.so libDetectors.so
# Each subdirectory must supply rules for building sources it contributes
%.o: %.cpp
@echo 'Building file: $@'
@echo 'Invoking: GCC C++ Compiler'
# $(CC) -I$(ROOTINCS) -O0 -g3 -Wall -c -fmessage-length=0 -fPIC -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o"$@" "$<"
# $(CC) -I$(ROOTINCS) -DTHREADING -I$(BOOSTINCS) -O2 -Wall -mmmx -msse -msse2 -msse3 -mfpmath=sse,387 -march=nocona -mtune=nocona -c -fmessage-length=0 -fPIC -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o"$@" "$<"
$(CC) -I$(ROOTINCS) -O2 -Wall -mmmx -msse -msse2 -msse3 -mfpmath=sse,387 -march=nocona -mtune=nocona -c -fmessage-length=0 -fPIC -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o"$@" "$<"
@echo 'Finished building: $@'
@echo ' '
# fortran object files
$(TELOSS)/ELOSS.o:
@echo 'Building file: $@'
@echo 'Invoking: gfortran'
$(F90) -c -fPIC -o"$(TELOSS)/ELOSS.o" $(TELOSS)/ELOSS.f90
@echo 'Finished building target: $@'
@echo ' '
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment