LCFIVertex  0.7.2
SignificanceFit.cc
1 // First of all, make sure that AIDA was enabled by setting the environment
2 // // variable MARLIN_USE_AIDA when Marlin was compiled. The makefile will then
3 // // have done the setup and defined this macro.
4 //
5 #ifndef MARLIN_USE_AIDA
6 
7 #warning "--------------------------------------------------------------------------------"
8 #warning "- SignificanceFit requires MARLIN_USE_AIDA to be defined. Did you enable -"
9 #warning "- AIDA when compiling Marlin? SignificanceFit will not be compiled. -"
10 #warning "--------------------------------------------------------------------------------"
11 
12 // Can't do anything else.
13 #else
14 
15 //change USING_RAIDA to USING_JAIDA if you are interested in using this processor! In RAIDA the processor is only a place holder! This is really due to the limited functionality of RAIDA!!!
16 #define USING_RAIDA
17 
18 #ifdef USING_RAIDA
19 #warning "USING_RAIDA defined"
20 #else
21 #define USING_JAIDA
22 #warning "USING_JAIDA defined"
23 #endif
24 
25 
26 
27 #include"SignificanceFit.h"
28 #include <iostream>
29 #include <vector>
30 #include <algorithm>
31 #include <set>
32 #include <string>
33 
34 
35 // LCIO includes
36 #include <EVENT/LCCollection.h>
37 #include "EVENT/LCFloatVec.h"
38 
39 #include <EVENT/Vertex.h>
40 
41 //LCFI includes
42 #include <inc/jet.h>
43 #include <inc/track.h>
44 #include <inc/event.h>
45 #include <util/inc/memorymanager.h>
46 #include <inc/lciointerface.h>
47 #include <util/inc/vector3.h>
48 #include <util/inc/projection.h>
49 
50 
51 // AIDA includes...
52 
53 
54 #include <marlin/AIDAProcessor.h>
55 #include <AIDA/IHistogramFactory.h>
56 #include <AIDA/IHistogram1D.h>
57 #include <AIDA/IAxis.h>
58 #include <AIDA/IFitFactory.h>
59 #include <AIDA/IFitter.h>
60 #include <AIDA/IPlotter.h>
61 #include <AIDA/IPlotterFactory.h>
62 #include <AIDA/IPlotterRegion.h>
63 #include <AIDA/IFitResult.h>
64 #include <AIDA/IAnalysisFactory.h>
65 #include <AIDA/IFunctionFactory.h>
66 #include <AIDA/IFunction.h>
67 #include <AIDA/IFitParameterSettings.h>
68 
69 
70 #include <UTIL/LCRelationNavigator.h>
71 
72 
73 #define PI 3.14159265
74 
75 using namespace lcio ;
76 using namespace marlin ;
77 using namespace vertex_lcfi;
78 
79 SignificanceFit aSignificanceFit;
80 
81 
82 SignificanceFit::SignificanceFit() : Processor("SignificanceFit")
83 {
84 
85  registerInputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
86  "JetRPCollection" ,
87  "Name of the ReconstructedParticle collection that represents jets",
88  _JetRPColName ,
89  std::string("Jets") ) ;
90  registerInputCollection( lcio::LCIO::VERTEX,
91  "IPVertexCollection" ,
92  "Name of the Vertex collection that contains the primary vertex (Optional)" ,
93  _IPVertexCollectionName ,
94  std::string("IPVertex") ) ;
95  registerOptionalParameter( "Cutoffgauss" ,
96  "Maximum value for the local fitting of the gaussian, starting point 0" ,
97  _Cutoffgauss,
98  double(5));
99  registerOptionalParameter( "Cutoffexp1" ,
100  "Maximum value for the local fitting of the first exponential, starting point Cutoffgauss" ,
101  _Cutoffexp1,
102  double(40));
103 
104  registerOptionalParameter( "Cutoffexp2" ,
105  "Maximum value for the local fitting of the second exponential, starting point Cutoffexp1",
106  _Cutoffexp2,
107  double(200));
108  registerOptionalParameter( "GaussAmlitudeInit" ,
109  "Initialization value for Gauss amplitude in the Gaussian local fit",
110  _GaussAmp,
111  double(5000));
112  registerOptionalParameter( "GaussSigmaInit" ,
113  "Initialization value for Gauss sigma in the Gaussian local fit",
114  _GaussSigma,
115  double(2));
116 
117  registerOptionalParameter( "ExpAmlitudeInit" ,
118  "Initialization value for amplitude of both exponentials in the exponentials local fit",
119  _ExpAmp,
120  double(100));
121  registerOptionalParameter( "ExpLambdaInit" ,
122  "Initialization value for decay constant of both exponentials in the exponentials local fit",
123  _ExpLambda,
124  double(0));
125 }
126 
127 
128 void SignificanceFit::init()
129 {
130 
131 #ifdef USING_JAIDA
132 
133  AIDA::IHistogramFactory* histofact = AIDAProcessor::histogramFactory(this);
134 
135  rphisto =histofact->createHistogram1D("-d_0 over sigma(d_0)",400,0,_Cutoffexp2);
136  zhisto = histofact->createHistogram1D("-z_0 over sigma(z_0)",400,0,_Cutoffexp2);
137 
138 
139  // A work around doing local fitting given that AIDAJNI does not yet support this. AIDA 3.3 tough does.
140 
141  rphistogauss =histofact->createHistogram1D("-d_0 over sigma(d_0)",10,0,_Cutoffgauss);
142  zhistogauss =histofact->createHistogram1D("-z_0 over sigma(z_0)",10,0,_Cutoffgauss);
143 
144  rphistoexp1 =histofact->createHistogram1D("-d_0 over sigma(d_0)",10,_Cutoffgauss,_Cutoffexp1);
145  zhistoexp1 =histofact->createHistogram1D("-z_0 over sigma(z_0)",10,_Cutoffgauss,_Cutoffexp1);
146 
147  rphistoexp2 =histofact->createHistogram1D("-d_0 over sigma(d_0)",10,_Cutoffexp1,_Cutoffexp2);
148  zhistoexp2 =histofact->createHistogram1D("-z_0 over sigma(z_0)",10,_Cutoffexp1,_Cutoffexp2);
149 
150  printParameters() ;
151 #endif
152 
153  _nRun = 0 ;
154  _nEvt = 0 ;
155 
156 }
157 
158 
159 void SignificanceFit::processEvent( LCEvent * evt )
160 {
161 #ifdef USING_JAIDA
162 
163  LCCollection* JetRPCol = evt->getCollection( _JetRPColName );
164  LCCollection* VertexCol;
165  VertexCol = evt->getCollection( _IPVertexCollectionName );
166 
167  Vector3 IPPos;
168  Matrix3x3 IPErr;
169 
170  int nVerts = VertexCol->getNumberOfElements() ;
171  bool done = 0;
172 
173  for(int i=0; i< nVerts ; i++)
174  {
175  lcio::Vertex* iVertex = dynamic_cast<lcio::Vertex*>(VertexCol->getElementAt(i));
176  if (iVertex->isPrimary())
177  {
178  IPPos.x() = iVertex->getPosition()[0];
179  IPPos.y() = iVertex->getPosition()[1];
180  IPPos.z() = iVertex->getPosition()[2];
181  IPErr(0,0) = iVertex->getCovMatrix()[0];
182  IPErr(1,0) = iVertex->getCovMatrix()[1];
183  IPErr(1,1) = iVertex->getCovMatrix()[2];
184  IPErr(2,0) = iVertex->getCovMatrix()[3];
185  IPErr(2,1) = iVertex->getCovMatrix()[4];
186  IPErr(2,2) = iVertex->getCovMatrix()[5];
187  done = 1;
188  }
189  if (done) break;
190  }
191  if (!done) done = 1;//TODO Throw something
192 
193 
194  Event* MyEvent = new Event(IPPos,IPErr);
195 
196  //Event* MyEvent = new Event(Vector3(),SymMatrix3x3());
197 
198  MemoryManager<Event>::Event()->registerObject(MyEvent);
199  std::vector< vertex_lcfi::Jet*> ThisVector;
200 
201  int nRCP = JetRPCol->getNumberOfElements() ;
202 
203  if(nRCP ==0 ) std::cerr<<"Warning: SignificanceFit:336 : NO jets present "<<std::endl;
204 
205  for(int i=0; i< nRCP ; i++)
206  {
207  ReconstructedParticle* JetRP = dynamic_cast<ReconstructedParticle*>(JetRPCol->getElementAt(i));
208  vertex_lcfi::Jet* ThisJet = jetFromLCIORP(MyEvent,JetRP);
209  ThisVector.push_back(ThisJet);
210  }
211 
212  if(ThisVector.size()>0)
213  {
214  for(std::vector<Jet*>::iterator iJet2=ThisVector.begin();iJet2 != ThisVector.end();++iJet2)
215  {
216 
217  for (std::vector<vertex_lcfi::Track*>::const_iterator iTrack = (*iJet2)->tracks().begin(); iTrack != (*iJet2)->tracks().end() ;++iTrack)
218  {
219 
220  double SignificanceRPhi = (*iTrack)->signedSignificance(RPhi, *iJet2);
221  double SignificanceZ = (*iTrack)->signedSignificance(Z, *iJet2);
222  if(SignificanceRPhi<0)
223  {
224 
225  rphistogauss->fill(fabs(SignificanceRPhi));
226  rphistoexp1->fill(fabs(SignificanceRPhi));
227  rphistoexp2->fill(fabs(SignificanceRPhi));
228  rphisto->fill(fabs(SignificanceRPhi));
229 
230  }
231  if(SignificanceZ<0)
232  {
233  zhistogauss->fill(fabs(SignificanceZ));
234  zhistoexp1->fill(fabs(SignificanceZ));
235  zhistoexp2->fill(fabs(SignificanceZ));
236  zhisto->fill(fabs(SignificanceZ));
237 
238  }
239  }
240 
241  }
242 
243  }
244  MetaMemoryManager::Event()->delAllObjects();
245 #endif
246 
247  _nEvt++;
248 
249 }
250 
251 
252 void SignificanceFit::processRunHeader( LCRunHeader* run)
253 {
254  _nRun++ ;
255 }
256 
257 void SignificanceFit::check( LCEvent * evt )
258 {
259 
260 
261 
262 }
263 
264 
265 
266 void SignificanceFit::end(){
267 
268 #ifdef USING_JAIDA
269 
270  _analysisFactory = AIDAProcessor::GetIAnalysisFactory(this);
271  AIDA::IFitFactory* fitFactory =_analysisFactory->createFitFactory();
272  AIDA::IPlotter* plotter = _analysisFactory->createPlotterFactory()->create("Plot");
273  AIDA::IPlotter* plotterZ = _analysisFactory->createPlotterFactory()->create("Plot");
274  AIDA::IFitter* fitter = fitFactory->createFitter("Chi2");
275  AIDA::IFunctionFactory* funFactory =_analysisFactory->createFunctionFactory(*AIDAProcessor::tree(this));
276 
277  std::vector<double> initialParsgauss;
278  std::vector<double> initialPars;
279  std::vector<double> globalPars;
280  std::vector<double> outPars;
281 
282  std::vector<double> zinitialParsgauss;
283  std::vector<double> zinitialPars;
284  std::vector<double> zglobalPars;
285  std::vector<double> zoutPars;
286 
287  AIDA::IFunction* functiongauss = funFactory->createFunctionByName("g","g");
288 
289  AIDA::IFunction* functionexp1 = funFactory->createFunctionByName("e","e");
290  AIDA::IFunction* functionexp2 = funFactory->createFunctionByName("e","e");
291 
292  AIDA::IFunction* functionfinal = funFactory->createFunctionByName("g","g+e+e");
293 
294 
295  initialPars.push_back(_ExpAmp);
296  initialPars.push_back(_ExpLambda);
297 
298 
299  initialParsgauss.push_back(_GaussAmp);
300  initialParsgauss.push_back(0);
301  initialParsgauss.push_back(_GaussSigma);
302 
303  functiongauss->setParameters( initialParsgauss );
304 
305  functionexp1->setParameters( initialPars );
306  functionexp2->setParameters( initialPars );
307 
308 
309  fitter->fitParameterSettings("mean").setBounds(-0.0000000001,0.0000000001);
310 
311  AIDA::IFitResult* resultgauss = fitter->fit(*rphistogauss, *functiongauss);
312  AIDA::IFitResult* resultexp1 = fitter->fit(*rphistoexp1, *functionexp1);
313  AIDA::IFitResult* resultexp2 = fitter->fit(*rphistoexp2, *functionexp2);
314 
315  std::vector<std::string> fParNamesgauss = resultgauss->fittedParameterNames();
316  std::vector<std::string> fParNamesexp1 = resultexp1->fittedParameterNames();
317  std::vector<std::string> fParNamesexp2 = resultexp2->fittedParameterNames();
318  std::vector<double> fParsgauss = resultgauss->fittedParameters();
319  std::vector<double> fParsexp1 = resultexp1->fittedParameters();
320  std::vector<double> fParsexp2 = resultexp2->fittedParameters();
321 
322  std::cout<<" RPHI - Plane "<<std::endl;
323  std::cout<<std::endl;
324 
325  std::cout<<" Local Fit of Gaussian "<<std::endl;
326  for(unsigned int i=0; i<fParsgauss.size();i++)
327  {
328  std::cout<<" Parameter "<< fParNamesgauss[i]<<" "<<fParsgauss[i]<<std::endl;
329  globalPars.push_back(fParsgauss[i]);
330  }
331  std::cout<<"Chi^2/ndf: "<<resultgauss->quality()<<std::endl;
332 
333  std::cout<<" Local Fit of first Exponential "<<std::endl;
334  for(unsigned int i=0; i<fParsexp1.size();i++)
335  {
336  std::cout<<" Parameter "<< fParNamesexp1[i]<<" "<<fParsexp1[i]<<std::endl;
337  globalPars.push_back(fParsexp1[i]);
338  }
339  std::cout<<"Chi^2/ndf: "<<resultexp1->quality()<<std::endl;
340 
341  std::cout<<" Local Fit of second Exponential "<<std::endl;
342  for(unsigned int i=0; i<fParsexp2.size();i++)
343  {
344  std::cout<<" Parameter "<< fParNamesexp1[i]<<" "<<fParsexp2[i]<<std::endl;
345  globalPars.push_back(fParsexp2[i]);
346  }
347 
348  std::cout<<"Chi^2/ndf: "<<resultexp2->quality()<<std::endl;
349 
350 
351  functionfinal->setParameters(globalPars);
352  AIDA::IFitResult* resultfinal = fitter->fit(*rphisto, *functionfinal);
353 
354  std::vector<std::string> fParNamesglobal = resultfinal->fittedParameterNames();
355  std::vector<double> fParsglobal = resultfinal->fittedParameters();
356 
357  std::cout<<" Global Fit "<<std::endl;
358  for(unsigned int i=0; i<fParsglobal.size();i++)
359  {
360  std::cout<<" Parameter "<< fParNamesglobal[i]<<" "<<fParsglobal[i]<<std::endl;
361  }
362  std::cout<<"Chi^2/ndf: "<<resultfinal->quality()<<std::endl;
363 
364  plotter->createRegions(2,2,0);
365  plotter->region(0)->plot(*rphistogauss);
366  plotter->region(1)->plot(*rphistoexp1);
367  plotter->region(2)->plot(*rphistoexp2);
368  plotter->region(3)->plot(*rphisto);
369  plotter->region(0)->plot(resultgauss->fittedFunction());
370  plotter->region(1)->plot(resultexp1->fittedFunction());
371  plotter->region(2)->plot(resultexp2->fittedFunction());
372  plotter->region(3)->plot(resultfinal->fittedFunction());
373 
374  outPars.push_back(fParsglobal[2]);
375  outPars.push_back(-sqrt(2/(PI))*fParsglobal[3]/(fParsglobal[0]*fParsglobal[2]*fParsglobal[4]));
376  outPars.push_back(-(fParsglobal[4]));
377  outPars.push_back(-sqrt(2/(PI))*fParsglobal[5]/(fParsglobal[0]*fParsglobal[2]*fParsglobal[6]));
378  outPars.push_back(-(fParsglobal[6]));
379 
380  plotter->show();
381 
382 
383  zinitialPars.push_back(100);
384  zinitialPars.push_back(0);
385 
386  zinitialParsgauss.push_back(5000);
387  zinitialParsgauss.push_back(0);
388  zinitialParsgauss.push_back(2);
389 
390  functiongauss->setParameters( zinitialParsgauss );
391 
392  functionexp1->setParameters( zinitialPars );
393  functionexp2->setParameters( zinitialPars );
394 
395 
396  fitter->fitParameterSettings("mean").setBounds(-0.0000000001,0.0000000001);
397 
398  AIDA::IFitResult* zresultgauss = fitter->fit(*zhistogauss, *functiongauss);
399  AIDA::IFitResult* zresultexp1 = fitter->fit(*zhistoexp1, *functionexp1);
400  AIDA::IFitResult* zresultexp2 = fitter->fit(*zhistoexp2, *functionexp2);
401 
402 
403  std::vector<std::string> zfParNamesgauss = zresultgauss->fittedParameterNames();
404  std::vector<std::string> zfParNamesexp1 = zresultexp1->fittedParameterNames();
405  std::vector<std::string> zfParNamesexp2 = zresultexp2->fittedParameterNames();
406  std::vector<double> zfParsgauss = zresultgauss->fittedParameters();
407  std::vector<double> zfParsexp1 = zresultexp1->fittedParameters();
408  std::vector<double> zfParsexp2 = zresultexp2->fittedParameters();
409 
410  std::cout<<std::endl;
411  std::cout<<std::endl;
412  std::cout<<std::endl;
413 
414  std::cout<<" Z - axis "<<std::endl;
415  std::cout<<std::endl;
416 
417  std::cout<<" Local Fit of Gaussian "<<std::endl;
418  for(unsigned int i=0; i<fParsgauss.size();i++)
419  {
420  std::cout<<" Parameter "<< zfParNamesgauss[i]<<" "<<zfParsgauss[i]<<std::endl;
421  zglobalPars.push_back(zfParsgauss[i]);
422  }
423  std::cout<<"Chi^2/ndf: "<<zresultgauss->quality()<<std::endl;
424 
425  std::cout<<" Local Fit of first Exponential "<<std::endl;
426  for(unsigned int i=0; i<zfParsexp1.size();i++)
427  {
428  std::cout<<" Parameter "<< zfParNamesexp1[i]<<" "<<zfParsexp1[i]<<std::endl;
429  zglobalPars.push_back(zfParsexp1[i]);
430  }
431  std::cout<<"Chi^2/ndf: "<<zresultexp1->quality()<<std::endl;
432 
433  std::cout<<" Local Fit of second Exponential "<<std::endl;
434  for(unsigned int i=0; i<zfParsexp2.size();i++)
435  {
436  std::cout<<" Parameter "<< zfParNamesexp1[i]<<" "<<zfParsexp2[i]<<std::endl;
437  globalPars.push_back(zfParsexp2[i]);
438  }
439 
440  std::cout<<"Chi^2/ndf: "<<zresultexp2->quality()<<std::endl;
441 
442 
443  functionfinal->setParameters(zglobalPars);
444  AIDA::IFitResult* zresultfinal = fitter->fit(*zhisto, *functionfinal);
445 
446  std::vector<std::string> zfParNamesglobal = zresultfinal->fittedParameterNames();
447  std::vector<double> zfParsglobal = zresultfinal->fittedParameters();
448 
449  std::cout<<" Global Fit "<<std::endl;
450  for(unsigned int i=0; i<zfParsglobal.size();i++)
451  {
452  std::cout<<" Parameter "<< zfParNamesglobal[i]<<" "<<zfParsglobal[i]<<std::endl;
453  }
454  std::cout<<"Chi^2/ndf: "<<zresultfinal->quality()<<std::endl;
455 
456  plotterZ->createRegions(2,2,0);
457  plotterZ->region(0)->plot(*zhistogauss);
458  plotterZ->region(1)->plot(*zhistoexp1);
459  plotterZ->region(2)->plot(*zhistoexp2);
460  plotterZ->region(3)->plot(*zhisto);
461  plotterZ->region(0)->plot(zresultgauss->fittedFunction());
462  plotterZ->region(1)->plot(zresultexp1->fittedFunction());
463  plotterZ->region(2)->plot(zresultexp2->fittedFunction());
464  plotterZ->region(3)->plot(zresultfinal->fittedFunction());
465 
466  plotterZ->show();
467 
468  zoutPars.push_back(zfParsglobal[2]);
469  zoutPars.push_back(-sqrt(2/(PI))*zfParsglobal[3]/(zfParsglobal[0]*zfParsglobal[2]*zfParsglobal[4]));
470  zoutPars.push_back(-1*zfParsglobal[4]);
471  zoutPars.push_back(-sqrt(2/(PI))*zfParsglobal[5]/(zfParsglobal[0]*zfParsglobal[2]*zfParsglobal[6]));
472  zoutPars.push_back(-1*zfParsglobal[6]);
473 
474 
475  std::cout<<"PARAMETERS FOR RPHI Joint Probability"<<std::endl;
476  for(unsigned int i=0;i<5;i++)
477  {
478  std::cout<<outPars[i]<<std::endl;
479  }
480  std::cout<<"PARAMETERS FOR Z Joint Probability"<<std::endl;
481  for(unsigned int i=0;i<5;i++)
482  {
483  std::cout<<zoutPars[i]<<std::endl;
484  }
485 
486  MetaMemoryManager::Run()->delAllObjects();
487 
488 #endif
489  std::cout << "SignificanceFitProcessor::end() " << name()
490  << " processed " << _nEvt << " events in " << _nRun << " runs "
491  << std::endl ;
492 
493 }
494 
495 
496 
497 #endif // end if ifdef MARLIN_USE_AIDA
498 
499 
vertex_lcfi::MemoryManager::registerObject
void registerObject(T *pointer)
Register an object for memory management.
Definition: include/vertex_lcfi/util/inc/memorymanager.h:109
vertex_lcfi::Jet
Simple Jet Class.
Definition: include/vertex_lcfi/inc/jet.h:24
vertex_lcfi::Event
Event.
Definition: include/vertex_lcfi/inc/event.h:23
vertex_lcfi::MemoryManager
Memory management.
Definition: include/vertex_lcfi/util/inc/memorymanager.h:58
Generated on Fri Dec 2 2016 12:36:49 for LCFIVertex by   doxygen 1.8.6