Jpp  17.2.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JTwist.cc File Reference
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TRandom3.h"
#include "JLang/JException.hh"
#include "JMath/JConstants.hh"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JGeometry3D/JEigen3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JCompass/JHit.hh"
#include "JCompass/JModel.hh"
#include "JFit/JSimplex.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 Test program for fit of quaternions. More...
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Test program for fit of quaternions.

Definition at line 27 of file JTwist.cc.

28 {
29  using namespace std;
30  using namespace JPP;
31 
32  string outputFile;
33  int numberOfEvents;
34  double X;
35  double Y;
36  double Z;
37  double alpha;
38  JPosition3D sigma;
39  UInt_t seed;
40  bool fit;
41  int debug;
42 
43  try {
44 
45  JParser<> zap("Test program for fit of quaternions.");
46 
47  zap['o'] = make_field(outputFile) = "twist.root";
48  zap['n'] = make_field(numberOfEvents) = 1000;
49  zap['x'] = make_field(X, "tilt angle around x-axis [deg]") = 0.0;
50  zap['y'] = make_field(Y, "tilt angle around y-axis [deg]") = 0.0;
51  zap['z'] = make_field(Z, "tilt angle around z-axis [deg]") = 0.0;
52  zap['a'] = make_field(alpha, "twist angle [deg/m]");
53  zap['s'] = make_field(sigma, "(x,y,z) resolution [deg]");
54  zap['S'] = make_field(seed) = 0;
55  zap['F'] = make_field(fit, "enable fit");
56  zap['d'] = make_field(debug) = 2;
57 
58  zap(argc, argv);
59  }
60  catch(const exception &error) {
61  FATAL(error.what() << endl);
62  }
63 
64 
65  gRandom->SetSeed(seed);
66 
67 
68  TFile out(outputFile.c_str(), "recreate");
69 
70  TH1D hc("chi2/NDF", NULL, 100, 0.0, 2.0);
71  TH1D h0("tilt", NULL, 100, 0.0, 2.0);
72  TH1D h1("twist", NULL, 100, 0.0, 0.1);
73 
74 
75  // model
76 
77  const JModel model(JQuaternion3D(JQuaternion3X(X * PI / 180.0),
78  JQuaternion3Y(Y * PI / 180.0),
79  JQuaternion3Z(Z * PI / 180.0)),
80  JQuaternion3D(JQuaternion3Z(alpha * PI / 180.0)));
81 
82  STATUS("Q0 " << FIXED(12,8) << model.Q0 << endl);
83  STATUS("Q1 " << FIXED(12,8) << model.Q1 << endl);
84 
85 
86  JSimplex<JModel> simplex;
87 
88  JSimplex<JModel>::MAXIMUM_ITERATIONS = 10000;
89 
90  simplex.debug = debug;
91 
92  const JChi2 getChi2(EM_NORMAL);
93 
94 
95  for (int counter = 0; counter != numberOfEvents; ++counter) {
96 
97  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
98 
99  // generate data
100 
101  const int N = 18; // number of floors
102  const double z1 = 30.0; // height of first floor
103  const double dz = 10.0; // distance between floors
104  const double eps = 1.0e-6; // minimal sigma
105 
106  vector<JHit> data;
107 
108  for (int i = 1; i <= N; ++i) {
109 
110  const double z = z1 + (i - 1) * dz; // height of floor i
111 
112  data.push_back(JHit(i, z, model(z), max(eps, sigma.getLength())));
113  }
114 
115  // smearing
116 
117  for (vector<JHit>::iterator i = data.begin(); i != data.end(); ++i) {
118  i->mul(JQuaternion3D(JQuaternion3X(gRandom->Gaus(0.0, sigma.getX()) * PI / 180.0),
119  JQuaternion3Y(gRandom->Gaus(0.0, sigma.getY()) * PI / 180.0),
120  JQuaternion3Z(gRandom->Gaus(0.0, sigma.getZ()) * PI / 180.0)));
121  }
122 
123  for (vector<JHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
124  DEBUG(FIXED(4,1) << i->getZ() << ' ' << FIXED(12,8) << i->getQuaternion() << endl);
125  }
126 
127 
128  JModel result(data.begin(), data.end()); // prefit
129 
130  DEBUG("R0 " << FIXED(12,8) << result.Q0 << ' ' << getAngle(model.Q0, result.Q0) << endl);
131  DEBUG("R1 " << FIXED(12,8) << result.Q1 << ' ' << getAngle(model.Q1, result.Q1) << endl);
132 
133  if (fit) {
134 
135  simplex.value = result; // start value
136 
137  simplex.step.resize(4);
138 
139  simplex.step[0] = JModel(JQuaternion3X(0.50 * PI / 180.0), JQuaternion3D::getIdentity());
140  simplex.step[1] = JModel(JQuaternion3Y(0.50 * PI / 180.0), JQuaternion3D::getIdentity());
141  simplex.step[2] = JModel(JQuaternion3Z(0.50 * PI / 180.0), JQuaternion3D::getIdentity());
142  simplex.step[3] = JModel(JQuaternion3D::getIdentity(), JQuaternion3Z(0.05 * PI / 180.0));
143 
144  const double chi2 = simplex(getChi2, data.begin(), data.end());
145  const int ndf = data.size() * 4 - simplex.step.size();
146 
147  DEBUG("number of iterations " << simplex.numberOfIterations << endl);
148 
149  result = simplex.value; // final value
150 
151  DEBUG("chi2 " << chi2 << '/' << ndf << endl);
152  DEBUG("R0 " << FIXED(12,8) << simplex.value.Q0 << ' ' << getAngle(model.Q0, result.Q0) << endl);
153  DEBUG("R1 " << FIXED(12,8) << simplex.value.Q1 << ' ' << getAngle(model.Q1, result.Q1) << endl);
154 
155  hc.Fill(chi2 / ndf);
156  }
157 
158  //h0.Fill(getAngle(model.Q0, result.Q0));
159  h0.Fill(getAngle(JDirection3D(0.0, 0.0, 1.0).rotate(model .Q0),
160  JDirection3D(0.0, 0.0, 1.0).rotate(result.Q0)));
161  h1.Fill(getAngle(model.Q1, result.Q1));
162  }
163  STATUS(endl);
164 
165  out.Write();
166  out.Close();
167 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1517
JGEOMETRY3D::getAngle
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
Definition: JQuaternion3D.hh:1165
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
N
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
Definition: JShowerPostfit.sh:95
Y
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Definition: detector-XY:fit1d.sh:33
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
GAUSS_HERMITE::sigma
const double sigma[]
Definition: JQuadrature.cc:74
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
JTOOLS::result
return result
Definition: JPolint.hh:764
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
chi2
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
X
no fit printf nominal n $STRING awk v X
Definition: detector-XY:fit.sh:102
Z
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
Definition: JAlignDetector.sh:81
JFIT::JSimplex< JModel >::MAXIMUM_ITERATIONS
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JSimplex.hh:237
JFIT::getChi2
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
debug
int debug
debug level
Definition: archive-put-Detectors.sh:95
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5