5 #ifndef MARLIN_USE_AIDA
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 "--------------------------------------------------------------------------------"
19 #warning "USING_RAIDA defined"
22 #warning "USING_JAIDA defined"
27 #include"SignificanceFit.h"
36 #include <EVENT/LCCollection.h>
37 #include "EVENT/LCFloatVec.h"
39 #include <EVENT/Vertex.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>
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>
70 #include <UTIL/LCRelationNavigator.h>
75 using namespace lcio ;
76 using namespace marlin ;
77 using namespace vertex_lcfi;
82 SignificanceFit::SignificanceFit() : Processor(
"SignificanceFit")
85 registerInputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
87 "Name of the ReconstructedParticle collection that represents jets",
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" ,
99 registerOptionalParameter(
"Cutoffexp1" ,
100 "Maximum value for the local fitting of the first exponential, starting point Cutoffgauss" ,
104 registerOptionalParameter(
"Cutoffexp2" ,
105 "Maximum value for the local fitting of the second exponential, starting point Cutoffexp1",
108 registerOptionalParameter(
"GaussAmlitudeInit" ,
109 "Initialization value for Gauss amplitude in the Gaussian local fit",
112 registerOptionalParameter(
"GaussSigmaInit" ,
113 "Initialization value for Gauss sigma in the Gaussian local fit",
117 registerOptionalParameter(
"ExpAmlitudeInit" ,
118 "Initialization value for amplitude of both exponentials in the exponentials local fit",
121 registerOptionalParameter(
"ExpLambdaInit" ,
122 "Initialization value for decay constant of both exponentials in the exponentials local fit",
128 void SignificanceFit::init()
133 AIDA::IHistogramFactory* histofact = AIDAProcessor::histogramFactory(
this);
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);
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);
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);
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);
159 void SignificanceFit::processEvent( LCEvent * evt )
163 LCCollection* JetRPCol = evt->getCollection( _JetRPColName );
164 LCCollection* VertexCol;
165 VertexCol = evt->getCollection( _IPVertexCollectionName );
170 int nVerts = VertexCol->getNumberOfElements() ;
173 for(
int i=0; i< nVerts ; i++)
175 lcio::Vertex* iVertex =
dynamic_cast<lcio::Vertex*
>(VertexCol->getElementAt(i));
176 if (iVertex->isPrimary())
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];
199 std::vector< vertex_lcfi::Jet*> ThisVector;
201 int nRCP = JetRPCol->getNumberOfElements() ;
203 if(nRCP ==0 ) std::cerr<<
"Warning: SignificanceFit:336 : NO jets present "<<std::endl;
205 for(
int i=0; i< nRCP ; i++)
207 ReconstructedParticle* JetRP =
dynamic_cast<ReconstructedParticle*
>(JetRPCol->getElementAt(i));
209 ThisVector.push_back(ThisJet);
212 if(ThisVector.size()>0)
214 for(std::vector<Jet*>::iterator iJet2=ThisVector.begin();iJet2 != ThisVector.end();++iJet2)
217 for (std::vector<vertex_lcfi::Track*>::const_iterator iTrack = (*iJet2)->tracks().begin(); iTrack != (*iJet2)->tracks().end() ;++iTrack)
220 double SignificanceRPhi = (*iTrack)->signedSignificance(RPhi, *iJet2);
221 double SignificanceZ = (*iTrack)->signedSignificance(Z, *iJet2);
222 if(SignificanceRPhi<0)
225 rphistogauss->fill(fabs(SignificanceRPhi));
226 rphistoexp1->fill(fabs(SignificanceRPhi));
227 rphistoexp2->fill(fabs(SignificanceRPhi));
228 rphisto->fill(fabs(SignificanceRPhi));
233 zhistogauss->fill(fabs(SignificanceZ));
234 zhistoexp1->fill(fabs(SignificanceZ));
235 zhistoexp2->fill(fabs(SignificanceZ));
236 zhisto->fill(fabs(SignificanceZ));
244 MetaMemoryManager::Event()->delAllObjects();
252 void SignificanceFit::processRunHeader( LCRunHeader* run)
257 void SignificanceFit::check( LCEvent * evt )
266 void SignificanceFit::end(){
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));
277 std::vector<double> initialParsgauss;
278 std::vector<double> initialPars;
279 std::vector<double> globalPars;
280 std::vector<double> outPars;
282 std::vector<double> zinitialParsgauss;
283 std::vector<double> zinitialPars;
284 std::vector<double> zglobalPars;
285 std::vector<double> zoutPars;
287 AIDA::IFunction* functiongauss = funFactory->createFunctionByName(
"g",
"g");
289 AIDA::IFunction* functionexp1 = funFactory->createFunctionByName(
"e",
"e");
290 AIDA::IFunction* functionexp2 = funFactory->createFunctionByName(
"e",
"e");
292 AIDA::IFunction* functionfinal = funFactory->createFunctionByName(
"g",
"g+e+e");
295 initialPars.push_back(_ExpAmp);
296 initialPars.push_back(_ExpLambda);
299 initialParsgauss.push_back(_GaussAmp);
300 initialParsgauss.push_back(0);
301 initialParsgauss.push_back(_GaussSigma);
303 functiongauss->setParameters( initialParsgauss );
305 functionexp1->setParameters( initialPars );
306 functionexp2->setParameters( initialPars );
309 fitter->fitParameterSettings(
"mean").setBounds(-0.0000000001,0.0000000001);
311 AIDA::IFitResult* resultgauss = fitter->fit(*rphistogauss, *functiongauss);
312 AIDA::IFitResult* resultexp1 = fitter->fit(*rphistoexp1, *functionexp1);
313 AIDA::IFitResult* resultexp2 = fitter->fit(*rphistoexp2, *functionexp2);
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();
322 std::cout<<
" RPHI - Plane "<<std::endl;
323 std::cout<<std::endl;
325 std::cout<<
" Local Fit of Gaussian "<<std::endl;
326 for(
unsigned int i=0; i<fParsgauss.size();i++)
328 std::cout<<
" Parameter "<< fParNamesgauss[i]<<
" "<<fParsgauss[i]<<std::endl;
329 globalPars.push_back(fParsgauss[i]);
331 std::cout<<
"Chi^2/ndf: "<<resultgauss->quality()<<std::endl;
333 std::cout<<
" Local Fit of first Exponential "<<std::endl;
334 for(
unsigned int i=0; i<fParsexp1.size();i++)
336 std::cout<<
" Parameter "<< fParNamesexp1[i]<<
" "<<fParsexp1[i]<<std::endl;
337 globalPars.push_back(fParsexp1[i]);
339 std::cout<<
"Chi^2/ndf: "<<resultexp1->quality()<<std::endl;
341 std::cout<<
" Local Fit of second Exponential "<<std::endl;
342 for(
unsigned int i=0; i<fParsexp2.size();i++)
344 std::cout<<
" Parameter "<< fParNamesexp1[i]<<
" "<<fParsexp2[i]<<std::endl;
345 globalPars.push_back(fParsexp2[i]);
348 std::cout<<
"Chi^2/ndf: "<<resultexp2->quality()<<std::endl;
351 functionfinal->setParameters(globalPars);
352 AIDA::IFitResult* resultfinal = fitter->fit(*rphisto, *functionfinal);
354 std::vector<std::string> fParNamesglobal = resultfinal->fittedParameterNames();
355 std::vector<double> fParsglobal = resultfinal->fittedParameters();
357 std::cout<<
" Global Fit "<<std::endl;
358 for(
unsigned int i=0; i<fParsglobal.size();i++)
360 std::cout<<
" Parameter "<< fParNamesglobal[i]<<
" "<<fParsglobal[i]<<std::endl;
362 std::cout<<
"Chi^2/ndf: "<<resultfinal->quality()<<std::endl;
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());
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]));
383 zinitialPars.push_back(100);
384 zinitialPars.push_back(0);
386 zinitialParsgauss.push_back(5000);
387 zinitialParsgauss.push_back(0);
388 zinitialParsgauss.push_back(2);
390 functiongauss->setParameters( zinitialParsgauss );
392 functionexp1->setParameters( zinitialPars );
393 functionexp2->setParameters( zinitialPars );
396 fitter->fitParameterSettings(
"mean").setBounds(-0.0000000001,0.0000000001);
398 AIDA::IFitResult* zresultgauss = fitter->fit(*zhistogauss, *functiongauss);
399 AIDA::IFitResult* zresultexp1 = fitter->fit(*zhistoexp1, *functionexp1);
400 AIDA::IFitResult* zresultexp2 = fitter->fit(*zhistoexp2, *functionexp2);
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();
410 std::cout<<std::endl;
411 std::cout<<std::endl;
412 std::cout<<std::endl;
414 std::cout<<
" Z - axis "<<std::endl;
415 std::cout<<std::endl;
417 std::cout<<
" Local Fit of Gaussian "<<std::endl;
418 for(
unsigned int i=0; i<fParsgauss.size();i++)
420 std::cout<<
" Parameter "<< zfParNamesgauss[i]<<
" "<<zfParsgauss[i]<<std::endl;
421 zglobalPars.push_back(zfParsgauss[i]);
423 std::cout<<
"Chi^2/ndf: "<<zresultgauss->quality()<<std::endl;
425 std::cout<<
" Local Fit of first Exponential "<<std::endl;
426 for(
unsigned int i=0; i<zfParsexp1.size();i++)
428 std::cout<<
" Parameter "<< zfParNamesexp1[i]<<
" "<<zfParsexp1[i]<<std::endl;
429 zglobalPars.push_back(zfParsexp1[i]);
431 std::cout<<
"Chi^2/ndf: "<<zresultexp1->quality()<<std::endl;
433 std::cout<<
" Local Fit of second Exponential "<<std::endl;
434 for(
unsigned int i=0; i<zfParsexp2.size();i++)
436 std::cout<<
" Parameter "<< zfParNamesexp1[i]<<
" "<<zfParsexp2[i]<<std::endl;
437 globalPars.push_back(zfParsexp2[i]);
440 std::cout<<
"Chi^2/ndf: "<<zresultexp2->quality()<<std::endl;
443 functionfinal->setParameters(zglobalPars);
444 AIDA::IFitResult* zresultfinal = fitter->fit(*zhisto, *functionfinal);
446 std::vector<std::string> zfParNamesglobal = zresultfinal->fittedParameterNames();
447 std::vector<double> zfParsglobal = zresultfinal->fittedParameters();
449 std::cout<<
" Global Fit "<<std::endl;
450 for(
unsigned int i=0; i<zfParsglobal.size();i++)
452 std::cout<<
" Parameter "<< zfParNamesglobal[i]<<
" "<<zfParsglobal[i]<<std::endl;
454 std::cout<<
"Chi^2/ndf: "<<zresultfinal->quality()<<std::endl;
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());
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]);
475 std::cout<<
"PARAMETERS FOR RPHI Joint Probability"<<std::endl;
476 for(
unsigned int i=0;i<5;i++)
478 std::cout<<outPars[i]<<std::endl;
480 std::cout<<
"PARAMETERS FOR Z Joint Probability"<<std::endl;
481 for(
unsigned int i=0;i<5;i++)
483 std::cout<<zoutPars[i]<<std::endl;
486 MetaMemoryManager::Run()->delAllObjects();
489 std::cout <<
"SignificanceFitProcessor::end() " << name()
490 <<
" processed " << _nEvt <<
" events in " << _nRun <<
" runs "
497 #endif // end if ifdef MARLIN_USE_AIDA
void registerObject(T *pointer)
Register an object for memory management.