WCSimVertexFinder.cc

Go to the documentation of this file.
00001 #include "WCSimVertexFinder.hh"
00002 #include "WCSimVertexGeometry.hh"
00003 
00004 #include "WCSimInterface.hh"
00005 #include "WCSimGeometry.hh"
00006 #include "WCSimParameters.hh"
00007 
00008 #include "WCSimRecoDigit.hh"
00009 #include "WCSimRecoEvent.hh"
00010 #include "WCSimTrueEvent.hh"
00011 
00012 #include "TMath.h"
00013 
00014 #include <cmath>
00015 #include <iostream>
00016 #include <cassert>
00017 
00018 ClassImp(WCSimVertexFinder)
00019 
00020 static WCSimVertexFinder* fgVertexFinder = 0;
00021 
00022 static void point_position_chi2(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00023 {
00024   Bool_t printDebugMessages = 0;
00025 
00026   Double_t vtxX = par[0]; // centimetres
00027   Double_t vtxY = par[1]; // centimetres
00028   Double_t vtxZ = par[2]; // centimetres
00029 
00030   Double_t vtime  = 0.0;
00031   Double_t fom = 0.0;
00032 
00033   fgVertexFinder->point_position_itr();
00034   fgVertexFinder->PointPositionChi2(vtxX,vtxY,vtxZ,
00035                                     vtime,fom);
00036 
00037   f = -fom; // note: need to maximize this fom
00038 
00039   if( printDebugMessages ){
00040     std::cout << " [point_position_chi2] [" << fgVertexFinder->point_position_iterations() << "] (x,y,z)=(" << vtxX << "," << vtxY << "," << vtxZ << ") vtime=" << vtime << " fom=" << fom << std::endl;
00041   }
00042 
00043   return;
00044 }
00045 
00046 static void point_direction_chi2(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00047 {
00048   Bool_t printDebugMessages = 0;
00049   
00050   Double_t vtxX = fgVertexFinder->fVtxX;
00051   Double_t vtxY = fgVertexFinder->fVtxY;
00052   Double_t vtxZ = fgVertexFinder->fVtxZ;
00053 
00054   Double_t dirTheta = par[0]; // radians
00055   Double_t dirPhi   = par[1]; // radians
00056   
00057   Double_t dirX = sin(dirTheta)*cos(dirPhi);
00058   Double_t dirY = sin(dirTheta)*sin(dirPhi);
00059   Double_t dirZ = cos(dirTheta);
00060 
00061   Double_t vangle = 0.0;
00062   Double_t fom = 0.0;
00063 
00064   fgVertexFinder->point_direction_itr();
00065   fgVertexFinder->PointDirectionChi2(vtxX,vtxY,vtxZ,
00066                                      dirX,dirY,dirZ,
00067                                      vangle,fom);
00068 
00069   f = -fom; // note: need to maximize this fom
00070 
00071   if( printDebugMessages ){
00072     std::cout << " [point_direction_chi2] [" << fgVertexFinder->point_direction_iterations() << "] (px,py,pz)=(" << dirX << "," << dirY << "," << dirZ << ") fom=" << fom << std::endl;
00073   }
00074 
00075   return; 
00076 }
00077 
00078 static void point_vertex_chi2(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00079 {
00080   Bool_t printDebugMessages = 0;
00081   
00082   Double_t vtxX     = par[0]; // centimetres
00083   Double_t vtxY     = par[1]; // centimetres
00084   Double_t vtxZ     = par[2]; // centimetres  
00085 
00086   Double_t dirTheta = par[3]; // radians
00087   Double_t dirPhi   = par[4]; // radians
00088 
00089   Double_t dirX = sin(dirTheta)*cos(dirPhi);
00090   Double_t dirY = sin(dirTheta)*sin(dirPhi);
00091   Double_t dirZ = cos(dirTheta);
00092 
00093   Double_t vangle = 0.0;
00094   Double_t vtime  = 0.0;
00095   Double_t fom = 0.0;
00096 
00097   fgVertexFinder->point_vertex_itr();
00098   fgVertexFinder->PointVertexChi2(vtxX,vtxY,vtxZ,
00099                                   dirX,dirY,dirZ,
00100                                   vangle,vtime,fom);
00101 
00102   f = -fom; // note: need to maximize this fom
00103 
00104   if( printDebugMessages ){
00105     std::cout << " [point_vertex_chi2] [" << fgVertexFinder->point_vertex_iterations() << "] (x,y,z)=(" << vtxX << "," << vtxY << "," << vtxZ << ") (px,py,pz)=(" << dirX << "," << dirY << "," << dirZ << ") vtime=" << vtime << " fom=" << fom << std::endl;
00106   }
00107 
00108   return;
00109 }
00110 
00111 static void extended_vertex_chi2(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00112 {
00113   Bool_t printDebugMessages = 0;
00114   
00115   Double_t vtxX     = par[0]; // centimetres
00116   Double_t vtxY     = par[1]; // centimetres
00117   Double_t vtxZ     = par[2]; // centimetres  
00118 
00119   Double_t dirTheta = par[3]; // radians
00120   Double_t dirPhi   = par[4]; // radians
00121 
00122   Double_t dirX = sin(dirTheta)*cos(dirPhi);
00123   Double_t dirY = sin(dirTheta)*sin(dirPhi);
00124   Double_t dirZ = cos(dirTheta);
00125 
00126   Double_t vangle = 0.0;
00127   Double_t vtime  = 0.0;
00128   Double_t fom = 0.0;
00129 
00130   fgVertexFinder->extended_vertex_itr();
00131   fgVertexFinder->ExtendedVertexChi2(vtxX,vtxY,vtxZ,
00132                                      dirX,dirY,dirZ,
00133                                      vangle,vtime,fom);
00134 
00135   f = -fom; // note: need to maximize this fom
00136 
00137   if( printDebugMessages ){
00138     std::cout << " [extended_vertex_chi2] [" << fgVertexFinder->extended_vertex_iterations() << "] (x,y,z)=(" << vtxX << "," << vtxY << "," << vtxZ << ") (px,py,pz)=(" << dirX << "," << dirY << "," << dirZ << ") vtime=" << vtime << " fom=" << fom << std::endl;
00139   }
00140 
00141   return;
00142 }
00143 
00144 static void vertex_time_fom(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00145 {   
00146   Double_t vtxTime = par[0]; // nanoseconds
00147 
00148   Double_t fom = 0.0;
00149 
00150   fgVertexFinder->TimePropertiesFoM(vtxTime,fom);
00151 
00152   f = -fom; // note: need to maximize this fom
00153 
00154   return;
00155 }
00156 
00157 static void vertex_time_lnl(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00158 {  
00159   Bool_t printDebugMessages = 0;
00160   
00161   Double_t vtxTime = par[0]; // nanoseconds
00162   Double_t vtxParam0 = fgVertexFinder->fFixTimeParam0;
00163 
00164   if( fgVertexFinder->fFitTimeParams ){
00165     vtxParam0 = par[1]; 
00166   }
00167 
00168   Double_t fom = 0.0;
00169 
00170   fgVertexFinder->time_fit_itr();
00171   fgVertexFinder->TimePropertiesLnL(vtxTime,vtxParam0,fom);
00172 
00173   f = -fom; // note: need to maximize this fom
00174 
00175   if( printDebugMessages ){
00176     std::cout << "  [vertex_time_lnl] [" << fgVertexFinder->time_fit_iterations() << "] vtime=" << vtxTime << " vparam=" << vtxParam0 << " fom=" << fom << std::endl;
00177   }
00178 
00179   return;
00180 }
00181 
00182 static void vertex_cone_lnl(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00183 {
00184   Bool_t printDebugMessages = 0;
00185   
00186   Double_t vtxParam0 = fgVertexFinder->fFixConeParam0;
00187   Double_t vtxParam1 = fgVertexFinder->fFixConeParam1;
00188   Double_t vtxParam2 = fgVertexFinder->fFixConeParam2;
00189 
00190   if( fgVertexFinder->fFitConeParams ){
00191     vtxParam0 = par[0];
00192     vtxParam1 = par[1];
00193     vtxParam2 = par[2];
00194   }
00195  
00196   Double_t vangle = 0.0;
00197   Double_t fom = 0.0;
00198 
00199   fgVertexFinder->cone_fit_itr();
00200   fgVertexFinder->ConePropertiesLnL(vtxParam0,vtxParam1,vtxParam2,vangle,fom);
00201 
00202   f = -fom; // note: need to maximize this fom
00203 
00204   if( printDebugMessages ){
00205     std::cout << "  [vertex_cone_lnl] [" << fgVertexFinder->cone_fit_iterations() << "] vparam0=" << vtxParam0 << " vparam1=" << vtxParam1 << " vtxParam2=" << vtxParam2 << " fom=" << fom << std::endl;
00206   }
00207 
00208   return;
00209 }
00210 
00211 WCSimVertexFinder* WCSimVertexFinder::Instance()
00212 {
00213   if( !fgVertexFinder ){
00214     fgVertexFinder = new WCSimVertexFinder();
00215   }
00216 
00217   if( !fgVertexFinder ){
00218     assert(fgVertexFinder);
00219   }
00220 
00221   if( fgVertexFinder ){
00222 
00223   }
00224 
00225   return fgVertexFinder;
00226 }
00227 
00228 void WCSimVertexFinder::PointFitOnly(Bool_t yesno)
00229 {
00230   WCSimVertexFinder::Instance()->SetPointFitOnly(yesno);
00231 }
00232 
00233 void WCSimVertexFinder::UseTrueVertex(Bool_t yesno)
00234 {
00235   WCSimVertexFinder::Instance()->UsingTrueVertex(yesno);
00236 }
00237 
00238 void WCSimVertexFinder::UseTruePosition(Bool_t yesno)
00239 {
00240   WCSimVertexFinder::Instance()->UsingTruePosition(yesno);
00241 }
00242 
00243 void WCSimVertexFinder::UseTrueDirection(Bool_t yesno)
00244 {
00245   WCSimVertexFinder::Instance()->UsingTrueDirection(yesno);
00246 }
00247 
00248 void WCSimVertexFinder::SeedWithTrueVertex(Bool_t yesno)
00249 {
00250   WCSimVertexFinder::Instance()->SeedingWithTrueVertex(yesno);
00251 }
00252 
00253 void WCSimVertexFinder::SeedWithSimpleVertex(Bool_t yesno)
00254 {
00255   WCSimVertexFinder::Instance()->SeedingWithSimpleVertex(yesno);
00256 }
00257 
00258 void WCSimVertexFinder::SeedWithQuadruples(Bool_t yesno)
00259 {
00260   WCSimVertexFinder::Instance()->SeedingWithQuadruples(yesno);
00261 }
00262 
00263 void WCSimVertexFinder::NumSeeds(Int_t nseeds)
00264 {
00265   WCSimVertexFinder::Instance()->SetNumSeeds(nseeds);
00266 }
00267 
00268 void WCSimVertexFinder::FitWeights(Double_t tw, Double_t cw)
00269 {
00270   WCSimVertexFinder::Instance()->SetFitWeights(tw,cw);
00271 }
00272 
00273 void WCSimVertexFinder::FixTimeParams(Double_t param0)
00274 {
00275   WCSimVertexFinder::Instance()->SetFixTimeParams(param0);
00276 }
00277 
00278 void WCSimVertexFinder::FitTimeParams()
00279 {
00280   WCSimVertexFinder::Instance()->SetFitTimeParams();
00281 }  
00282 
00283 void WCSimVertexFinder::FixConeParams(Double_t param0, Double_t param1, Double_t param2)
00284 {
00285   WCSimVertexFinder::Instance()->SetFixConeParams(param0,param1,param2);
00286 }
00287 
00288 void WCSimVertexFinder::FitConeParams()
00289 {
00290   WCSimVertexFinder::Instance()->SetFitConeParams();
00291 }  
00292 
00293 void WCSimVertexFinder::FixVertexBias(Double_t bias)
00294 {
00295   WCSimVertexFinder::Instance()->SetVertexBias(bias);
00296 }
00297 
00298 void WCSimVertexFinder::PrintParameters()
00299 {
00300   WCSimVertexFinder::Instance()->RunPrintParameters();
00301 }
00302 
00303 WCSimVertexFinder::WCSimVertexFinder()
00304 {
00305   // default configuration
00306   fBaseFOM = 100.0;
00307   fPointFitOnly = 0;
00308   fUseTrueVertex = 0; 
00309   fUseTruePosition = 0;
00310   fUseTrueDirection = 0;
00311   fSeedWithTrueVertex = 0;
00312   fSeedWithSimpleVertex = 0;
00313   fSeedWithQuadruples = 1;
00314   fNumSeeds = 200;
00315 
00316   fFitTimeParams = 0;     // don't fit by default
00317   fFixTimeParam0 = 0.20;  // scattering parameter
00318 
00319   fFitConeParams = 1;     // do fit by default  
00320   fFixConeParam0 = 0.25;  // track length parameter
00321   fFixConeParam1 = 0.50;  // track length parameter
00322   fFixConeParam2 = 0.75;  // particle ID:  0.0[electron]->1.0[muon]
00323 
00324   fTimeFitWeight = 0.50;  // nominal time weight
00325   fConeFitWeight = 0.50;  // nominal cone weight
00326 
00327   fFixVertexBias = 1;     // fix vertex bias
00328   fVertexBias = 25.0;     // size of vertex bias [cm]
00329 
00330   fIntegralsDone = 0;
00331 
00332   fVtxX = 0.0;
00333   fVtxY = 0.0;
00334   fVtxZ = 0.0;
00335   fDirX = 0.0;
00336   fDirY = 0.0;
00337   fDirZ = 0.0;
00338   fVtxFOM = 0.0;
00339 
00340   fSimplePosition = 0;
00341   fSimpleDirection = 0;
00342   fPointPosition = 0;
00343   fPointDirection = 0;
00344   fPointVertex = 0;
00345   fExtendedVertex = 0;
00346 
00347   fMinuitPointPosition = new TMinuit();
00348   fMinuitPointPosition->SetPrintLevel(-1);
00349   fMinuitPointPosition->SetMaxIterations(5000);
00350 
00351   fMinuitPointDirection = new TMinuit();
00352   fMinuitPointDirection->SetPrintLevel(-1);
00353   fMinuitPointDirection->SetMaxIterations(5000);
00354 
00355   fMinuitPointVertex = new TMinuit();
00356   fMinuitPointVertex->SetPrintLevel(-1);
00357   fMinuitPointVertex->SetMaxIterations(5000);
00358 
00359   fMinuitExtendedVertex = new TMinuit();
00360   fMinuitExtendedVertex->SetPrintLevel(-1);
00361   fMinuitExtendedVertex->SetMaxIterations(5000); 
00362 
00363   fMinuitTimeFit = new TMinuit();
00364   fMinuitTimeFit->SetPrintLevel(-1);
00365   fMinuitTimeFit->SetMaxIterations(5000);   
00366 
00367   fMinuitConeFit = new TMinuit();
00368   fMinuitConeFit->SetPrintLevel(-1);
00369   fMinuitConeFit->SetMaxIterations(5000);   
00370 
00371   fPass = 0;
00372   fItr = 0; 
00373 
00374   fTimeFitItr = 0;
00375   fConeFitItr = 0;
00376   fPointPosItr = 0;
00377   fPointDirItr = 0;
00378   fPointVtxItr = 0;
00379   fExtendedVtxItr = 0;  
00380 
00381   // --- debug ---
00382   fTimeParam0 = 0.0;
00383   fConeParam0 = 0.0;
00384   fConeParam1 = 0.0; 
00385   fConeParam2 = 0.0;
00386 }
00387 
00388 WCSimVertexFinder::~WCSimVertexFinder()
00389 {
00390   // delete fitter
00391   // =============
00392   delete fMinuitPointPosition;
00393   delete fMinuitPointDirection;
00394   delete fMinuitPointVertex;
00395   delete fMinuitExtendedVertex;
00396 
00397   delete fMinuitTimeFit;
00398   delete fMinuitConeFit;
00399 
00400   // clear vertices
00401   // ==============
00402   for(UInt_t i=0; i<vVertexList.size(); i++ ){
00403     delete (WCSimRecoVertex*)(vVertexList.at(i));
00404   }
00405 
00406   vVertexList.clear();
00407 }
00408 
00409 void WCSimVertexFinder::RunPrintParameters()
00410 {
00411   std::cout << " *** WCSimVertexFinder::PrintParameters() *** " << std::endl;
00412 
00413   std::cout << "  Vertex Finding Parameters: " << std::endl
00414             << "   PointFitOnly = " << fPointFitOnly << std::endl
00415             << "   UseTrueVertex = " << fUseTrueVertex << std::endl
00416             << "   UseTruePosition = " << fUseTruePosition << std::endl
00417             << "   UseTrueDirection = " << fUseTrueDirection << std::endl
00418             << "   SeedWithTrueVertex = " << fSeedWithTrueVertex << std::endl
00419             << "   SeedWithSimpleVertex = " << fSeedWithSimpleVertex << std::endl
00420             << "   SeedWithQuadruples = " << fSeedWithQuadruples << std::endl
00421             << "   NumSeeds = " << fNumSeeds << std::endl
00422             << "  Vertex Fitting Parameters: " << std::endl
00423             << "   BaseFOM = " << fBaseFOM << std::endl
00424             << "   FitTimeParams = " << fFitTimeParams << " (" << fFixTimeParam0 << ") " << std::endl
00425             << "   FitConeParams = " << fFitConeParams << " (" << fFixConeParam0 << "," << fFixConeParam1 << "," << fFixConeParam2 << ") " << std::endl
00426             << "   Weights = (" << fTimeFitWeight << "," << fConeFitWeight << ") " << std::endl
00427             << "   FixVertexBias = " << fFixVertexBias << " (" << fVertexBias << ") " << std::endl;
00428 
00429   return;
00430 }
00431 
00432 void WCSimVertexFinder::Reset()
00433 {
00434   return this->Clear();
00435 }
00436 
00437 void WCSimVertexFinder::Clear()
00438 {
00439   // clear vertices
00440   // ==============
00441   for(UInt_t i=0; i<vVertexList.size(); i++ ){
00442     delete (WCSimRecoVertex*)(vVertexList.at(i));
00443   }
00444 
00445   vVertexList.clear();
00446 
00447   // clear vertices
00448   // ==============
00449   fSimplePosition = 0;
00450   fSimpleDirection = 0;
00451   fPointPosition = 0;
00452   fPointDirection = 0; 
00453   fPointVertex = 0;
00454   fExtendedVertex = 0;
00455 
00456   return;
00457 }
00458 
00459 
00460 WCSimRecoVertex* WCSimVertexFinder::Run(WCSimRecoEvent* myEvent)
00461 {
00462   this->Clear();
00463 
00464   if( fPointFitOnly ){
00465     return (WCSimRecoVertex*)(this->RunPointFit(myEvent));
00466   }
00467   else{
00468     return (WCSimRecoVertex*)(this->RunExtendedFit(myEvent));
00469   }
00470 }
00471 
00472 WCSimRecoVertex* WCSimVertexFinder::RunPointFit(WCSimRecoEvent* myEvent)
00473 {
00474   // return true vertex
00475   // ==================
00476   if( fUseTrueVertex ){
00477     return this->RunPointFitFromTruth(myEvent);
00478   }
00479 
00480   // load event
00481   // ==========
00482   WCSimVertexGeometry::Instance()->LoadEvent(myEvent);
00483 
00484   std::cout << " *** WCSimVertexFinder::RunPointFit(...) *** " << std::endl;
00485 
00486   // point positiobn fit
00487   // ===================
00488   WCSimRecoVertex* simplePos = (WCSimRecoVertex*)(this->FindSimplePosition());
00489   WCSimRecoVertex* pointPos  = (WCSimRecoVertex*)(this->FitPointPosition(simplePos));
00490 
00491   // point direction fit
00492   // ===================
00493   WCSimRecoVertex* simpleDir = (WCSimRecoVertex*)(this->FindSimpleDirection(pointPos));
00494   WCSimRecoVertex* pointDir  = (WCSimRecoVertex*)(this->FitPointDirection(simpleDir));  
00495 
00496   // point vertex fit
00497   // ================
00498   WCSimRecoVertex* pointVtx  = (WCSimRecoVertex*)(this->FitPointVertex(pointDir)); 
00499 
00500   // return vertex
00501   // =============
00502   return pointVtx;
00503 }
00504 
00505 WCSimRecoVertex* WCSimVertexFinder::RunExtendedFit(WCSimRecoEvent* myEvent)
00506 {    
00507   // return true vertex
00508   // ==================
00509   if( fUseTrueVertex ){  
00510     return this->RunExtendedFitFromTruth(myEvent);
00511   }
00512 
00513   // point fit
00514   // =========
00515   WCSimRecoVertex* pointVtx = (WCSimRecoVertex*)(this->RunPointFit(myEvent));
00516 
00517   std::cout << " *** WCSimVertexFinder::RunExtendedFit(...) *** " << std::endl;
00518 
00519   // extended fit to vertex and direction
00520   // ====================================
00521   WCSimRecoVertex* extendedVtx = (WCSimRecoVertex*)(this->FitExtendedVertex(pointVtx));
00522 
00523   // correct bias in vertex and direction
00524   // ====================================
00525   WCSimRecoVertex* extendedVtxCorrected = (WCSimRecoVertex*)(this->CorrectExtendedVertex(extendedVtx));
00526 
00527   // return vertex
00528   // =============
00529   return extendedVtxCorrected;
00530 }
00531 
00532 WCSimRecoVertex* WCSimVertexFinder::RunPointFitFromTruth(WCSimRecoEvent* myEvent)
00533 {
00534   std::cout << " *** WCSimVertexFinder::RunPointFit(...) *** " << std::endl;
00535   std::cout << "  --- reconstruct vertex from truth --- " << std::endl;
00536   
00537   WCSimTrueEvent* myTrueEvent = (WCSimTrueEvent*)(WCSimInterface::TrueEvent());
00538 
00539   if( myTrueEvent->GetNTracks()>0 ){
00540     Double_t vtxX = myTrueEvent->GetVtxX();
00541     Double_t vtxY = myTrueEvent->GetVtxY();
00542     Double_t vtxZ = myTrueEvent->GetVtxZ();
00543     Double_t dirX = myTrueEvent->GetDirX();
00544     Double_t dirY = myTrueEvent->GetDirY();
00545     Double_t dirZ = myTrueEvent->GetDirZ();
00546 
00547     return (WCSimRecoVertex*)(this->RunPointFit(myEvent, 
00548                                                  vtxX,vtxY,vtxZ, 
00549                                                   dirX,dirY,dirZ));
00550   }
00551   else{
00552     return (WCSimRecoVertex*)(this->BuildDummyVertex());
00553   }
00554 }
00555 
00556 WCSimRecoVertex* WCSimVertexFinder::RunExtendedFitFromTruth(WCSimRecoEvent* myEvent)
00557 {
00558   std::cout << " *** WCSimVertexFinder::RunExtendedFit(...) *** " << std::endl;
00559   std::cout << "  --- reconstruct vertex from truth --- " << std::endl;
00560   
00561   WCSimTrueEvent* myTrueEvent = (WCSimTrueEvent*)(WCSimInterface::TrueEvent());
00562 
00563   if( myTrueEvent->GetNTracks()>0 ){
00564     Double_t vtxX = myTrueEvent->GetVtxX();
00565     Double_t vtxY = myTrueEvent->GetVtxY();
00566     Double_t vtxZ = myTrueEvent->GetVtxZ();
00567     Double_t dirX = myTrueEvent->GetDirX();
00568     Double_t dirY = myTrueEvent->GetDirY();
00569     Double_t dirZ = myTrueEvent->GetDirZ();
00570 
00571     return (WCSimRecoVertex*)(this->RunExtendedFit(myEvent, 
00572                                                     vtxX,vtxY,vtxZ, 
00573                                                      dirX,dirY,dirZ));
00574   }
00575   else{
00576     return (WCSimRecoVertex*)(this->BuildDummyVertex());
00577   }
00578 }
00579 
00580 WCSimRecoVertex* WCSimVertexFinder::Run(WCSimRecoEvent* myEvent, Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00581 {
00582   this->Clear();
00583 
00584   if( fPointFitOnly ){
00585     return (WCSimRecoVertex*)(this->RunPointFit(myEvent, 
00586                                                  vtxX,vtxY,vtxZ, 
00587                                                   dirX,dirY,dirZ));
00588   }
00589   else{
00590     return (WCSimRecoVertex*)(this->RunExtendedFit(myEvent, 
00591                                                     vtxX,vtxY,vtxZ, 
00592                                                      dirX,dirY,dirZ));
00593   }
00594 }
00595 
00596 WCSimRecoVertex* WCSimVertexFinder::RunPointFit(WCSimRecoEvent* myEvent, Double_t vtxX, Double_t vtxY, Double_t vtxZ)
00597 {   
00598   // load event
00599   // ==========
00600   WCSimVertexGeometry::Instance()->LoadEvent(myEvent);
00601  
00602   // fix point vertex
00603   // ================
00604   WCSimRecoVertex* simpleVtx = (WCSimRecoVertex*)(this->FixSimplePosition(vtxX,vtxY,vtxZ));
00605   WCSimRecoVertex* pointPos = (WCSimRecoVertex*)(this->FixPointPosition(vtxX,vtxY,vtxZ));
00606 
00607   // return vertex
00608   // =============
00609   return pointPos;
00610 }
00611 
00612 WCSimRecoVertex* WCSimVertexFinder::RunPointFit(WCSimRecoEvent* myEvent, Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00613 {   
00614   // fix point vertex
00615   // ================
00616   WCSimRecoVertex* pointPos = (WCSimRecoVertex*)(this->RunPointFit(myEvent,
00617                                                                     vtxX,vtxY,vtxZ));
00618 
00619   // fix point direction
00620   // ===================
00621   WCSimRecoVertex* simpleDir = (WCSimRecoVertex*)(this->FixSimpleDirection(vtxX,vtxY,vtxZ,
00622                                                                             dirX,dirY,dirZ));
00623   WCSimRecoVertex* pointDir = (WCSimRecoVertex*)(this->FixPointDirection(vtxX,vtxY,vtxZ,
00624                                                                           dirX,dirY,dirZ));
00625 
00626   // fix point vertex
00627   // ================
00628   WCSimRecoVertex* pointVtx = (WCSimRecoVertex*)(this->FixPointVertex(vtxX,vtxY,vtxZ,
00629                                                                        dirX,dirY,dirZ));
00630 
00631   // return vertex
00632   // =============
00633   return pointVtx;
00634 }
00635 
00636 WCSimRecoVertex* WCSimVertexFinder::RunExtendedFit(WCSimRecoEvent* myEvent, Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00637 {     
00638   // fix point direction
00639   // ===================
00640   WCSimRecoVertex* pointVtx = (WCSimRecoVertex*)(this->RunPointFit(myEvent,
00641                                                                     vtxX,vtxY,vtxZ,
00642                                                                      dirX,dirY,dirZ));
00643 
00644   // fix extended vertex
00645   // ===================
00646   WCSimRecoVertex* extendedVtx = (WCSimRecoVertex*)(this->FixExtendedVertex(vtxX,vtxY,vtxZ,
00647                                                                              dirX,dirY,dirZ));
00648 
00649   // return vertex
00650   // =============
00651   return extendedVtx;
00652 }
00653 
00654 WCSimRecoVertex* WCSimVertexFinder::BuildTrueVertex()
00655 {
00656   WCSimTrueEvent* myTrueEvent = (WCSimTrueEvent*)(WCSimInterface::TrueEvent());
00657 
00658   Double_t vtxX = myTrueEvent->GetVtxX();
00659   Double_t vtxY = myTrueEvent->GetVtxY();
00660   Double_t vtxZ = myTrueEvent->GetVtxZ();  
00661 
00662   return (WCSimRecoVertex*)(this->FixSimplePosition(vtxX,vtxY,vtxZ));
00663 }
00664 
00665 WCSimRecoVertex* WCSimVertexFinder::BuildDummyVertex()
00666 {
00667   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00668   vVertexList.push_back(newVertex);
00669 
00670   fSimplePosition  = newVertex;
00671   fSimpleDirection = newVertex;
00672   fPointPosition   = newVertex;
00673   fPointDirection  = newVertex;
00674   fPointVertex     = newVertex;
00675   fExtendedVertex  = newVertex;
00676 
00677   std::cout << "  --- no true vertex, building dummy vertex --- " << std::endl;
00678 
00679   return newVertex;
00680 }
00681 
00682 WCSimRecoVertex* WCSimVertexFinder::FindSimplePosition()
00683 {
00684   if( fSeedWithTrueVertex ){
00685     std::cout << "  --- seed vertex with truth --- " << std::endl;
00686     return (WCSimRecoVertex*)(this->BuildTrueVertex());
00687   }
00688 
00689   if( fSeedWithSimpleVertex ){
00690     fSimplePosition = (WCSimRecoVertex*)(WCSimVertexGeometry::Instance()->CalcSimpleVertex()); 
00691   }
00692   else if( fSeedWithQuadruples ){
00693     fSimplePosition = (WCSimRecoVertex*)(this->FindSeedPosition());
00694   }
00695   else{
00696     fSimplePosition = (WCSimRecoVertex*)(WCSimVertexGeometry::Instance()->CalcSimpleVertex());  
00697   }
00698 
00699   std::cout << "  simple vertex: " << std::endl
00700             << "    (vx,vy,vz)=(" << fSimplePosition->GetX() << "," << fSimplePosition->GetY() << "," << fSimplePosition->GetZ() << ") " << std::endl
00701             << "      vtime=" << fSimplePosition->GetTime() << " itr=" << fSimplePosition->GetIterations() << " fom=" << fSimplePosition->GetFOM() << std::endl;
00702 
00703   if( fSimplePosition->GetPass()==0 ) std::cout << "   <warning> simple vertex calculation failed! " << std::endl;
00704 
00705   return fSimplePosition;
00706 }
00707 
00708 WCSimRecoVertex* WCSimVertexFinder::FindSimpleDirection(WCSimRecoVertex* myVertex)
00709 {
00710   fSimpleDirection = (WCSimRecoVertex*)(this->FindSeedDirection(myVertex));
00711 
00712   std::cout << "  simple direction: " << std::endl
00713             << "    (vx,vy,vz)=(" << fSimpleDirection->GetX() << "," << fSimpleDirection->GetY() << "," << fSimpleDirection->GetZ() << ") " << std::endl
00714             << "    (px,py,pz)=(" << fSimpleDirection->GetDirX() << "," << fSimpleDirection->GetDirY() << "," << fSimpleDirection->GetDirZ() << ") " << std::endl
00715             << "      vtime=" << fSimpleDirection->GetTime() << " itr=" << fSimpleDirection->GetIterations() << " fom=" << fSimpleDirection->GetFOM() << std::endl;
00716 
00717   if( fSimpleDirection->GetPass()==0 ) std::cout << "   <warning> simple direction calculation failed! " << std::endl;
00718 
00719   return fSimpleDirection;
00720 }
00721 
00722 WCSimRecoVertex* WCSimVertexFinder::FindSeedPosition()
00723 {
00724   Double_t vtxX = 0.0;
00725   Double_t vtxY = 0.0;
00726   Double_t vtxZ = 0.0;
00727   Double_t vtxTime = 950.0;
00728   Double_t vtxFOM = 0.0;
00729 
00730   Int_t bestSeed = -1;
00731   Double_t bestFOM = -1.0;
00732 
00733   WCSimVertexGeometry::Instance()->CalcVertexSeeds(fNumSeeds);
00734   UInt_t nlast = WCSimVertexGeometry::Instance()->GetNSeeds();
00735 
00736   for( UInt_t n=0; n<nlast; n++ ){
00737     vtxX = WCSimVertexGeometry::Instance()->GetSeedVtxX(n);
00738     vtxY = WCSimVertexGeometry::Instance()->GetSeedVtxY(n);
00739     vtxZ = WCSimVertexGeometry::Instance()->GetSeedVtxZ(n);
00740 
00741     this->PointPositionChi2(vtxX,vtxY,vtxZ,
00742                             vtxTime,vtxFOM);
00743 
00744     if( vtxFOM>bestFOM ){
00745       bestSeed = n; 
00746       bestFOM = vtxFOM;
00747     }
00748   }
00749   
00750   if( bestSeed>=0 ){
00751     vtxX = WCSimVertexGeometry::Instance()->GetSeedVtxX(bestSeed);
00752     vtxY = WCSimVertexGeometry::Instance()->GetSeedVtxY(bestSeed);
00753     vtxZ = WCSimVertexGeometry::Instance()->GetSeedVtxZ(bestSeed);
00754   }
00755 
00756   WCSimRecoVertex* newVertex = (WCSimRecoVertex*)(this->FixSimplePosition(vtxX,vtxY,vtxZ));
00757   newVertex->SetFOM(fBaseFOM,nlast,1); 
00758 
00759   return newVertex;
00760 }
00761 
00762 WCSimRecoVertex* WCSimVertexFinder::FindSeedDirection(WCSimRecoVertex* myVertex)
00763 {
00764   WCSimRecoVertex* seedVertex = (WCSimRecoVertex*)(WCSimVertexGeometry::Instance()->CalcSimpleDirection(myVertex));  
00765 
00766   Double_t vtxX = seedVertex->GetX();
00767   Double_t vtxY = seedVertex->GetY();
00768   Double_t vtxZ = seedVertex->GetZ();
00769 
00770   Double_t dirX = seedVertex->GetDirX();
00771   Double_t dirY = seedVertex->GetDirY();
00772   Double_t dirZ = seedVertex->GetDirZ();
00773 
00774   WCSimRecoVertex* newVertex = (WCSimRecoVertex*)(this->FixSimpleDirection(vtxX,vtxY,vtxZ,
00775                                                                            dirX,dirY,dirZ));
00776   return newVertex;
00777 }
00778 
00779 WCSimRecoVertex* WCSimVertexFinder::FitPointPosition(WCSimRecoVertex* myVertex)
00780 {  
00781   return (WCSimRecoVertex*)(this->FitPointPositionWithMinuit(myVertex));
00782 }
00783 
00784 WCSimRecoVertex* WCSimVertexFinder::FitPointDirection(WCSimRecoVertex* myVertex)
00785 {  
00786   return (WCSimRecoVertex*)(this->FitPointDirectionWithMinuit(myVertex));
00787 }
00788 
00789 WCSimRecoVertex* WCSimVertexFinder::FitPointVertex(WCSimRecoVertex* myVertex)
00790 {      
00791   return (WCSimRecoVertex*)(this->FitPointVertexWithMinuit(myVertex));
00792 }
00793 
00794 WCSimRecoVertex* WCSimVertexFinder::FitExtendedVertex(WCSimRecoVertex* myVertex)
00795 {
00796   return (WCSimRecoVertex*)(this->FitExtendedVertexWithMinuit(myVertex));
00797 }
00798 
00799 WCSimRecoVertex* WCSimVertexFinder::FixSimplePosition(Double_t vtxX, Double_t vtxY, Double_t vtxZ)
00800 {  
00801   // initialization
00802   // ==============
00803   Double_t vtxTime = 950.0;
00804   Double_t vtxFOM = fBaseFOM;
00805 
00806   // create new vertex
00807   // =================
00808   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00809   vVertexList.push_back(newVertex);
00810   fSimplePosition = newVertex;
00811 
00812   // set vertex
00813   // ==========
00814   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00815   newVertex->SetFOM(vtxFOM,1,1); 
00816  
00817   // print vertex
00818   // ============
00819   std::cout << "  set simple position: " << std::endl
00820             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
00821             << "      vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
00822 
00823   // return vertex
00824   // =============
00825   return newVertex;
00826 }
00827 
00828 WCSimRecoVertex* WCSimVertexFinder::FixPointPosition(Double_t vtxX, Double_t vtxY, Double_t vtxZ)
00829 {  
00830   // initialization
00831   // ==============
00832   Double_t vtxTime = 950.0;
00833   Double_t vtxFOM = 0.0;
00834 
00835   // create new vertex
00836   // =================
00837   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00838   vVertexList.push_back(newVertex);
00839   fPointPosition = newVertex;
00840 
00841   // calculate vertex
00842   // ================
00843   this->PointPositionChi2(vtxX,vtxY,vtxZ,
00844                           vtxTime,vtxFOM);
00845 
00846   // set vertex
00847   // ==========
00848   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00849   newVertex->SetFOM(vtxFOM,1,1);
00850 
00851   // set status
00852   // ==========
00853   newVertex->SetStatus(WCSimRecoVertex::kOK);
00854 
00855   // print vertex
00856   // ============
00857   std::cout << "  set point position: " << std::endl
00858             <<"     (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
00859             << "      vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
00860 
00861   // return vertex
00862   // =============
00863   return newVertex;
00864 }
00865 
00866 WCSimRecoVertex* WCSimVertexFinder::FixSimpleDirection(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00867 {
00868   // initialization
00869   // ==============
00870   Double_t vtxTime = 950.0;
00871   Double_t dirFOM = fBaseFOM;
00872 
00873   // create new vertex
00874   // =================
00875   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00876   vVertexList.push_back(newVertex);
00877   fSimpleDirection = newVertex; 
00878 
00879   // set vertex
00880   // ==========
00881   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00882   newVertex->SetDirection(dirX,dirY,dirZ);  
00883   newVertex->SetFOM(dirFOM,1,1);
00884 
00885   // set status
00886   // ==========
00887   newVertex->SetStatus(WCSimRecoVertex::kOK);
00888 
00889   // print vertex
00890   // ============
00891   std::cout << "  set simple direction: " << std::endl
00892             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
00893             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
00894             << "      vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
00895 
00896   // return vertex
00897   // =============
00898   return newVertex;
00899 }
00900 
00901 WCSimRecoVertex* WCSimVertexFinder::FixPointDirection(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00902 {
00903   // initialization
00904   // ==============
00905   Double_t vtxTime = 950.0;
00906   Double_t vtxAngle = 42.0;
00907   Double_t vtxFOM = 0.0;
00908   Double_t dirFOM = 0.0;
00909 
00910   // create new vertex
00911   // =================
00912   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00913   vVertexList.push_back(newVertex);
00914   fPointDirection = newVertex;
00915 
00916   // figure of merit
00917   // ===============
00918   this->PointPositionChi2(vtxX,vtxY,vtxZ,
00919                           vtxTime,vtxFOM);
00920 
00921   this->PointDirectionChi2(vtxX,vtxY,vtxZ,
00922                            dirX,dirY,dirZ,
00923                            vtxAngle,dirFOM);
00924 
00925   // set vertex and direction
00926   // ========================
00927   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00928   newVertex->SetDirection(dirX,dirY,dirZ);
00929   newVertex->SetConeAngle(vtxAngle);
00930   newVertex->SetFOM(dirFOM,1,1);
00931 
00932   // set status
00933   // ==========
00934   newVertex->SetStatus(WCSimRecoVertex::kOK);
00935 
00936   // print vertex
00937   // ============
00938   std::cout << "  set point direction: " << std::endl
00939             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
00940             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
00941             << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
00942 
00943   // return vertex
00944   // =============
00945   return newVertex;
00946 }
00947 
00948 WCSimRecoVertex* WCSimVertexFinder::FixPointVertex(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00949 {
00950   // initialization
00951   // ==============
00952   Double_t vtxTime = 950.0;
00953   Double_t vtxAngle = 42.0;
00954   Double_t vtxFOM = 0.0;
00955 
00956   // create new vertex
00957   // =================
00958   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00959   vVertexList.push_back(newVertex);
00960   fPointVertex = newVertex;
00961 
00962   // figure of merit
00963   // ===============
00964   this->PointVertexChi2(vtxX,vtxY,vtxZ,
00965                         dirX,dirY,dirZ,
00966                         vtxAngle,vtxTime,vtxFOM);
00967 
00968   // set vertex and direction
00969   // ========================
00970   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00971   newVertex->SetDirection(dirX,dirY,dirZ);
00972   newVertex->SetConeAngle(vtxAngle);
00973   newVertex->SetFOM(vtxFOM,1,1);
00974 
00975   // set status
00976   // ==========
00977   newVertex->SetStatus(WCSimRecoVertex::kOK);
00978 
00979   // print vertex
00980   // ============
00981   std::cout << "  set point vertex: " << std::endl
00982             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
00983             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
00984             << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
00985   
00986   // return vertex
00987   // =============
00988   return newVertex;
00989 }
00990 
00991 WCSimRecoVertex* WCSimVertexFinder::FixExtendedVertex(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
00992 {
00993   // initialization
00994   // ==============
00995   Double_t vtxTime = 950.0;
00996   Double_t vtxAngle = 42.0;
00997   Double_t vtxFOM = 0.0;
00998 
00999   // create new vertex
01000   // =================
01001   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
01002   vVertexList.push_back(newVertex);
01003   fExtendedVertex = newVertex;
01004 
01005   // figure of merit
01006   // ===============
01007   this->ExtendedVertexChi2(vtxX,vtxY,vtxZ,
01008                            dirX,dirY,dirZ,
01009                            vtxAngle,vtxTime,vtxFOM); 
01010 
01011 
01012   // set vertex and direction
01013   // ========================
01014   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
01015   newVertex->SetDirection(dirX,dirY,dirZ);
01016   newVertex->SetConeAngle(vtxAngle);
01017   newVertex->SetFOM(vtxFOM,1,1);
01018 
01019   // set status
01020   // ==========
01021   newVertex->SetStatus(WCSimRecoVertex::kOK);
01022 
01023   // print vertex
01024   // ============
01025   std::cout << "  set extended vertex: " << std::endl
01026             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
01027             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
01028             << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
01029   
01030   // return vertex
01031   // =============
01032   return newVertex;
01033 }
01034 
01035 WCSimRecoVertex* WCSimVertexFinder::CorrectExtendedVertex(WCSimRecoVertex* myVertex)
01036 {
01037   // get vertex and direction
01038   // ========================
01039   Double_t vtxX = myVertex->GetX();
01040   Double_t vtxY = myVertex->GetY();
01041   Double_t vtxZ = myVertex->GetZ();
01042   Double_t vtxTime = myVertex->GetTime();
01043 
01044   Double_t dirX = myVertex->GetDirX();
01045   Double_t dirY = myVertex->GetDirY();
01046   Double_t dirZ = myVertex->GetDirZ();
01047 
01048   // get other vertex parameters
01049   // ===========================
01050   Double_t angle = myVertex->GetConeAngle();
01051   Double_t length = myVertex->GetTrackLength();
01052 
01053   Double_t fom = myVertex->GetFOM();
01054   Int_t nsteps = myVertex->GetIterations();
01055   Bool_t pass = myVertex->GetPass();
01056   Int_t status = myVertex->GetStatus();
01057 
01058   // fix vertex bias
01059   // ===============
01060   if( fFixVertexBias ){
01061     vtxX -= fVertexBias*dirX;
01062     vtxY -= fVertexBias*dirY;
01063     vtxZ -= fVertexBias*dirZ;
01064     vtxTime -= fVertexBias/WCSimParameters::SpeedOfLight();
01065   }
01066 
01067   // create new vertex
01068   // =================
01069   WCSimRecoVertex* newVertex = new WCSimRecoVertex(vtxX,vtxY,vtxZ,vtxTime,
01070                                                    dirX,dirY,dirZ,
01071                                                    angle,length,
01072                                                    fom,nsteps,pass,status);
01073   vVertexList.push_back(newVertex);
01074   fExtendedVertex = newVertex;
01075 
01076   // print vertex
01077   // ============
01078   if( fFixVertexBias ){
01079     std::cout << "  corrected extended vertex: " << std::endl
01080               << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
01081               << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
01082               << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
01083   }
01084 
01085   // return vertex
01086   // =============
01087   return newVertex;
01088 }
01089 
01090 WCSimRecoVertex* WCSimVertexFinder::FitPointPositionWithMinuit(WCSimRecoVertex* myVertex)
01091 {  
01092   // initialization
01093   // ==============
01094   Double_t vtxX = 0.0;
01095   Double_t vtxY = 0.0;
01096   Double_t vtxZ = 0.0;
01097   Double_t vtxTime = 950.0;
01098 
01099   Double_t dirX = 0.0;
01100   Double_t dirY = 0.0;
01101   Double_t dirZ = 0.0;
01102 
01103   Double_t vtxFOM = 0.0;
01104 
01105   // seed vertex
01106   // ===========
01107   Bool_t foundSeed = myVertex->FoundVertex();
01108   Double_t seedX = myVertex->GetX();
01109   Double_t seedY = myVertex->GetY();
01110   Double_t seedZ = myVertex->GetZ();
01111 
01112   // current status
01113   // ==============
01114   Int_t status = myVertex->GetStatus();
01115 
01116   // reset counter
01117   // =============
01118   point_position_reset_itr();
01119 
01120   // create new vertex
01121   // =================
01122   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
01123   vVertexList.push_back(newVertex);
01124   fPointPosition = newVertex;
01125 
01126   // abort if necessary
01127   // ==================
01128   if( foundSeed==0 ){
01129     std::cout << "   <warning> point position fit failed to find input vertex " << std::endl;    
01130     status |= WCSimRecoVertex::kFailPointPosition;
01131     newVertex->SetStatus(status);
01132     return newVertex;
01133   }
01134 
01135   // run Minuit
01136   // ==========  
01137   // three-parameter fit to vertex coordinates
01138 
01139   Int_t err = 0;
01140   Int_t flag = 0;
01141 
01142   Double_t fitXpos = 0.0;
01143   Double_t fitYpos = 0.0;
01144   Double_t fitZpos = 0.0;
01145 
01146   Double_t fitXposErr = 0.0;
01147   Double_t fitYposErr = 0.0;
01148   Double_t fitZposErr = 0.0;
01149 
01150   Double_t* arglist = new Double_t[10];
01151   arglist[0]=1;  // 1: standard minimization
01152                  // 2: try to improve minimum
01153 
01154   // re-initialize everything...
01155   fMinuitPointPosition->mncler();
01156   fMinuitPointPosition->SetFCN(point_position_chi2);
01157   fMinuitPointPosition->mnexcm("SET STR",arglist,1,err);
01158   fMinuitPointPosition->mnparm(0,"x",seedX,250.0,-5000.0,+5000.0,err);
01159   fMinuitPointPosition->mnparm(1,"y",seedY,250.0,-5000.0,+5000.0,err);
01160   fMinuitPointPosition->mnparm(2,"z",seedZ,250.0,-5000.0,+5000.0,err);
01161 
01162   flag = fMinuitPointPosition->Migrad();
01163   fMinuitPointPosition->GetParameter(0,fitXpos,fitXposErr);
01164   fMinuitPointPosition->GetParameter(1,fitYpos,fitYposErr);
01165   fMinuitPointPosition->GetParameter(2,fitZpos,fitZposErr);
01166 
01167   delete [] arglist;
01168 
01169   // sort results
01170   // ============
01171   vtxX = fitXpos; 
01172   vtxY = fitYpos;
01173   vtxZ = fitZpos;
01174   vtxTime = 950.0;
01175 
01176   vtxFOM = 0.0;
01177   
01178   fPass = 0;               // flag = 0: normal termination
01179   if( flag==0 ) fPass = 1; // anything else: abnormal termination 
01180 
01181   fItr = point_position_iterations();
01182 
01183   // calculate vertex
01184   // ================
01185   this->PointPositionChi2(vtxX,vtxY,vtxZ,
01186                           vtxTime,vtxFOM);
01187 
01188   // set vertex and direction
01189   // ========================
01190   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
01191   newVertex->SetDirection(dirX,dirY,dirZ);
01192   newVertex->SetFOM(vtxFOM,fItr,fPass);
01193 
01194   // set status
01195   // ==========
01196   if( !fPass ) status |= WCSimRecoVertex::kFailPointPosition;
01197   newVertex->SetStatus(status);
01198 
01199   // print vertex
01200   // ============
01201   std::cout << "  fitted point position: " << std::endl
01202             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
01203             << "      vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
01204 
01205   if( !fPass ) std::cout << "   <warning> point position fit failed to converge! Error code: " << flag << std::endl;
01206 
01207   // return vertex
01208   // =============  
01209   return newVertex;
01210 }
01211 
01212 WCSimRecoVertex* WCSimVertexFinder::FitPointDirectionWithMinuit(WCSimRecoVertex* myVertex)
01213 {
01214   // initialization
01215   // ==============
01216   Bool_t foundSeed = ( myVertex->FoundVertex() 
01217                     && myVertex->FoundDirection() );
01218   
01219   Double_t vtxX = myVertex->GetX();
01220   Double_t vtxY = myVertex->GetY();
01221   Double_t vtxZ = myVertex->GetZ();
01222   Double_t vtxTime = myVertex->GetTime();
01223 
01224   Double_t dirX = 0.0;
01225   Double_t dirY = 0.0;
01226   Double_t dirZ = 0.0;  
01227 
01228   Double_t vtxAngle = 42.0;
01229 
01230   Double_t vtxFOM = 0.0;
01231 
01232   // seed direction
01233   // ==============
01234   Double_t seedDirX = myVertex->GetDirX();
01235   Double_t seedDirY = myVertex->GetDirY();
01236   Double_t seedDirZ = myVertex->GetDirZ();
01237   
01238   Double_t seedTheta = acos(seedDirZ);
01239   Double_t seedPhi = 0.0;
01240 
01241   if( seedDirX!=0.0 ){
01242     seedPhi = atan(seedDirY/seedDirX);
01243   }
01244   if( seedDirX<=0.0 ){
01245     if( seedDirY>0.0 ) seedPhi += TMath::Pi();
01246     if( seedDirY<0.0 ) seedPhi -= TMath::Pi();
01247   }  
01248 
01249   // current status
01250   // ==============
01251   Int_t status = myVertex->GetStatus();
01252 
01253   // set global vertex
01254   // =================
01255   fVtxX = vtxX;
01256   fVtxY = vtxY;
01257   fVtxZ = vtxZ;
01258 
01259   // reset counter
01260   // =============
01261   point_direction_reset_itr();
01262 
01263   // create new vertex
01264   // =================
01265   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
01266   vVertexList.push_back(newVertex);
01267   fPointDirection = newVertex;
01268 
01269   // abort if necessary
01270   // ==================
01271   if( foundSeed==0 ){
01272     std::cout << "   <warning> point direction fit failed to find input vertex " << std::endl;
01273     status |= WCSimRecoVertex::kFailPointDirection;
01274     newVertex->SetStatus(status);
01275     return newVertex;
01276   }
01277 
01278   // run Minuit
01279   // ==========  
01280   // two-parameter fit to direction coordinates
01281   
01282   Int_t err = 0;
01283   Int_t flag = 0;
01284 
01285   Double_t dirTheta;
01286   Double_t dirPhi;
01287 
01288   Double_t dirThetaErr;
01289   Double_t dirPhiErr;
01290 
01291   Double_t* arglist = new Double_t[10];
01292   arglist[0]=1;  // 1: standard minimization
01293                  // 2: try to improve minimum
01294 
01295   // re-initialize everything...
01296   fMinuitPointDirection->mncler();
01297   fMinuitPointDirection->SetFCN(point_direction_chi2);
01298   fMinuitPointDirection->mnexcm("SET STR",arglist,1,err);
01299   fMinuitPointDirection->mnparm(0,"theta",seedTheta,0.125*TMath::Pi(),0.0,TMath::Pi(),err);
01300   fMinuitPointDirection->mnparm(1,"phi",seedPhi,0.25*TMath::Pi(),-1.0*TMath::Pi(),+3.0*TMath::Pi(),err);
01301 
01302   flag = fMinuitPointDirection->Migrad();
01303   fMinuitPointDirection->GetParameter(0,dirTheta,dirThetaErr);
01304   fMinuitPointDirection->GetParameter(1,dirPhi,dirPhiErr);
01305 
01306   delete [] arglist;
01307 
01308   // sort results
01309   // ============
01310   dirX = sin(dirTheta)*cos(dirPhi);
01311   dirY = sin(dirTheta)*sin(dirPhi);
01312   dirZ = cos(dirTheta);
01313 
01314   vtxFOM = 0.0;
01315   
01316   fPass = 0;               // flag = 0: normal termination
01317   if( flag==0 ) fPass = 1; // anything else: abnormal termination 
01318                            
01319   fItr = point_direction_iterations();
01320 
01321   // calculate vertex
01322   // ================
01323   this->PointDirectionChi2(vtxX,vtxY,vtxZ,
01324                            dirX,dirY,dirZ,
01325                            vtxAngle,vtxFOM);
01326 
01327   // set vertex and direction
01328   // ========================
01329   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
01330   newVertex->SetDirection(dirX,dirY,dirZ);
01331   newVertex->SetConeAngle(vtxAngle);
01332   newVertex->SetFOM(vtxFOM,fItr,fPass);
01333 
01334   // set status
01335   // ==========
01336   if( !fPass ) status |= WCSimRecoVertex::kFailPointDirection;
01337   newVertex->SetStatus(status);
01338 
01339   // print vertex
01340   // ============
01341   std::cout << "  fitted point direction: " << std::endl
01342             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
01343             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
01344             << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
01345 
01346   if( !fPass ) std::cout << "   <warning> point direction fit failed to converge! Error code: " << flag << std::endl;
01347 
01348   // return vertex
01349   // ============= 
01350   return newVertex;
01351 }
01352 
01353 WCSimRecoVertex* WCSimVertexFinder::FitPointVertexWithMinuit(WCSimRecoVertex* myVertex)
01354 {
01355   // initialization
01356   // ==============
01357   Double_t vtxX = 0.0;
01358   Double_t vtxY = 0.0;
01359   Double_t vtxZ = 0.0;
01360   Double_t vtxTime = 950.0;
01361 
01362   Double_t dirX = 0.0;
01363   Double_t dirY = 0.0;
01364   Double_t dirZ = 0.0;
01365     
01366   Double_t vtxAngle = 42.0;
01367 
01368   Double_t vtxFOM = 0.0;
01369 
01370   // seed vertex
01371   // ===========  
01372   Bool_t foundSeed = ( myVertex->FoundVertex() 
01373                     && myVertex->FoundDirection() );
01374 
01375   Double_t seedX = myVertex->GetX();
01376   Double_t seedY = myVertex->GetY();
01377   Double_t seedZ = myVertex->GetZ();
01378 
01379   Double_t seedDirX = myVertex->GetDirX();
01380   Double_t seedDirY = myVertex->GetDirY();
01381   Double_t seedDirZ = myVertex->GetDirZ();
01382 
01383   Double_t seedTheta = acos(seedDirZ);
01384   Double_t seedPhi = 0.0;
01385 
01386   if( seedDirX!=0.0 ){
01387     seedPhi = atan(seedDirY/seedDirX);
01388   }
01389   if( seedDirX<=0.0 ){
01390     if( seedDirY>0.0 ) seedPhi += TMath::Pi();
01391     if( seedDirY<0.0 ) seedPhi -= TMath::Pi();
01392   }  
01393 
01394   // current status
01395   // ==============
01396   Int_t status = myVertex->GetStatus();
01397 
01398   // reset counter
01399   // =============
01400   point_vertex_reset_itr();
01401 
01402   // create new vertex
01403   // =================
01404   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
01405   vVertexList.push_back(newVertex);
01406   fPointVertex = newVertex;
01407 
01408   // abort if necessary
01409   // ==================
01410   if( foundSeed==0 ){
01411     std::cout << "   <warning> point vertex fit failed to find input vertex " << std::endl;   
01412     status |= WCSimRecoVertex::kFailPointVertex;
01413     newVertex->SetStatus(status);
01414     return newVertex;
01415   }
01416 
01417   // run Minuit
01418   // ==========  
01419   // five-parameter fit to vertex and direction
01420 
01421   Int_t err = 0;
01422   Int_t flag = 0;
01423 
01424   Double_t fitXpos = 0.0;
01425   Double_t fitYpos = 0.0;
01426   Double_t fitZpos = 0.0;
01427   Double_t fitTheta = 0.0;
01428   Double_t fitPhi = 0.0;
01429 
01430   Double_t fitXposErr = 0.0;
01431   Double_t fitYposErr = 0.0;
01432   Double_t fitZposErr = 0.0;
01433   Double_t fitThetaErr = 0.0;
01434   Double_t fitPhiErr = 0.0;
01435 
01436   Double_t* arglist = new Double_t[10];
01437   arglist[0]=2;  // 1: standard minimization
01438                  // 2: try to improve minimum
01439 
01440   // re-initialize everything...
01441   fMinuitPointVertex->mncler();
01442   fMinuitPointVertex->SetFCN(point_vertex_chi2);
01443   fMinuitPointVertex->mnexcm("SET STR",arglist,1,err);
01444   fMinuitPointVertex->mnparm(0,"x",seedX,250.0,-5000.0,+5000.0,err);
01445   fMinuitPointVertex->mnparm(1,"y",seedY,250.0,-5000.0,+5000.0,err);
01446   fMinuitPointVertex->mnparm(2,"z",seedZ,250.0,-5000.0,+5000.0,err);
01447   fMinuitPointVertex->mnparm(3,"theta",seedTheta,0.125*TMath::Pi(),0.0,TMath::Pi(),err);
01448   fMinuitPointVertex->mnparm(4,"phi",seedPhi,0.25*TMath::Pi(),-1.0*TMath::Pi(),+3.0*TMath::Pi(),err);
01449   
01450   flag = fMinuitPointVertex->Migrad();
01451   fMinuitPointVertex->GetParameter(0,fitXpos,fitXposErr);
01452   fMinuitPointVertex->GetParameter(1,fitYpos,fitYposErr);
01453   fMinuitPointVertex->GetParameter(2,fitZpos,fitZposErr);
01454   fMinuitPointVertex->GetParameter(3,fitTheta,fitThetaErr);
01455   fMinuitPointVertex->GetParameter(4,fitPhi,fitPhiErr);
01456 
01457   delete [] arglist;
01458 
01459   // sort results
01460   // ============
01461   vtxX = fitXpos; 
01462   vtxY = fitYpos;
01463   vtxZ = fitZpos;
01464   vtxTime = 950.0;
01465 
01466   dirX = sin(fitTheta)*cos(fitPhi);
01467   dirY = sin(fitTheta)*sin(fitPhi);
01468   dirZ = cos(fitTheta);  
01469 
01470   vtxFOM = 0.0;
01471   
01472   fPass = 0;               // flag = 0: normal termination
01473   if( flag==0 ) fPass = 1; // anything else: abnormal termination 
01474 
01475   fItr = point_vertex_iterations();
01476 
01477   // calculate vertex
01478   // ================
01479   this->PointVertexChi2(vtxX,vtxY,vtxZ,
01480                         dirX,dirY,dirZ, 
01481                         vtxAngle,vtxTime,vtxFOM);
01482 
01483   // set vertex and direction
01484   // ========================
01485   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
01486   newVertex->SetDirection(dirX,dirY,dirZ);
01487   newVertex->SetConeAngle(vtxAngle);
01488   newVertex->SetFOM(vtxFOM,fItr,fPass);
01489 
01490   // set status
01491   // ==========
01492   if( !fPass ) status |= WCSimRecoVertex::kFailPointVertex;
01493   newVertex->SetStatus(status);
01494 
01495   // print vertex
01496   // ============
01497   std::cout << "  fitted point vertex: " << std::endl
01498             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
01499             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
01500             << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
01501 
01502   if( !fPass ) std::cout << "   <warning> point vertex fit failed to converge! Error code: " << flag << std::endl;
01503 
01504   // return vertex
01505   // =============  
01506   return newVertex;
01507 }
01508 
01509 WCSimRecoVertex* WCSimVertexFinder::FitExtendedVertexWithMinuit(WCSimRecoVertex* myVertex)
01510 {
01511   // initialization
01512   // ==============
01513   Double_t vtxX = 0.0;
01514   Double_t vtxY = 0.0;
01515   Double_t vtxZ = 0.0;
01516   Double_t vtxTime = 950.0;
01517 
01518   Double_t dirX = 0.0;
01519   Double_t dirY = 0.0;
01520   Double_t dirZ = 0.0;
01521   
01522   Double_t vtxAngle = 42.0;
01523 
01524   Double_t vtxFOM = 0.0;
01525 
01526   // seed vertex
01527   // ===========  
01528   Bool_t foundSeed = ( myVertex->FoundVertex() 
01529                     && myVertex->FoundDirection() );
01530 
01531   Double_t seedX = myVertex->GetX();
01532   Double_t seedY = myVertex->GetY();
01533   Double_t seedZ = myVertex->GetZ();
01534 
01535   Double_t seedDirX = myVertex->GetDirX();
01536   Double_t seedDirY = myVertex->GetDirY();
01537   Double_t seedDirZ = myVertex->GetDirZ();
01538 
01539   Double_t seedTheta = acos(seedDirZ);
01540   Double_t seedPhi = 0.0;
01541 
01542   if( seedDirX!=0.0 ){
01543     seedPhi = atan(seedDirY/seedDirX);
01544   }
01545   if( seedDirX<=0.0 ){
01546     if( seedDirY>0.0 ) seedPhi += TMath::Pi();
01547     if( seedDirY<0.0 ) seedPhi -= TMath::Pi();
01548   }  
01549 
01550   // current status
01551   // ==============
01552   Int_t status = myVertex->GetStatus();
01553 
01554   // reset counter
01555   // =============
01556   extended_vertex_reset_itr();
01557 
01558   // create new vertex
01559   // =================
01560   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
01561   vVertexList.push_back(newVertex);
01562   fExtendedVertex = newVertex;
01563 
01564   // abort if necessary
01565   // ==================
01566   if( foundSeed==0 ){
01567     std::cout << "   <warning> extended vertex fit failed to find input vertex " << std::endl;   
01568     status |= WCSimRecoVertex::kFailExtendedVertex;
01569     newVertex->SetStatus(status);
01570     return newVertex;
01571   }
01572 
01573   // run Minuit
01574   // ==========  
01575   // five-parameter fit to vertex and direction
01576 
01577   Int_t err = 0;
01578   Int_t flag = 0;
01579 
01580   Double_t fitXpos = 0.0;
01581   Double_t fitYpos = 0.0;
01582   Double_t fitZpos = 0.0;
01583   Double_t fitTheta = 0.0;
01584   Double_t fitPhi = 0.0;
01585 
01586   Double_t fitXposErr = 0.0;
01587   Double_t fitYposErr = 0.0;
01588   Double_t fitZposErr = 0.0;
01589   Double_t fitThetaErr = 0.0;
01590   Double_t fitPhiErr = 0.0;
01591 
01592   Double_t* arglist = new Double_t[10];
01593   arglist[0]=2;  // 1: standard minimization
01594                  // 2: try to improve minimum
01595 
01596   // re-initialize everything...
01597   fMinuitExtendedVertex->mncler();
01598   fMinuitExtendedVertex->SetFCN(extended_vertex_chi2);
01599   fMinuitExtendedVertex->mnexcm("SET STR",arglist,1,err);
01600   fMinuitExtendedVertex->mnparm(0,"x",seedX,250.0,-5000.0,+5000.0,err);
01601   fMinuitExtendedVertex->mnparm(1,"y",seedY,250.0,-5000.0,+5000.0,err);
01602   fMinuitExtendedVertex->mnparm(2,"z",seedZ,250.0,-5000.0,+5000.0,err);
01603   fMinuitExtendedVertex->mnparm(3,"theta",seedTheta,0.125*TMath::Pi(),0.0,TMath::Pi(),err);
01604   fMinuitExtendedVertex->mnparm(4,"phi",seedPhi,0.25*TMath::Pi(),-1.0*TMath::Pi(),+3.0*TMath::Pi(),err);
01605   
01606   flag = fMinuitExtendedVertex->Migrad();
01607   fMinuitExtendedVertex->GetParameter(0,fitXpos,fitXposErr);
01608   fMinuitExtendedVertex->GetParameter(1,fitYpos,fitYposErr);
01609   fMinuitExtendedVertex->GetParameter(2,fitZpos,fitZposErr);
01610   fMinuitExtendedVertex->GetParameter(3,fitTheta,fitThetaErr);
01611   fMinuitExtendedVertex->GetParameter(4,fitPhi,fitPhiErr);
01612 
01613   delete [] arglist;
01614 
01615   // sort results
01616   // ============
01617   vtxX = fitXpos; 
01618   vtxY = fitYpos;
01619   vtxZ = fitZpos;
01620   vtxTime = 950.0;
01621 
01622   dirX = sin(fitTheta)*cos(fitPhi);
01623   dirY = sin(fitTheta)*sin(fitPhi);
01624   dirZ = cos(fitTheta);  
01625 
01626   vtxFOM = 0.0;
01627   
01628   fPass = 0;               // flag = 0: normal termination
01629   if( flag==0 ) fPass = 1; // anything else: abnormal termination 
01630 
01631   fItr = extended_vertex_iterations();
01632 
01633   // calculate vertex
01634   // ================
01635   this->ExtendedVertexChi2(vtxX,vtxY,vtxZ,
01636                            dirX,dirY,dirZ, 
01637                            vtxAngle,vtxTime,vtxFOM);
01638 
01639   // set vertex and direction
01640   // ========================
01641   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
01642   newVertex->SetDirection(dirX,dirY,dirZ);
01643   newVertex->SetConeAngle(vtxAngle);
01644   newVertex->SetFOM(vtxFOM,fItr,fPass);
01645 
01646   // set status
01647   // ==========
01648   if( !fPass ) status |= WCSimRecoVertex::kFailExtendedVertex;
01649   newVertex->SetStatus(status);
01650 
01651   // print vertex
01652   // ============
01653   std::cout << "  fitted extended vertex: " << std::endl
01654             << "    (vx,vy,vz)=(" << newVertex->GetX() << "," << newVertex->GetY() << "," << newVertex->GetZ() << ") " << std::endl
01655             << "    (px,py,pz)=(" << newVertex->GetDirX() << "," << newVertex->GetDirY() << "," << newVertex->GetDirZ() << ") " << std::endl
01656             << "      angle=" << newVertex->GetConeAngle() << " vtime=" << newVertex->GetTime() << " itr=" << newVertex->GetIterations() << " fom=" << newVertex->GetFOM() << std::endl;
01657 
01658   if( !fPass ) std::cout << "   <warning> extended vertex fit failed to converge! Error code: " << flag << std::endl;
01659 
01660   // return vertex
01661   // =============  
01662   return newVertex;
01663 }
01664   
01665 void WCSimVertexFinder::PointPositionChi2(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t& vtxTime, Double_t& fom)
01666 {  
01667   // figure of merit
01668   // ===============
01669   Double_t vtxFOM = 0.0;
01670   Double_t penaltyFOM = 0.0;
01671   Double_t fixPositionFOM = 0.0;
01672 
01673   // calculate residuals
01674   // ===================
01675   this->PointResiduals(vtxX,vtxY,vtxZ);
01676 
01677   // calculate figure of merit
01678   // =========================
01679   this->PointPositionChi2(vtxTime,vtxFOM);
01680 
01681   // calculate penalty
01682   // =================
01683   this->PenaltyChi2(vtxX,vtxY,vtxZ,penaltyFOM);
01684 
01685   // fix true position
01686   // =================
01687   if( fUseTruePosition ){
01688     this->FixPositionChi2(vtxX,vtxY,vtxZ,fixPositionFOM);
01689   }
01690 
01691   // calculate overall figure of merit
01692   // =================================
01693   fom = vtxFOM + penaltyFOM + fixPositionFOM;
01694 
01695   // truncate
01696   if( fom<-999.999*fBaseFOM ) fom = -999.999*fBaseFOM;
01697 
01698   return;
01699 }
01700 
01701 void WCSimVertexFinder::PointDirectionChi2(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ, Double_t& vtxAngle, Double_t& fom)
01702 {  
01703   // figure of merit
01704   // ===============
01705   Double_t vtxFOM = 0.0;
01706   Double_t fixPositionFOM = 0.0;
01707   Double_t fixDirectionFOM = 0.0;
01708   
01709   // calculate residuals
01710   // ===================
01711   this->PointResiduals(vtxX,vtxY,vtxZ,
01712                        dirX,dirY,dirZ);
01713 
01714   // calculate figure of merit
01715   // =========================
01716   this->PointDirectionChi2(vtxAngle,vtxFOM);
01717 
01718   // fix true position
01719   // =================
01720   if( fUseTruePosition ){
01721     this->FixPositionChi2(vtxX,vtxY,vtxZ,fixPositionFOM);
01722   }
01723 
01724   // fix true direction
01725   // ==================
01726   if( fUseTrueDirection ){
01727     this->FixDirectionChi2(dirX,dirY,dirZ,fixDirectionFOM);
01728   }
01729 
01730   // calculate overall figure of merit
01731   // =================================
01732   fom = vtxFOM + fixDirectionFOM;
01733 
01734   // truncate
01735   if( fom<-999.999*fBaseFOM ) fom = -999.999*fBaseFOM;
01736 
01737   return;
01738 }
01739 
01740 void WCSimVertexFinder::PointVertexChi2(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ, Double_t& vtxAngle, Double_t& vtxTime, Double_t& fom)
01741 {  
01742   // figure of merit
01743   // ===============
01744   Double_t vtxFOM = 0.0;
01745   Double_t penaltyFOM = 0.0;
01746   Double_t fixPositionFOM = 0.0;
01747   Double_t fixDirectionFOM = 0.0;
01748 
01749   // calculate residuals
01750   // ===================
01751   this->PointResiduals(vtxX,vtxY,vtxZ,
01752                        dirX,dirY,dirZ);
01753 
01754   // calculate figure of merit
01755   // =========================
01756   this->PointVertexChi2(vtxAngle,vtxTime,vtxFOM);
01757 
01758   // calculate penalty
01759   // =================
01760   this->PenaltyChi2(vtxX,vtxY,vtxZ,penaltyFOM);
01761 
01762   // fix true position
01763   // =================
01764   if( fUseTruePosition ){
01765     this->FixPositionChi2(vtxX,vtxY,vtxZ,fixPositionFOM);
01766   }
01767 
01768   // fix true direction
01769   // ==================
01770   if( fUseTrueDirection ){
01771     this->FixDirectionChi2(dirX,dirY,dirZ,fixDirectionFOM);
01772   }
01773 
01774   // calculate overall figure of merit
01775   // =================================
01776   fom = vtxFOM + penaltyFOM + fixPositionFOM + fixDirectionFOM;
01777 
01778   // truncate
01779   if( fom<-999.999*fBaseFOM ) fom = -999.999*fBaseFOM;
01780 
01781   return;
01782 }
01783 
01784 void WCSimVertexFinder::ExtendedVertexChi2(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ, Double_t& vtxAngle, Double_t& vtxTime, Double_t& fom)
01785 {  
01786   // figure of merit
01787   // ===============
01788   Double_t vtxFOM = 0.0;
01789   Double_t penaltyFOM = 0.0;
01790   Double_t fixPositionFOM = 0.0;
01791   Double_t fixDirectionFOM = 0.0;
01792 
01793   // calculate residuals
01794   // ===================
01795   this->ExtendedResiduals(vtxX,vtxY,vtxZ,
01796                           dirX,dirY,dirZ);
01797 
01798   // calculate figure of merit
01799   // =========================
01800   this->ExtendedVertexChi2(vtxAngle,vtxTime,vtxFOM);
01801 
01802   // calculate penalty
01803   // =================
01804   this->PenaltyChi2(vtxX,vtxY,vtxZ,penaltyFOM);
01805 
01806   // fix true position
01807   // =================
01808   if( fUseTruePosition ){
01809     this->FixPositionChi2(vtxX,vtxY,vtxZ,fixPositionFOM);
01810   }
01811 
01812   // fix true direction
01813   // ==================
01814   if( fUseTrueDirection ){
01815     this->FixDirectionChi2(dirX,dirY,dirZ,fixDirectionFOM);
01816   }
01817 
01818   // calculate overall figure of merit
01819   // =================================
01820   fom = vtxFOM + penaltyFOM + fixPositionFOM + fixDirectionFOM;
01821 
01822   // truncate
01823   if( fom<-999.999*fBaseFOM ) fom = -999.999*fBaseFOM;
01824 
01825   return;
01826 }
01827 
01828 void WCSimVertexFinder::PointPositionChi2(Double_t& vtxTime, Double_t& fom)
01829 {
01830   // calculate figure of merit
01831   // =========================
01832   this->FitPointTimePropertiesLnL(vtxTime,fom);
01833 
01834   return;
01835 }
01836 
01837 void WCSimVertexFinder::PointDirectionChi2(Double_t& vtxAngle, Double_t& fom)
01838 {
01839   // calculate figure of merit
01840   // =========================  
01841   this->FitPointConePropertiesLnL(vtxAngle,fom);
01842 
01843   return;
01844 }
01845 
01846 void WCSimVertexFinder::PointVertexChi2(Double_t& vtxAngle, Double_t& vtxTime, Double_t& fom)
01847 {
01848   // calculate figure of merit
01849   // =========================
01850   Double_t timeFOM = 0.0;
01851   Double_t coneFOM = 0.0;
01852 
01853   this->FitPointConePropertiesLnL(vtxAngle,coneFOM);
01854   this->FitPointTimePropertiesLnL(vtxTime,timeFOM);
01855   
01856   fom = (fTimeFitWeight*timeFOM+fConeFitWeight*coneFOM)/(fTimeFitWeight+fConeFitWeight);
01857 
01858   return;
01859 }
01860 
01861 void WCSimVertexFinder::ExtendedVertexChi2(Double_t& vtxAngle, Double_t& vtxTime, Double_t& fom)
01862 {
01863   // calculate figure of merit
01864   // =========================
01865   Double_t timeFOM = 0.0;
01866   Double_t coneFOM = 0.0;
01867 
01868   this->FitExtendedConePropertiesLnL(vtxAngle,coneFOM);
01869   this->FitExtendedTimePropertiesLnL(vtxTime,timeFOM);
01870   
01871   fom = (fTimeFitWeight*timeFOM+fConeFitWeight*coneFOM)/(fTimeFitWeight+fConeFitWeight);
01872 
01873   return;
01874 }
01875 
01876 void WCSimVertexFinder::PointResiduals(Double_t vtxX, Double_t vtxY, Double_t vtxZ)
01877 {
01878   return WCSimVertexGeometry::Instance()->CalcPointResiduals(vtxX,vtxY,vtxZ,0.0, 
01879                                                               0.0, 0.0, 0.0);
01880 }
01881 
01882 void WCSimVertexFinder::PointResiduals(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
01883 {  
01884   return WCSimVertexGeometry::Instance()->CalcPointResiduals(vtxX,vtxY,vtxZ,0.0,
01885                                                               dirX,dirY,dirZ);
01886 }
01887 
01888 void WCSimVertexFinder::ExtendedResiduals(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t dirX, Double_t dirY, Double_t dirZ)
01889 {
01890   return WCSimVertexGeometry::Instance()->CalcExtendedResiduals(vtxX,vtxY,vtxZ,0.0,
01891                                                                  dirX,dirY,dirZ);
01892 }
01893 
01894 void WCSimVertexFinder::FitTimePropertiesFoM(Double_t& vtxTime, Double_t& vtxFOM)
01895 { 
01896   // calculate mean and rms
01897   // ====================== 
01898   Double_t meanTime = 950.0;
01899 
01900   this->FindSimpleTimeProperties(meanTime); 
01901 
01902   // reset counter
01903   // =============
01904   time_fit_reset_itr();
01905 
01906   // run Minuit
01907   // ==========  
01908   // one-parameter fit to vertex time
01909 
01910   Int_t err = 0;
01911   Int_t flag = 0;
01912 
01913   Double_t seedTime = meanTime;
01914 
01915   Double_t fitTime = 0.0;
01916   Double_t fitTimeErr = 0.0;
01917   
01918   Double_t* arglist = new Double_t[10];
01919   arglist[0]=1;  // 1: standard minimization
01920                  // 2: try to improve minimum
01921 
01922   // re-initialize everything...
01923   fMinuitTimeFit->mncler();
01924   fMinuitTimeFit->SetFCN(vertex_time_fom);
01925   fMinuitTimeFit->mnexcm("SET STR",arglist,1,err);
01926   fMinuitTimeFit->mnparm(0,"vtxTime",seedTime,50.0,0.0,2000.0,err);
01927   
01928   flag = fMinuitTimeFit->Migrad();
01929   fMinuitTimeFit->GetParameter(0,fitTime,fitTimeErr);
01930   
01931   delete [] arglist;
01932 
01933 
01934   // calculate figure of merit
01935   // =========================
01936   Double_t fom = -999.999*fBaseFOM;
01937 
01938   this->TimePropertiesFoM(fitTime,fom);  
01939 
01940 
01941   // return figure of merit
01942   // ======================
01943   vtxTime = fitTime;
01944   vtxFOM = fom;
01945 
01946   return;
01947 }
01948 
01949 void WCSimVertexFinder::FitPointTimePropertiesLnL(Double_t& vtxTime, Double_t& vtxFOM)
01950 { 
01951   // calculate mean and rms
01952   // ====================== 
01953   Double_t meanTime = 950.0;
01954 
01955   this->FindSimpleTimeProperties(meanTime); 
01956 
01957   // reset counter
01958   // =============
01959   time_fit_reset_itr();
01960 
01961   // run Minuit
01962   // ==========  
01963   // one-parameter fit to time profile
01964 
01965   Int_t err = 0;
01966   Int_t flag = 0;
01967 
01968   Double_t seedTime = meanTime;
01969 
01970   Double_t fitParam = fFixTimeParam0;
01971 
01972   Double_t fitTime = 0.0;
01973   Double_t fitTimeErr = 0.0;  
01974   
01975   Double_t* arglist = new Double_t[10];
01976   arglist[0]=1;  // 1: standard minimization
01977                  // 2: try to improve minimum
01978 
01979   // re-initialize everything...
01980   fMinuitTimeFit->mncler();
01981   fMinuitTimeFit->SetFCN(vertex_time_lnl);
01982   fMinuitTimeFit->mnexcm("SET STR",arglist,1,err);
01983   fMinuitTimeFit->mnparm(0,"vtxTime",seedTime,50.0,0.0,2000.0,err);
01984   
01985   flag = fMinuitTimeFit->Migrad();
01986   fMinuitTimeFit->GetParameter(0,fitTime,fitTimeErr);
01987    
01988   delete [] arglist;
01989 
01990 
01991   // calculate figure of merit
01992   // =========================
01993   Double_t fom = -999.999*fBaseFOM;
01994 
01995   this->TimePropertiesLnL(fitTime,fitParam,fom);  
01996 
01997   //
01998   // std::cout << "   TimeFit: itr=" << time_fit_iterations() << " seedTime=" << seedTime << " fitTime=" << fitTime << " fitParam=" << fitParam << " fom=" << fom << std::endl;
01999   //
02000 
02001   // --- debug ---
02002   fTimeParam0 = fitParam;
02003 
02004   // return figure of merit
02005   // ======================
02006   vtxTime = fitTime;
02007   vtxFOM = fom;
02008 
02009   return;
02010 }
02011 
02012 void WCSimVertexFinder::FitExtendedTimePropertiesLnL(Double_t& vtxTime, Double_t& vtxFOM)
02013 {   
02014   // return result from point fit
02015   // ============================
02016   if( fFitTimeParams==0 ){
02017     return this->FitPointTimePropertiesLnL(vtxTime,vtxFOM);
02018   }
02019 
02020   // calculate mean and rms
02021   // ====================== 
02022   Double_t meanTime = 950.0;
02023 
02024   this->FindSimpleTimeProperties(meanTime); 
02025 
02026   // reset counter
02027   // =============
02028   time_fit_reset_itr();
02029 
02030   // run Minuit
02031   // ==========  
02032   // two-parameter fit to time profile
02033 
02034   Int_t err = 0;
02035   Int_t flag = 0;
02036 
02037   Double_t seedTime = meanTime;
02038   Double_t seedParam = fFixTimeParam0;
02039 
02040   Double_t fitTime = seedTime;
02041   Double_t fitTimeErr = 0.0;  
02042 
02043   Double_t fitParam = seedParam;
02044   Double_t fitParamErr = 0.0;
02045   
02046   Double_t* arglist = new Double_t[10];
02047   arglist[0]=1;  // 1: standard minimization
02048                  // 2: try to improve minimum
02049 
02050   // re-initialize everything...
02051   fMinuitTimeFit->mncler();
02052   fMinuitTimeFit->SetFCN(vertex_time_lnl);
02053   fMinuitTimeFit->mnexcm("SET STR",arglist,1,err);
02054   fMinuitTimeFit->mnparm(0,"vtxTime",seedTime,50.0,0.0,2000.0,err);
02055   fMinuitTimeFit->mnparm(1,"vtxParam",seedParam,0.05,0.0,1.0,err);
02056   
02057   flag = fMinuitTimeFit->Migrad();
02058   fMinuitTimeFit->GetParameter(0,fitTime,fitTimeErr);
02059   fMinuitTimeFit->GetParameter(1,fitParam,fitParamErr);
02060   
02061   delete [] arglist;
02062 
02063 
02064   // calculate figure of merit
02065   // =========================
02066   Double_t fom = -999.999*fBaseFOM;
02067 
02068   this->TimePropertiesLnL(fitTime,fitParam,fom);  
02069 
02070   //
02071   // std::cout << "   TimeFit: itr=" << time_fit_iterations() << " seedTime=" << seedTime << " fitTime=" << fitTime << " fitParam=" << fitParam << " fom=" << fom << std::endl;
02072   //
02073 
02074   // --- debug ---
02075   fTimeParam0 = fitParam;
02076 
02077   // return figure of merit
02078   // ======================
02079   vtxTime = fitTime;
02080   vtxFOM = fom;
02081 
02082   return;
02083 }
02084 
02085 void WCSimVertexFinder::FitConePropertiesFoM(Double_t& coneAngle, Double_t& coneFOM)
02086 {
02087   coneAngle = 42.0; // nominal cone angle
02088 
02089   this->ConePropertiesFoM(coneFOM);
02090 
02091   return;
02092 }
02093 
02094 void WCSimVertexFinder::FitPointConePropertiesLnL(Double_t& coneAngle, Double_t& coneFOM)
02095 {  
02096   coneAngle  = 42.0; // nominal cone angle
02097     
02098   this->ConePropertiesLnL(fFixConeParam0,fFixConeParam1,fFixConeParam2,
02099                           coneAngle,coneFOM);  
02100 
02101   // --- debug ---
02102   fConeParam0 = fFixConeParam0;
02103   fConeParam1 = fFixConeParam1;
02104   fConeParam2 = fFixConeParam2;
02105 
02106   return;
02107 }
02108 
02109 void WCSimVertexFinder::FitExtendedConePropertiesLnL(Double_t& coneAngle, Double_t& coneFOM)
02110 {  
02111   // return result from point fit
02112   // ============================
02113   if( fFitConeParams==0 ){
02114     return this->FitPointConePropertiesLnL(coneAngle,coneFOM);
02115   }
02116 
02117   // reset counter
02118   // =============
02119   cone_fit_reset_itr();
02120 
02121   // run Minuit
02122   // ==========  
02123   // one-parameter fit to angular distribution
02124 
02125   Int_t err = 0;
02126   Int_t flag = 0;
02127 
02128   Double_t seedParam0 = fFixConeParam0;
02129   Double_t seedParam1 = fFixConeParam1;
02130   Double_t seedParam2 = fFixConeParam2;
02131 
02132   Double_t fitParam0 = seedParam0;
02133   Double_t fitParam1 = seedParam1;
02134   Double_t fitParam2 = seedParam2;
02135 
02136   Double_t fitParam0Err = 0.0;
02137   Double_t fitParam1Err = 0.0; 
02138   Double_t fitParam2Err = 0.0;
02139 
02140   Double_t fitAngle  = 42.0;
02141   
02142   Double_t* arglist = new Double_t[10];
02143   arglist[0]=1;  // 1: standard minimization
02144                  // 2: try to improve minimum
02145 
02146   // re-initialize everything...
02147   fMinuitConeFit->mncler();
02148   fMinuitConeFit->SetFCN(vertex_cone_lnl);
02149   fMinuitConeFit->mnexcm("SET STR",arglist,1,err);
02150   fMinuitConeFit->mnparm(0,"vtxParam0",seedParam0,0.25,0.0,1.0,err);
02151   fMinuitConeFit->mnparm(1,"vtxParam1",seedParam1,0.25,0.0,1.0,err);
02152   fMinuitConeFit->mnparm(2,"vtxParam2",seedParam2,0.25,0.0,1.0,err);
02153 
02154   flag = fMinuitConeFit->Migrad();
02155   fMinuitConeFit->GetParameter(0,fitParam0,fitParam0Err);
02156   fMinuitConeFit->GetParameter(1,fitParam1,fitParam1Err);
02157   fMinuitConeFit->GetParameter(2,fitParam2,fitParam2Err);
02158   
02159   delete [] arglist;
02160 
02161 
02162   // calculate figure of merit
02163   // =========================
02164   Double_t fom = -999.999*fBaseFOM;
02165   
02166   this->ConePropertiesLnL(fitParam0,fitParam1,fitParam2,
02167                           fitAngle,fom);  
02168 
02169   //
02170   // std::cout << "   ConeFit: itr=" << cone_fit_iterations() << " fitParam0=" << fitParam0 << " fitParam1=" << fitParam1 << " fitParam2=" << fitParam2 << " fom=" << fom << std::endl;
02171   // 
02172 
02173   // --- debug ---
02174   fConeParam0 = fitParam0;
02175   fConeParam1 = fitParam1;
02176   fConeParam2 = fitParam2;
02177 
02178   // return figure of merit
02179   // ======================
02180   coneAngle = fitAngle;
02181   coneFOM = fom;
02182 
02183   return;
02184 }
02185 
02186 void WCSimVertexFinder::FindSimpleTimeProperties(Double_t& vtxTime)
02187 {
02188   // reset mean and rms
02189   // ================== 
02190   Double_t meanTime = 950.0;
02191 
02192   // calculate mean and rms of hits inside cone
02193   // ==========================================
02194   Double_t Swx = 0.0;
02195   Double_t Sw = 0.0;
02196 
02197   Double_t delta = 0.0;
02198   Double_t sigma = 0.0;
02199   Double_t weight = 0.0;
02200   Double_t deweight = 0.0;
02201   Double_t deltaAngle = 0.0;
02202 
02203   Double_t myConeEdge = 42.0;      // [degrees]
02204   Double_t myConeEdgeSigma = 7.0;  // [degrees]
02205 
02206   for( Int_t idigit=0; idigit<WCSimVertexGeometry::Instance()->GetNDigits(); idigit++ ){
02207 
02208     if( WCSimVertexGeometry::Instance()->IsFiltered(idigit) ){
02209       delta = WCSimVertexGeometry::Instance()->GetDelta(idigit);
02210       sigma = WCSimVertexGeometry::Instance()->GetDeltaSigma(idigit);
02211       weight = 1.0/(sigma*sigma);
02212 
02213       // profile in angle
02214       deltaAngle = WCSimVertexGeometry::Instance()->GetAngle(idigit) - myConeEdge;
02215      
02216       // deweight hits outside cone
02217       if( deltaAngle<=0.0 ){
02218         deweight = 1.0;
02219       }
02220       else{
02221         deweight = 1.0/(1.0+(deltaAngle*deltaAngle)/(myConeEdgeSigma*myConeEdgeSigma));
02222       }
02223 
02224       // add to running total
02225       Swx += deweight*weight*delta;
02226       Sw  += deweight*weight;
02227     }
02228   }
02229 
02230   if( Sw>0.0 ){
02231     meanTime = Swx/Sw;
02232   }
02233 
02234   // return mean and rms
02235   // ===================
02236   vtxTime = meanTime;
02237 
02238   return;
02239 }
02240 
02241 void WCSimVertexFinder::TimePropertiesLnL(Double_t vtxTime, Double_t vtxParam, Double_t& vtxFOM)
02242 { 
02243   // nuisance parameters
02244   // ===================
02245   Double_t scatter = vtxParam; 
02246 
02247   // internal variables
02248   // ==================
02249   Double_t delta = 0.0;       // time residual of each hit
02250   Double_t sigma = 0.0;       // time resolution of each hit
02251 
02252   Double_t sigmaEarly = 1.5;  // width of early light
02253   Double_t sigmaLate  = 3.5;  // width of late light
02254   Double_t deltaShort = 5.0;  // decay time for scattered light [short]
02255   Double_t deltaLong  = 10.0; // decay time for scattered light [long]
02256 
02257   Double_t alpha = scatter;   // size of scattering (Gaussian+Exponential)
02258   Double_t beta  = 0.15;      // size of late light (second Gaussian)
02259   Double_t gamma = 0.66;      // coefficeient for double-exponential
02260 
02261   Double_t A = 0.0;           // normalisation of first Gaussian
02262   Double_t Ascat = 0.0;       // normalisation of second Gaussian
02263   Double_t Bshort = 0.0;      // normalisation of scattering (Gaussian+Exponential)
02264   Double_t Blong = 0.0;       // normalisation of scattering (Gaussian+Exponential)
02265 
02266   Double_t Preal = 0.0;       // probability of real hit
02267   Double_t P = 0.0;           // probability of hit
02268 
02269   Double_t chi2 = 0.0;        // log-likelihood: chi2 = -2.0*log(L)
02270   Double_t ndof = 0.0;        // total number of hits
02271   Double_t fom = 0.0;         // figure of merit
02272 
02273   // tuning parameters
02274   // =================
02275   Double_t fTimeFitNoiseRate = 0.10;  // hits/ns [0.40 for electrons, 0.02 for muons]
02276   
02277   // add noise to model
02278   // ==================
02279   Double_t nFilterDigits = WCSimVertexGeometry::Instance()->GetNFilterDigits();
02280   Double_t Pnoise = fTimeFitNoiseRate/nFilterDigits;
02281 
02282   // loop over digits
02283   // ================
02284   for( Int_t idigit=0; idigit<WCSimVertexGeometry::Instance()->GetNDigits(); idigit++ ){
02285 
02286     if( WCSimVertexGeometry::Instance()->IsFiltered(idigit) ){
02287       delta = WCSimVertexGeometry::Instance()->GetDelta(idigit) - vtxTime;
02288       sigma = WCSimVertexGeometry::Instance()->GetDeltaSigma(idigit);
02289 
02290       A      = 1.0 / ( 2.0*sigma*sqrt(0.5*TMath::Pi()) );
02291       Ascat  = 1.0 / ( (sigmaEarly+sigmaLate)*sqrt(0.5*TMath::Pi()) );
02292       Bshort = 1.0 / ( sigmaEarly*sqrt(0.5*TMath::Pi()) + (9.0/8.0)*deltaShort );
02293       Blong  = 1.0 / ( sigmaEarly*sqrt(0.5*TMath::Pi()) + (9.0/8.0)*deltaLong );
02294       
02295       if( delta<=0 ){
02296         Preal = (1.0-beta)*(1.0-alpha)*A*exp(-(delta*delta)/(2.0*sigma*sigma))
02297                 + beta*(1.0-alpha)*Ascat*exp(-(delta*delta)/(2.0*sigmaEarly*sigmaEarly))
02298                     + alpha*gamma*Bshort*exp(-(delta*delta)/(2.0*sigmaEarly*sigmaEarly))
02299                + alpha*(1.0-gamma)*Blong*exp(-(delta*delta)/(2.0*sigmaEarly*sigmaEarly));
02300         P = (1.0-Pnoise)*Preal + Pnoise;
02301       }
02302       else{
02303         Preal = (1.0-beta)*(1.0-alpha)*A*exp(-(delta*delta)/(2.0*sigma*sigma)) 
02304                 + beta*(1.0-alpha)*Ascat*exp(-(delta*delta)/(2.0*sigmaLate*sigmaLate)) 
02305                 + alpha*gamma*Bshort*( (delta/deltaShort)/pow((1.0+4.0*(delta/deltaShort)*(delta/deltaShort)),2.0) + exp(-delta/deltaShort) ) 
02306              + alpha*(1.0-gamma)*Blong*( (delta/deltaLong)/pow((1.0+4.0*(delta/deltaLong)*(delta/deltaLong)),2.0) + exp(-delta/deltaLong) );
02307         P = (1.0-Pnoise)*Preal + Pnoise;
02308       }
02309 
02310       chi2 += -2.0*log(P);
02311       ndof += 1.0;
02312     }
02313   }
02314 
02315   // calculate figure of merit
02316   // =========================   
02317   if( ndof>0.0 ){
02318     fom = fBaseFOM - 5.0*chi2/ndof;
02319   }
02320 
02321   // return figure of merit
02322   // ======================
02323   vtxFOM = fom;
02324 
02325   return;
02326 }
02327 
02328 void WCSimVertexFinder::ConePropertiesLnL(Double_t coneParam0, Double_t coneParam1, Double_t coneParam2, Double_t& coneAngle, Double_t& coneFOM)
02329 {  
02330   // nuisance parameters
02331   // ===================
02332   Double_t alpha  = coneParam0;
02333   Double_t alpha0 = coneParam1;
02334   Double_t beta   = coneParam2;
02335 
02336   // internal variables
02337   // ==================
02338   Double_t deltaAngle = 0.0;
02339   Double_t sigmaAngle = 7.0;
02340   Double_t deltaAngle0 = 42.0*alpha0;
02341   
02342   Double_t digitQ = 0.0;
02343   Double_t sigmaQmin = 1.0;
02344   Double_t sigmaQmax = 10.0;
02345   Double_t sigmaQ = 0.0;
02346 
02347   Double_t A = 0.0;
02348   
02349   Double_t PconeA = 0.0;
02350   Double_t PconeB = 0.0;
02351   Double_t Pmu = 0.0;
02352   Double_t Pel = 0.0;
02353 
02354   Double_t Pcharge = 0.0;
02355   Double_t Pangle = 0.0;
02356   Double_t P = 0.0;
02357 
02358   Double_t chi2 = 0.0;
02359   Double_t ndof = 0.0;
02360 
02361   Double_t angle = 46.0;
02362   Double_t fom = 0.0;
02363 
02364   // hard-coded parameters: 200 kton (100 kton)
02365   // ==========================================
02366   Double_t lambdaMuShort = 0.5; //  0.5;
02367   Double_t lambdaMuLong  = 5.0; // 15.0;
02368   Double_t alphaMu =       1.0; //  4.5;
02369 
02370   Double_t lambdaElShort = 1.0; //  2.5;
02371   Double_t lambdaElLong =  7.5; // 15.0;
02372   Double_t alphaEl =       6.0; //  3.5;
02373 
02374   // numerical integrals
02375   // ===================
02376   fSconeA = 21.9938;  
02377   fSconeB =  0.0000;
02378 
02379   // inside cone
02380   Int_t nbinsInside = 420;
02381   for( Int_t n=0; n<nbinsInside; n++ ){
02382     deltaAngle = -42.0 + (n+0.5)*(42.0/(double)nbinsInside);
02383     fSconeB += 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) )
02384                            *( 1.0/(1.0+(deltaAngle*deltaAngle)/(deltaAngle0*deltaAngle0)) )
02385                            *( 42.0/(double)nbinsInside );
02386   }
02387 
02388   // outside cone
02389   if( fIntegralsDone == 0 ){
02390     fSmu = 0.0;
02391     fSel = 0.0;
02392 
02393     Int_t nbinsOutside = 1380;
02394     for( Int_t n=0; n<nbinsOutside; n++ ){
02395       deltaAngle = 0.0 + (n+0.5)*(138.0/(double)nbinsOutside);
02396 
02397       fSmu += 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) )
02398                           *( 1.0/(1.0+alphaMu*(lambdaMuShort/lambdaMuLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaMuShort*lambdaMuShort)) 
02399                                             + alphaMu*(lambdaMuShort/lambdaMuLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaMuLong*lambdaMuLong)) )
02400                           *( 138.0/(double)nbinsOutside );
02401 
02402       fSel += 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) )
02403                           *( 1.0/(1.0+alphaEl*(lambdaElShort/lambdaElLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaElShort*lambdaElShort)) 
02404                                           + alphaEl*(lambdaElShort/lambdaElLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaElLong*lambdaElLong)) )
02405                           *( 138.0/(double)nbinsOutside );
02406     }
02407 
02408     std::cout << " --- calculating integrals: Smu=" << fSmu << " Sel=" << fSel << std::endl;
02409 
02410     fIntegralsDone = 1;
02411   }
02412 
02413 
02414   // loop over digits
02415   // ================
02416   for( Int_t idigit=0; idigit<WCSimVertexGeometry::Instance()->GetNDigits(); idigit++ ){
02417 
02418     if( WCSimVertexGeometry::Instance()->IsFiltered(idigit) ){
02419       digitQ = WCSimVertexGeometry::Instance()->GetDigitQ(idigit);
02420       deltaAngle = WCSimVertexGeometry::Instance()->GetAngle(idigit) - 42.0;
02421 
02422       // pulse height distribution
02423       // =========================
02424       if( deltaAngle<=0 ){
02425         sigmaQ = sigmaQmax;
02426       }
02427       else{
02428         sigmaQ = sigmaQmin + (sigmaQmax-sigmaQmin)/(1.0+(deltaAngle*deltaAngle)/(sigmaAngle*sigmaAngle));
02429       }
02430 
02431       A = 1.0/(log(2.0)+0.5*TMath::Pi()*sigmaQ);
02432 
02433       if( digitQ<1.0 ){
02434         Pcharge = 2.0*A*digitQ/(1.0+digitQ*digitQ);
02435       }
02436       else{
02437         Pcharge = A/(1.0+(digitQ-1.0)*(digitQ-1.0)/(sigmaQ*sigmaQ));
02438       }
02439 
02440       // angular distribution
02441       // ====================
02442       A = 1.0/( alpha*fSconeA + (1.0-alpha)*fSconeB 
02443                + beta*fSmu + (1.0-beta)*fSel ); // numerical integrals 
02444 
02445       if( deltaAngle<=0 ){
02446 
02447         // pdfs inside cone:
02448         PconeA = 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) );
02449         PconeB = 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) )
02450                           *( 1.0/(1.0+(deltaAngle*deltaAngle)/(deltaAngle0*deltaAngle0)) );
02451 
02452         Pangle = A*( alpha*PconeA+(1.0-alpha)*PconeB );
02453       }         
02454       else{
02455 
02456         // pdfs outside cone
02457         Pmu = 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) )
02458                        *( 1.0/(1.0+alphaMu*(lambdaMuShort/lambdaMuLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaMuShort*lambdaMuShort)) 
02459                                           + alphaMu*(lambdaMuShort/lambdaMuLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaMuLong*lambdaMuLong)) );
02460 
02461         Pel = 1.4944765*sin( (42.0+deltaAngle)*(TMath::Pi()/180.0) )
02462                        *( 1.0/(1.0+alphaEl*(lambdaElShort/lambdaElLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaElShort*lambdaElShort)) 
02463                                           + alphaEl*(lambdaElShort/lambdaElLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaElLong*lambdaElLong)) );
02464 
02465         Pangle = A*( beta*Pmu+(1.0-beta)*Pel );
02466       }
02467 
02468       // overall probability
02469       // ===================
02470       P = Pcharge*Pangle;
02471       
02472       chi2 += -2.0*log(P);
02473       ndof += 1.0;
02474     }
02475   }
02476 
02477   // calculate figure of merit
02478   // =========================   
02479   if( ndof>0.0 ){
02480     fom = fBaseFOM - 5.0*chi2/ndof;
02481     angle = beta*43.0 + (1.0-beta)*49.0;
02482   }
02483 
02484   // return figure of merit
02485   // ======================
02486   coneAngle = angle;
02487   coneFOM = fom;
02488 
02489   return;
02490 }
02491 
02492 void WCSimVertexFinder::TimePropertiesFoM(Double_t vtxTime, Double_t& vtxFOM)
02493 { 
02494   // calculate figure of merit
02495   // =========================   
02496   Double_t Swx = 0.0;
02497   Double_t Sw = 0.0; 
02498 
02499   Double_t delta = 0.0;
02500   Double_t sigma = 0.0;
02501 
02502   Double_t weight = 0.0;
02503   Double_t deweight = 0.0;
02504   Double_t inweight = 0.0;
02505   Double_t outweight = 0.0;
02506   Double_t timeweight = 0.0;
02507 
02508   Double_t deltaAngle = 0.0;
02509   Double_t deltaEdge = 0;
02510 
02511   Double_t alpha = 0.0;
02512   Double_t penalty = 0.011; 
02513  
02514   // [ Weight = exp(-(delta*delta)/(2.0*sigma*sigma)) - Penalty,
02515   //   This drops below zero at Penalty = sqrt( 2.0*log(Weight) )
02516   //   To drop below zero at 2 sigma,   set Penalty = exp(-4/2)    = 0.135
02517   //   To drop below zero at 2.5 sigma, set Penalty = exp(-6.25/2) = 0.044
02518   //   To drop below zero at 3 sigma,   set Penalty = exp(-9/2)    = 0.011 ]
02519  
02520   Double_t sigmaRes = 2.0;         // time resolution
02521   Double_t sigmaLow = 2.0;         // time resolution inside cone
02522   Double_t sigmaHigh = 2.0;        // time resolution outside cone
02523 
02524                                                           // muons electrons
02525   Double_t myConeAngle = 42.0;     // cherenkov cone      // 42.0    42.0
02526   Double_t myConeAngleSigma = 2.0; //   [degrees]         //  2.0     2.0
02527                                                           
02528   Double_t myConeEdge = 42.0;      // overall window      // 42.0    48.0
02529   Double_t myConeEdgeSigma = 2.0;  //   [degrees]         //  2.0     8.0
02530 
02531   Double_t fom = 0.0;
02532 
02533   for( Int_t idigit=0; idigit<WCSimVertexGeometry::Instance()->GetNDigits(); idigit++ ){
02534 
02535     if( WCSimVertexGeometry::Instance()->IsFiltered(idigit) ){
02536       delta = WCSimVertexGeometry::Instance()->GetDelta(idigit) - vtxTime;
02537       sigma = WCSimVertexGeometry::Instance()->GetDeltaSigma(idigit);
02538       weight = 1.0/(sigma*sigma);
02539 
02540       deltaAngle = WCSimVertexGeometry::Instance()->GetAngle(idigit) - myConeAngle;
02541       deltaEdge  = WCSimVertexGeometry::Instance()->GetAngle(idigit) - myConeEdge;
02542 
02543       // time weights inside and outside cone
02544       inweight  = exp(-(delta*delta)/(2.0*sigmaRes*sigmaRes))-penalty;
02545 
02546       if( delta<=0 ){ 
02547         outweight = (2.0*sigmaRes/(sigmaLow+sigmaHigh))*(exp(-(delta*delta)/(2.0*sigmaLow*sigmaLow))-penalty);
02548       }
02549       else{
02550         outweight = (2.0*sigmaRes/(sigmaLow+sigmaHigh))*(exp(-(delta*delta)/(2.0*sigmaHigh*sigmaHigh))-penalty);
02551       }
02552       
02553       if( deltaAngle<=0 ){
02554         alpha = 1.0;
02555       }
02556       else{
02557         alpha = 1.0/(1.0+(deltaAngle*deltaAngle)/(myConeAngleSigma*myConeAngleSigma));
02558       }
02559 
02560       // deweight hits outside cone
02561       if( deltaEdge<=0.0 ){ 
02562         deweight = 1.0;
02563       } 
02564       else{
02565         deweight = 1.0/(1.0+(deltaEdge*deltaEdge)/(myConeEdgeSigma*myConeEdgeSigma));
02566       }
02567 
02568       // overall time weight
02569       timeweight = inweight*alpha + outweight*(1.0-alpha);
02570 
02571       // add to running total
02572       Swx += deweight*weight*timeweight;
02573       Sw  += deweight*weight;
02574     }
02575   }
02576 
02577   if( Sw>0.0 ){
02578     fom = fBaseFOM*Swx/Sw;
02579   }
02580 
02581   // return figure of merit
02582   // ======================
02583   vtxFOM = fom;
02584 
02585   return;
02586 }
02587 
02588 void WCSimVertexFinder::ConePropertiesFoM(Double_t& coneFOM)
02589 {  
02590   // calculate figure of merit
02591   // =========================
02592   Double_t coneEdge = 42.0;     // nominal cone angle
02593   Double_t coneEdgeLow = 21.0;  // cone edge (low side)      
02594   Double_t coneEdgeHigh = 3.0;  // cone edge (high side)   [muons: 3.0, electrons: 7.0]
02595 
02596   Double_t deltaAngle = 0.0;
02597   Double_t digitCharge = 0.0;
02598  
02599   Double_t coneCharge = 0.0;
02600   Double_t allCharge = 0.0;
02601 
02602   Double_t fom = 0.0;
02603 
02604   for( Int_t idigit=0; idigit<WCSimVertexGeometry::Instance()->GetNDigits(); idigit++ ){
02605 
02606     if( WCSimVertexGeometry::Instance()->IsFiltered(idigit) ){
02607       deltaAngle = WCSimVertexGeometry::Instance()->GetAngle(idigit) - coneEdge;
02608       digitCharge = WCSimVertexGeometry::Instance()->GetDigitQ(idigit);
02609 
02610       if( deltaAngle<=0.0 ){ 
02611         coneCharge += digitCharge*( 0.75 + 0.25/( 1.0 + (deltaAngle*deltaAngle)/(coneEdgeLow*coneEdgeLow) ) );
02612       }
02613       else{ 
02614         coneCharge += digitCharge*( 0.00 + 1.00/( 1.0 + (deltaAngle*deltaAngle)/(coneEdgeHigh*coneEdgeHigh) ) );
02615       }
02616 
02617       allCharge += digitCharge;
02618     }
02619   }
02620 
02621   if( allCharge>0.0 ){
02622     fom = fBaseFOM*coneCharge/allCharge;
02623   }
02624 
02625   // return figure of merit
02626   // ======================
02627   coneFOM = fom;
02628 
02629   return;
02630 }
02631 
02632 void WCSimVertexFinder::PenaltyChi2(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t& chi2)
02633 {  
02634   // default penalty
02635   // ===============
02636   Double_t deltaChi2 = 0;
02637 
02638   // check vertex position
02639   // =====================
02640   Bool_t inDetector = WCSimGeometry::Instance()->InsideDetector(vtxX,vtxY,vtxZ);
02641   if( inDetector==0 ){
02642     
02643     Double_t deltaR = WCSimGeometry::Instance()->DistanceToEdge(vtxX,vtxY,vtxZ);  // [cm]
02644     Double_t deltaRbase = 500.0; // [cm]
02645  
02646     deltaChi2 += -fBaseFOM*(deltaR/deltaRbase)*(deltaR/deltaRbase);
02647   }
02648 
02649   // return penalty
02650   // ==============
02651   chi2 = deltaChi2;
02652 
02653   return;
02654 }
02655 
02656 void WCSimVertexFinder::FixPositionChi2(Double_t vx, Double_t vy, Double_t vz, Double_t& chi2)
02657 {
02658   // default penalty
02659   // ===============
02660   Double_t deltaChi2 = 0;
02661 
02662   // get true event
02663   // ==============
02664   WCSimTrueEvent* myTrueEvent = (WCSimTrueEvent*)(WCSimInterface::TrueEvent());
02665 
02666   if( myTrueEvent->GetNTracks()>0 ){
02667     Double_t trueVtxX = myTrueEvent->GetVtxX();
02668     Double_t trueVtxY = myTrueEvent->GetVtxY();
02669     Double_t trueVtxZ = myTrueEvent->GetVtxZ();
02670   
02671     Double_t deltaRsq = (vx-trueVtxX)*(vx-trueVtxX)
02672                        +(vy-trueVtxY)*(vy-trueVtxY)
02673                        +(vz-trueVtxZ)*(vz-trueVtxZ);
02674     Double_t deltaRsqBase = 0.25*0.25; // [cm]
02675 
02676     deltaChi2 += -10.0*fBaseFOM*(1.0-1.0/(1.0+(deltaRsq/deltaRsqBase)));
02677   }
02678 
02679   // return penalty
02680   // ==============
02681   chi2 = deltaChi2;
02682 
02683   return;
02684 }
02685 
02686 void WCSimVertexFinder::FixDirectionChi2(Double_t px, Double_t py, Double_t pz, Double_t& chi2)
02687 {
02688   // default penalty
02689   // ===============
02690   Double_t deltaChi2 = 0;
02691 
02692   // get true event
02693   // ==============
02694   WCSimTrueEvent* myTrueEvent = (WCSimTrueEvent*)(WCSimInterface::TrueEvent());
02695 
02696   if( myTrueEvent->GetNTracks()>0 ){
02697     Double_t trueDirX = myTrueEvent->GetDirX();
02698     Double_t trueDirY = myTrueEvent->GetDirY();
02699     Double_t trueDirZ = myTrueEvent->GetDirZ();
02700   
02701     Double_t deltaTheta = (180.0/TMath::Pi())*(TMath::ACos(px*trueDirX+py*trueDirY+pz*trueDirZ));
02702     Double_t deltaThetaBase = 0.5; // [degrees]
02703 
02704     deltaChi2 += -10.0*fBaseFOM*(1.0-1.0/(1.0+(deltaTheta/deltaThetaBase)*(deltaTheta/deltaThetaBase)));
02705   }
02706 
02707   // return penalty
02708   // ==============
02709   chi2 = deltaChi2;
02710 
02711   return;
02712 }
02713