Newer
Older
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <getopt.h>
#include <fstream>
#include <sstream>
#include <string>
#include "Examples/Tutorial/LineShapes/BwShape.hh"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TMath.h"
#include "ErrLogger/ErrLogger.hh"
BwShape::BwShape(double MassRes, double MassWidth, double MassDec1, double MassDec2, unsigned int Lmax, double deltaMass) :
_theTFile(0)
{
if (MassRes <= MassDec1+MassDec2){
Alert << "Mass of the resonance smaller than the masses of the decay particles\n"
<< MassRes << " < " << MassDec1+MassDec2 << endmsg;
exit(0);
}
std::stringstream Lmaxstrstr;
Lmaxstrstr << Lmax;
std::string rootFileName="./BwShapeLmax"+Lmaxstrstr.str()+".root";
_theTFile=new TFile(rootFileName.c_str(),"recreate");
for (unsigned int i=0; i<=Lmax; ++i){
std::stringstream Lstrstr;
Lstrstr << i;
std::string histName="BreitWigner_L"+Lstrstr.str();
_histMap[i]= new TH1F(histName.c_str(),histName.c_str(),301, MassRes-deltaMass, MassRes+deltaMass);
histName="Argand_L"+Lstrstr.str();
TH2F* currentArgandHist=new TH2F(histName.c_str(),histName.c_str(),301, -1., 1., 301, 0., 1.3);
currentArgandHist->SetXTitle("Re(Bw)");
currentArgandHist->SetYTitle("Im(Bw)");
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
_argandHistMap[i]=currentArgandHist;
//initialize here the phase shift histogramms _phaseHistMap
}
double stepSize=2.*deltaMass/300;
for (unsigned int lIt=0; lIt<=Lmax; ++lIt){
TH1F* currentHist=_histMap[lIt];
TH2F* currentArgandHist=_argandHistMap[lIt];
for (double massIt=MassRes-deltaMass; massIt<MassRes+deltaMass; massIt+=stepSize){
Vector4<double> res4V(massIt, 0., 0., 0.);
complex<double> currentBW=BreitWignerBlattW(res4V, MassDec1, MassDec2, MassRes, MassWidth, lIt);
double weight=norm(currentBW);
currentHist->Fill(massIt,weight);
currentArgandHist->Fill(currentBW.real(),currentBW.imag());
//fill here the phase shift histogramms: hint "atan2(imag,real)"
}
}
}
BwShape::~BwShape()
{
_theTFile->Write();
_theTFile->Close();
}