LCFIVertex  0.7.2
jointprob.cpp
1 #include <algo/inc/jointprob.h>
2 #include <inc/track.h>
3 #include <inc/jet.h>
4 #include <util/inc/helixrep.h>
5 #include <util/inc/projection.h>
6 #include <vector>
7 #include <string>
8 #include <map>
9 #include <iostream>
10 #include <util/inc/string.h>
11 
12 namespace vertex_lcfi
13 {
14 
16  {
17  _MaxD0Significance = 200;
18  _MaxD0andZ0 = 5.0;
19 
20  temp.clear();
21  temp.push_back(1.01313412);
22  temp.push_back(0.0246350896);
23  temp.push_back(0.102197811);
24  temp.push_back(0.0411203019);
25  temp.push_back(0.0157710761);
26 
27  _ResolutionParameterRphi = temp;
28 
29  temp1.clear();
30  temp1.push_back(1.01629865);
31  temp1.push_back(0.0271386635);
32  temp1.push_back(0.0948112309);
33  temp1.push_back(0.0410759225);
34  temp1.push_back(0.0148685882);
35 
36  _ResolutionParameterZ = temp1;
37 
38  temp2.clear();
39  temp2.push_back(1.02015948);
40  temp2.push_back(0.0177643765);
41  temp2.push_back(0.144750029);
42  temp2.push_back(0.0288017225);
43  temp2.push_back(0.0237413906);
44 
45  _ResolutionParameter3D = temp2;
46 
47  _ParameterNames.push_back("MaxD0Significance");
48  _ParameterNames.push_back("MaxD0andZ0");
49  _ParameterNames.push_back("ResolutionParameterRphi");
50  _ParameterNames.push_back("ResolutionParameterZ");
51  //_ParameterNames.push_back("ResolutionParameter3D");
52  _ParameterValues.push_back(makeString(_MaxD0Significance));
53  _ParameterValues.push_back(makeString(_MaxD0andZ0));
54  _ParameterValues.push_back("Default Values");
55  _ParameterValues.push_back("Default Values");
56  _ParameterValues.push_back("Default Values");
57 
58  }
59 
60  string JointProb::name() const
61  {
62  return _Name;
63  }
64 
65  std::vector<string> JointProb::parameterNames() const
66  {
67  return _ParameterNames;
68  }
69 
70  std::vector<string> JointProb::parameterValues() const
71  {
72  _ParameterValues.clear();
73  _ParameterValues.push_back(makeString(_MaxD0Significance));
74  _ParameterValues.push_back(makeString(_MaxD0andZ0));
75  _ParameterValues.push_back(makeString(_ResolutionParameterRphi[0]));
76  _ParameterValues.push_back(makeString(_ResolutionParameterRphi[1]));
77  _ParameterValues.push_back(makeString(_ResolutionParameterRphi[2]));
78  _ParameterValues.push_back(makeString(_ResolutionParameterRphi[3]));
79  _ParameterValues.push_back(makeString(_ResolutionParameterRphi[4]));
80  _ParameterValues.push_back(makeString(_ResolutionParameterZ[0]));
81  _ParameterValues.push_back(makeString(_ResolutionParameterZ[1]));
82  _ParameterValues.push_back(makeString(_ResolutionParameterZ[2]));
83  _ParameterValues.push_back(makeString(_ResolutionParameterZ[3]));
84  _ParameterValues.push_back(makeString(_ResolutionParameterZ[4]));
85  //_ParameterValues.push_back(makeString(_ResolutionParameter3D));
86  return _ParameterValues;
87  }
88 
89  void JointProb::setStringParameter(const string & Parameter, const string & Value)
90  {
91  this->badParameter(Parameter);
92  }
93 
94  void JointProb::setDoubleParameter(const string & Parameter, const double Value)
95  {
96  if (Parameter == "MaxD0Significance")
97  {
98  _MaxD0Significance = Value;
99  return;
100  }
101  if (Parameter == "MaxD0andZ0")
102  {
103  _MaxD0andZ0 = Value;
104  return;
105  }
106  this->badParameter(Parameter);
107  }
108 
109  void JointProb::setPointerParameter(const string & Parameter, void * Value)
110  {
111  if (Parameter == "ResolutionParameterRphi")
112  {
113  _ResolutionParameterRphi = *(std::vector<double>*) Value;
114  return;
115  }
116  if (Parameter == "ResolutionParameterZ")
117  {
118  _ResolutionParameterZ = *(std::vector<double>*) Value;
119  return;
120  }
121  if (Parameter == "ResolutionParameter3D")
122  {
123  _ResolutionParameter3D = *(std::vector<double>*) Value;
124  return;
125  }
126  else this->badParameter(Parameter);
127  }
128 
129  std::map<Projection ,double> JointProb::calculateFor(Jet* MyJet) const
130  {
131 
132  std::map<Projection,double> ResultMap;
133 
134 
135  double totalprod[3] = {1,1,1};
136  double jprob[3] = {0,0,0};
137  int ntraks[3] = {0,0,0};
138  int j= 0;
139  int k= 0;
140  double maxdz = _MaxD0andZ0; // maximum cuts for d0 and z0 implemented below
141  double maxd0sig = _MaxD0Significance; // maximum cuts for d0 significance implemented below
142  double significancecompare = 0;
143  float sigma = 0 ;
144  float fact =0;
145 
146 
147 
148  for (std::vector<Track*>::const_iterator iTrack= (MyJet->tracks().begin()); iTrack != (MyJet->tracks().end()) ;++iTrack)
149  {
150  //series of cuts
151 
152  if( fabs( (*iTrack)->helixRep().d0() )< maxdz && fabs( (*iTrack)->helixRep().z0() ) < maxdz )
153  {
154 
155  //(d0,z0,3-d)
156  for(j = 0; j<3; j++ )
157  {
158  switch (j)
159  {
160  case 0:
161  {
162  significancecompare = fabs((*iTrack)->significance(RPhi));
163  break;
164  }
165  case 1:
166  {
167  significancecompare = fabs((*iTrack)->significance( Z/*,MyJet->momentum()*/ ));
168 
169  break;
170  }
171  case 2:
172  {
173  significancecompare = fabs((*iTrack)->significance(ThreeD));
174  break;
175  }
176  }
177 
178  if( significancecompare < maxd0sig )
179  {
180 
181  //calculate vertex parameter probability and multiply them, notice the normalization function
182  totalprod[j] *= probparam( significancecompare ,j , maxd0sig )/ probparam( 0 , j, maxd0sig );
183 
184 
185  ntraks[j]++;
186  }
187  }
188  }
189  }
190 
191 
192 
193  for( j = 0; j<3; j++ )
194  {
195 
196  if( ntraks[j]>0 )
197  {
198  fact =1;
199  sigma = 0;
200 
201  for( k =0; k < ntraks[j]; k++ )
202  {
203 
204  if( k > 0 ) fact *= k;
205  else fact = 1 ;
206 
207  //the summation over the tracks
208 
209  sigma += pow( ( -log( totalprod[j] ) ), k ) / fact;
210 
211  }
212  //and here we combine everything in the final probability
213  jprob[j] = totalprod[j] * sigma;
214  }
215 
216  else jprob [j] = -1; //default no chance value
217  }
218 
219  ResultMap[RPhi] = jprob[0];
220  ResultMap[Z] = jprob[1];
221  ResultMap[ThreeD] = jprob[2];
222 
223  return ResultMap;
224  }
225 
226  double JointProb::probparam(double parameter, int thecoord, double maxd0sig ) const
227  {
228 
229  double prob = 0;
230  double resolutionparameters[3][5];
231  if (_ResolutionParameterRphi.size() != 5 || _ResolutionParameterZ.size() != 5 || _ResolutionParameter3D.size() != 5)
232  std::cerr << "Warning jointprob.cpp:229 Parameters of wrong length" << std::endl;
233  for( int iii=0; iii<5; iii++)
234  {
235  resolutionparameters[0][iii] = _ResolutionParameterRphi[iii];
236  resolutionparameters[1][iii] = _ResolutionParameterZ[iii];
237  resolutionparameters[2][iii] = _ResolutionParameter3D[iii];
238  }
239 
240  // The if statement takes into account a different parametrization for different coordinates.
241  if ( thecoord < 2 )
242  {
243  // part one is the gaussian part
244  // to understand this part better one should look at the meaning of the complementary error function
245  prob = erfc( parameter / (sqrt( double(2) ) * resolutionparameters[thecoord][0] ))
246  -erfc( maxd0sig / ( sqrt( double(2) ) * resolutionparameters[thecoord][0] ) );
247 
248  // part 2 is the added exponential tails
249 
250  prob += resolutionparameters[thecoord][1]* ( exp(- resolutionparameters[thecoord][2] * parameter )
251  -exp(- resolutionparameters[thecoord][2] * maxd0sig ) )
252  +resolutionparameters[thecoord][3]* ( exp(- resolutionparameters[thecoord][4] * parameter )
253  -exp(- resolutionparameters[thecoord][4] * maxd0sig ) );
254 
255 
256  }
257  else
258  {
259 
260  //in this view the gaussian part is just squared.
261 
262  prob = exp(- ( parameter * parameter ) / ( resolutionparameters[thecoord][0] * resolutionparameters[thecoord][0] * double ( 2 ) ) )
263  -exp(- ( maxd0sig * maxd0sig ) / ( resolutionparameters[thecoord][0] * resolutionparameters[thecoord][0] * double ( 2 ) ) );
264 
265  prob += resolutionparameters[thecoord][1]* ( ( 1 + resolutionparameters[thecoord][2] * parameter )
266  *exp ( - resolutionparameters[thecoord][2] * parameter )
267  - ( 1 + resolutionparameters[thecoord][2] * maxd0sig )
268  *exp ( - resolutionparameters[thecoord][2] * maxd0sig ))
269  +resolutionparameters[thecoord][3]* ( ( 1 + resolutionparameters[thecoord][4] * parameter )
270  *exp ( - resolutionparameters[thecoord][4] * parameter )
271  - ( 1 + resolutionparameters[thecoord][4] * maxd0sig )
272  *exp ( - resolutionparameters[thecoord][4] * maxd0sig ));
273 
274  }
275 
276  return prob;
277 
278  }
279 
280 
281 }
std::map< Projection, double > calculateFor(Jet *MyJet) const
Run the algorithm on a jet.
Definition: jointprob.cpp:129
double prob(double ChiSquared, double DegreesOfFreedom)
Calculate probability from Chi Squared and Degrees of freedom.
void setStringParameter(const string &Parameter, const string &Value)
Set String Parameter.
Definition: jointprob.cpp:89
const std::vector< Track * > & tracks() const
Tracks.
Definition: jet.cpp:23
void setPointerParameter(const string &Parameter, void *Value)
Set Pointer Parameter.
Definition: jointprob.cpp:109
void setDoubleParameter(const string &Parameter, const double Value)
Set Double Parameter.
Definition: jointprob.cpp:94
std::vector< string > parameterValues() const
Parameter Values.
Definition: jointprob.cpp:70
JointProb()
Constructor.
Definition: jointprob.cpp:15
string name() const
Name.
Definition: jointprob.cpp:60
std::vector< string > parameterNames() const
Parameter Names.
Definition: jointprob.cpp:65