WCSimVertexGeometry.cc

Go to the documentation of this file.
00001 #include "WCSimVertexGeometry.hh"
00002 
00003 #include "WCSimRecoEvent.hh"
00004 #include "WCSimRecoVertex.hh"
00005 #include "WCSimRecoDigit.hh"
00006 
00007 #include "WCSimGeometry.hh"
00008 #include "WCSimParameters.hh"
00009 
00010 #include "TMath.h"
00011 #include "TRandom.h"
00012 
00013 #include <vector>
00014 
00015 #include <cmath>
00016 #include <iostream>
00017 #include <cassert>
00018 
00019 ClassImp(WCSimVertexGeometry)
00020 
00021 static WCSimVertexGeometry* fgVertexGeometry = 0;
00022 
00023 WCSimVertexGeometry* WCSimVertexGeometry::Instance()
00024 {
00025   if( !fgVertexGeometry ){
00026     fgVertexGeometry = new WCSimVertexGeometry();
00027   }
00028 
00029   if( !fgVertexGeometry ){
00030     assert(fgVertexGeometry);
00031   }
00032 
00033   if( fgVertexGeometry ){
00034 
00035   }
00036 
00037   return fgVertexGeometry;
00038 }
00039 
00040 WCSimVertexGeometry::WCSimVertexGeometry()
00041 {
00042   if( WCSimGeometry::TouchGeometry() ){
00043     fPMTs = WCSimGeometry::Instance()->GetNumPMTs();
00044   }
00045   else{
00046     fPMTs = 100000; // maximum number of digits
00047   }
00048   
00049   fRunNum = -1;
00050   fEventNum = -1;
00051   fTriggerNum = -1;
00052   
00053   fNDigits = 0;
00054   fNFilterDigits = 0;
00055 
00056   fThisDigit = 0;
00057   fLastEntry = 0; 
00058   fCounter = 0;
00059 
00060   fMeanQ = 0.0;
00061   fTotalQ = 0.0;
00062   fMeanFilteredQ = 0.0;
00063   fTotalFilteredQ = 0.0;
00064 
00065   fMinTime = 0.0;
00066   fMaxTime = 0.0;
00067   
00068   fVtxX1 = 0.0;
00069   fVtxY1 = 0.0;
00070   fVtxZ1 = 0.0;
00071   fVtxTime1 = 0.0;
00072 
00073   fVtxX2 = 0.0;
00074   fVtxY2 = 0.0;
00075   fVtxZ2 = 0.0;
00076   fVtxTime2 = 0.0;
00077 
00078   fIsFiltered = new Bool_t[fPMTs];
00079 
00080   fDigitX = new Double_t[fPMTs];
00081   fDigitY = new Double_t[fPMTs];
00082   fDigitZ = new Double_t[fPMTs];
00083   fDigitT = new Double_t[fPMTs];
00084   fDigitQ = new Double_t[fPMTs];
00085 
00086   fConeAngle = new Double_t[fPMTs];
00087   fZenith = new Double_t[fPMTs];
00088   fAzimuth = new Double_t[fPMTs];
00089   fSolidAngle = new Double_t[fPMTs];
00090 
00091   fDistPoint = new Double_t[fPMTs];
00092   fDistTrack = new Double_t[fPMTs];
00093   fDistPhoton = new Double_t[fPMTs];
00094   fDistScatter = new Double_t[fPMTs];
00095 
00096   fDeltaTime = new Double_t[fPMTs];
00097   fDeltaSigma = new Double_t[fPMTs];
00098 
00099   fDeltaAngle = new Double_t[fPMTs];
00100   fDeltaPoint = new Double_t[fPMTs];
00101   fDeltaTrack = new Double_t[fPMTs];
00102   fDeltaPhoton = new Double_t[fPMTs];
00103   fDeltaScatter = new Double_t[fPMTs];
00104 
00105   fPointPath = new Double_t[fPMTs];
00106   fExtendedPath = new Double_t[fPMTs];
00107 
00108   fPointResidual = new Double_t[fPMTs];
00109   fExtendedResidual = new Double_t[fPMTs];
00110 
00111   fDelta = new Double_t[fPMTs];
00112 
00113   for( Int_t n=0; n<fPMTs; n++ ){
00114     fIsFiltered[n] = 0.0;
00115 
00116     fDigitX[n] = 0.0;
00117     fDigitY[n] = 0.0;
00118     fDigitZ[n] = 0.0;
00119     fDigitT[n] = 0.0;
00120     fDigitQ[n] = 0.0;
00121 
00122     fConeAngle[n] = 0.0;
00123     fZenith[n] = 0.0;
00124     fAzimuth[n] = 0.0;
00125     fSolidAngle[n] = 0.0;
00126 
00127     fDistPoint[n] = 0.0;
00128     fDistTrack[n] = 0.0;
00129     fDistPhoton[n] = 0.0; 
00130     fDistScatter[n] = 0.0; 
00131 
00132     fDeltaTime[n] = 0.0;
00133     fDeltaSigma[n] = 0.0;
00134 
00135     fDeltaAngle[n] = 0.0;    
00136     fDeltaPoint[n] = 0.0;
00137     fDeltaTrack[n] = 0.0;
00138     fDeltaPhoton[n] = 0.0;
00139     fDeltaScatter[n] = 0.0;
00140 
00141     fPointPath[n] = 0.0;
00142     fExtendedPath[n] = 0.0;
00143     fPointResidual[n] = 0.0;
00144     fExtendedResidual[n] = 0.0;
00145 
00146     fDelta[n] = 0.0;
00147   }
00148 }
00149 
00150 WCSimVertexGeometry::~WCSimVertexGeometry()
00151 {
00152   if( fIsFiltered ) delete [] fIsFiltered;
00153 
00154   if( fDigitX ) delete [] fDigitX;
00155   if( fDigitY ) delete [] fDigitY;
00156   if( fDigitZ ) delete [] fDigitZ;
00157   if( fDigitT ) delete [] fDigitT;
00158   if( fDigitQ ) delete [] fDigitQ; 
00159 
00160   if( fConeAngle ) delete [] fConeAngle;
00161   if( fZenith ) delete [] fZenith;
00162   if( fAzimuth ) delete [] fAzimuth;
00163   if( fSolidAngle ) delete [] fSolidAngle;
00164 
00165   if( fDistPoint )  delete [] fDistPoint;
00166   if( fDistTrack )  delete [] fDistTrack;
00167   if( fDistPhoton ) delete [] fDistPhoton;   
00168   if( fDistScatter ) delete [] fDistScatter;
00169 
00170   if( fDeltaTime )   delete [] fDeltaTime;
00171   if( fDeltaSigma ) delete [] fDeltaSigma;
00172 
00173   if( fDeltaAngle )  delete [] fDeltaAngle;  
00174   if( fDeltaPoint )  delete [] fDeltaPoint;
00175   if( fDeltaTrack )  delete [] fDeltaTrack;
00176   if( fDeltaPhoton ) delete [] fDeltaPhoton;
00177   if( fDeltaScatter ) delete [] fDeltaScatter;
00178 
00179   if( fPointPath ) delete [] fPointPath;
00180   if( fExtendedPath ) delete [] fExtendedPath;
00181   if( fPointResidual )  delete [] fPointResidual;
00182   if( fExtendedResidual ) delete [] fExtendedResidual;
00183 
00184   if( fDelta ) delete [] fDelta;  
00185 }
00186 
00187 void WCSimVertexGeometry::LoadEvent(WCSimRecoEvent* myEvent)
00188 {
00189   // resetting arrays
00190   // ================
00191   if( myEvent==0 ) return;
00192 
00193   Int_t thisRunNum = myEvent->GetRun();
00194   Int_t thisEventNum = myEvent->GetEvent();
00195   Int_t thisTriggerNum = myEvent->GetTrigger();
00196 
00197   if( fRunNum==thisRunNum
00198    && fEventNum==thisEventNum
00199    && fTriggerNum==thisTriggerNum ){
00200     return;
00201   }
00202 
00203   fRunNum = thisRunNum;
00204   fEventNum = thisEventNum;
00205   fTriggerNum = thisTriggerNum;
00206 
00207   fNDigits = 0;
00208   fNFilterDigits = 0;
00209   
00210   fThisDigit = 0;
00211   fLastEntry = 0;
00212   fCounter = 0;
00213 
00214   // clear vertices
00215   // ==============
00216   for(UInt_t i=0; i<vVertexList.size(); i++ ){
00217     delete (WCSimRecoVertex*)(vVertexList.at(i));
00218   }
00219   vVertexList.clear();
00220 
00221   // clear seed vertices
00222   // ===================
00223   vSeedVtxX.clear();
00224   vSeedVtxY.clear();
00225   vSeedVtxZ.clear();
00226   vSeedVtxTime.clear();
00227   vSeedDigitList.clear();
00228 
00229   // get lists of digits
00230   // ==================
00231   std::vector<WCSimRecoDigit*>* myDigitList = (std::vector<WCSimRecoDigit*>*)(myEvent->GetDigitList());
00232   fNDigits = myDigitList->size();
00233 
00234   std::vector<WCSimRecoDigit*>* myFilterDigitList = (std::vector<WCSimRecoDigit*>*)(myEvent->GetFilterDigitList());
00235   fNFilterDigits = myFilterDigitList->size();
00236 
00237   // sanity checks
00238   // =============
00239   if( fNDigits<=0 ){
00240     std::cout << " *** WCSimVertexGeometry::LoadEvent(...) *** " << std::endl;
00241     std::cout << "   <warning> event has no digits! " << std::endl;
00242     return;
00243   }
00244 
00245   if( fNDigits>fPMTs ){
00246     std::cout << " *** WCSimVertexGeometry::LoadEvent(...) *** " << std::endl;
00247     std::cout << "   <warning> event will be truncated! " << std::endl;
00248     fNDigits = fPMTs;
00249   }
00250 
00251   if( fNFilterDigits>fNDigits ){
00252     fNFilterDigits = fNDigits;
00253   }
00254 
00255   // load digits
00256   // ===========
00257   fTotalQ = 0.0;
00258   fMeanQ = 0.0;
00259 
00260   fTotalFilteredQ = 0.0;
00261   fMeanFilteredQ = 0.0;
00262 
00263   fMinTime = -999999.9;
00264   fMaxTime = -999999.9;
00265 
00266   Double_t Swx = 0.0;
00267   Double_t Sw = 0.0;
00268 
00269   Double_t SFwx = 0.0;
00270   Double_t SFw = 0.0;
00271   
00272   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00273     WCSimRecoDigit* recoDigit = (WCSimRecoDigit*)(myDigitList->at(idigit));    
00274 
00275     fDigitX[idigit] = recoDigit->GetX();
00276     fDigitY[idigit] = recoDigit->GetY();
00277     fDigitZ[idigit] = recoDigit->GetZ();
00278     fDigitT[idigit] = recoDigit->GetTime();
00279     fDigitQ[idigit] = recoDigit->GetQPEs();
00280 
00281     fIsFiltered[idigit] = recoDigit->IsFiltered();
00282 
00283     fConeAngle[idigit] = 0.0;
00284     fZenith[idigit] = 0.0;
00285     fAzimuth[idigit] = 0.0;
00286     fSolidAngle[idigit] = 0.0;
00287 
00288     fDistPoint[idigit] = 0.0;
00289     fDistTrack[idigit] = 0.0;
00290     fDistPhoton[idigit] = 0.0;
00291     fDistScatter[idigit] = 0.0;
00292 
00293     fDeltaTime[idigit] = 0.0;
00294     fDeltaSigma[idigit] = 0.0;
00295     
00296     fDeltaAngle[idigit] = 0.0;
00297     fDeltaPoint[idigit] = 0.0;
00298     fDeltaTrack[idigit] = 0.0;
00299     fDeltaPhoton[idigit] = 0.0;
00300     fDeltaScatter[idigit] = 0.0;
00301 
00302     fPointPath[idigit] = 0.0;
00303     fExtendedPath[idigit] = 0.0;
00304     fPointResidual[idigit] = 0.0;
00305     fExtendedResidual[idigit] = 0.0;
00306 
00307     fDelta[idigit] = 0.0;
00308 
00309     Swx += recoDigit->GetQPEs();
00310     Sw += 1.0;
00311 
00312     if( fMinTime<0 
00313      || recoDigit->GetTime()<fMinTime ){
00314       fMinTime = recoDigit->GetTime();
00315     }
00316 
00317     if( fMaxTime<0 
00318      || recoDigit->GetTime()>fMaxTime ){
00319       fMaxTime = recoDigit->GetTime();
00320     }
00321 
00322     if( recoDigit->IsFiltered() ){
00323       SFwx += recoDigit->GetQPEs();
00324       SFw += 1.0;
00325     }
00326   }
00327 
00328   if( Sw>0.0 ){
00329     fTotalQ = Swx;
00330     fMeanQ = Swx/Sw;
00331   }
00332 
00333   if( SFw>0.0 ){
00334     fTotalFilteredQ = SFwx;
00335     fMeanFilteredQ = SFwx/SFw;
00336   }
00337 
00338   if( fMinTime<0 ) fMinTime = 0.0;
00339   if( fMaxTime<0 ) fMaxTime = 0.0;
00340 
00341   return;
00342 }
00343   
00344 WCSimRecoVertex* WCSimVertexGeometry::CalcSimpleVertex(WCSimRecoEvent* myEvent)
00345 {
00346   // load event
00347   // ==========
00348   this->LoadEvent(myEvent);
00349 
00350   // calculate simple vertex
00351   // =======================
00352   return this->CalcSimpleVertex(myEvent);
00353 }
00354 
00355 WCSimRecoVertex* WCSimVertexGeometry::CalcSimpleVertex()
00356 {
00357   // simple vertex
00358   // =============
00359   Double_t vtxX = 0.0;
00360   Double_t vtxY = 0.0;
00361   Double_t vtxZ = 0.0;
00362   Double_t vtxTime = 0.0;
00363 
00364   Int_t itr = 1;
00365   Bool_t pass = 1; 
00366   Double_t fom = 1.0;
00367 
00368   // create new vertex
00369   // =================
00370   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00371   vVertexList.push_back(newVertex);
00372 
00373   // calculate vertex
00374   // ================
00375   this->CalcSimpleVertex(vtxX,vtxY,vtxZ,vtxTime);
00376 
00377   // set vertex
00378   // ==========
00379   newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00380   newVertex->SetFOM(fom,itr,pass);  
00381 
00382   // set status
00383   // ==========
00384   newVertex->SetStatus(WCSimRecoVertex::kOK);
00385     
00386   // return vertex
00387   // =============
00388   return newVertex;
00389 }
00390 
00391 void WCSimVertexGeometry::CalcSimpleVertex(Double_t& vtxX, Double_t& vtxY, Double_t& vtxZ, Double_t& vtxTime)
00392 {
00393   // simple vertex
00394   // =============
00395   // just calculate average position of digits
00396 
00397   // default vertex
00398   // ==============
00399   vtxX = 0.0;
00400   vtxY = 0.0;
00401   vtxZ = 0.0;
00402   vtxTime = 0.0;  
00403 
00404   // loop over digits
00405   // ================
00406   Double_t Swx = 0.0;
00407   Double_t Swy = 0.0;
00408   Double_t Swz = 0.0;
00409   Double_t Swt = 0.0;
00410   Double_t Sw = 0.0;
00411 
00412   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00413     if( fIsFiltered[idigit] ){
00414       Swx += fDigitQ[idigit]*fDigitX[idigit];
00415       Swy += fDigitQ[idigit]*fDigitY[idigit];
00416       Swz += fDigitQ[idigit]*fDigitZ[idigit];
00417       Swt += fDigitQ[idigit]*fDigitT[idigit];
00418       Sw  += fDigitQ[idigit];
00419     }
00420   }
00421 
00422   // average position
00423   // ================
00424   if( Sw>0.0 ){
00425     vtxX = Swx/Sw;
00426     vtxY = Swy/Sw;
00427     vtxZ = Swz/Sw;
00428     vtxTime = Swt/Sw;
00429   }   
00430 
00431   return;
00432 }
00433   
00434 WCSimRecoVertex* WCSimVertexGeometry::CalcSimpleDirection(WCSimRecoEvent* myEvent, WCSimRecoVertex* myVertex)
00435 {
00436   // load event
00437   // ==========
00438   this->LoadEvent(myEvent);
00439 
00440   // calculate simple direction
00441   // ==========================
00442   return this->CalcSimpleDirection(myVertex);
00443 }
00444 
00445 WCSimRecoVertex* WCSimVertexGeometry::CalcSimpleDirection(WCSimRecoVertex* myVertex)
00446 {
00447   // load vertex
00448   // ===========
00449   Double_t vtxX = myVertex->GetX();
00450   Double_t vtxY = myVertex->GetY();
00451   Double_t vtxZ = myVertex->GetZ();
00452   Double_t vtxTime = myVertex->GetTime();
00453     
00454   // current status
00455   // ==============
00456   Int_t status = myVertex->GetStatus();
00457 
00458   // create new vertex
00459   // =================
00460   WCSimRecoVertex* newVertex = new WCSimRecoVertex();
00461   vVertexList.push_back(newVertex);
00462 
00463   // loop over digits
00464   // ================
00465   Double_t Swx = 0.0;
00466   Double_t Swy = 0.0;
00467   Double_t Swz = 0.0;
00468   Double_t Sw = 0.0;
00469 
00470   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00471     if( fIsFiltered[idigit] ){
00472       Double_t q = fDigitQ[idigit];
00473 
00474       Double_t dx = fDigitX[idigit] - vtxX;
00475       Double_t dy = fDigitY[idigit] - vtxY;
00476       Double_t dz = fDigitZ[idigit] - vtxZ;
00477       Double_t ds = sqrt(dx*dx+dy*dy+dz*dz);
00478 
00479       Double_t px = dx/ds;
00480       Double_t py = dy/ds;
00481       Double_t pz = dz/ds;
00482 
00483       Swx += q*px;
00484       Swy += q*py;
00485       Swz += q*pz;
00486       Sw  += q;
00487     }
00488   }
00489 
00490   // average direction
00491   // =================
00492   Double_t dirX = 0.0;
00493   Double_t dirY = 0.0;
00494   Double_t dirZ = 0.0;
00495   
00496   Int_t itr = 0;
00497   Bool_t pass = 0; 
00498   Double_t fom = 0.0;
00499 
00500   if( Sw>0.0 ){
00501     Double_t qx = Swx/Sw;
00502     Double_t qy = Swy/Sw;
00503     Double_t qz = Swz/Sw;
00504     Double_t qs = sqrt(qx*qx+qy*qy+qz*qz);
00505 
00506     dirX = qx/qs;
00507     dirY = qy/qs;
00508     dirZ = qz/qs;
00509 
00510     fom = 1.0;
00511     itr = 1;
00512     pass = 1; 
00513   }
00514 
00515   // set vertex and direction
00516   // ========================
00517   if( pass ){
00518     newVertex->SetVertex(vtxX,vtxY,vtxZ,vtxTime);
00519     newVertex->SetDirection(dirX,dirY,dirZ);
00520     newVertex->SetFOM(fom,itr,pass);
00521   }
00522 
00523   // set status
00524   // ==========
00525   if( !pass ) status |= WCSimRecoVertex::kFailSimpleDirection;
00526   newVertex->SetStatus(status);
00527 
00528   // return vertex
00529   // =============
00530   return newVertex;
00531 }
00532 
00533 void WCSimVertexGeometry::CalcResiduals(WCSimRecoEvent* myEvent, WCSimRecoVertex* myVertex)
00534 {
00535   // load event
00536   // ==========
00537   this->LoadEvent(myEvent);
00538 
00539   // calculate residuals
00540   // ===================
00541   return this->CalcResiduals(myVertex);
00542 }
00543 
00544 void WCSimVertexGeometry::CalcResiduals(WCSimRecoVertex* myVertex)
00545 {
00546   // sanity check
00547   // ============
00548   if( myVertex==0 ) return;
00549 
00550   // get vertex information
00551   // ======================
00552   Double_t vtxX = myVertex->GetX();
00553   Double_t vtxY = myVertex->GetY();
00554   Double_t vtxZ = myVertex->GetZ();
00555   Double_t vtxTime = myVertex->GetTime();
00556 
00557   Double_t dirX = myVertex->GetDirX();
00558   Double_t dirY = myVertex->GetDirY();
00559   Double_t dirZ = myVertex->GetDirZ();
00560 
00561   // calculate residuals
00562   // ===================
00563   return this->CalcResiduals( vtxX, vtxY, vtxZ, vtxTime,
00564                               dirX, dirY, dirZ );
00565 }
00566 
00567 void WCSimVertexGeometry::CalcPointResiduals(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t vtxTime, Double_t dirX, Double_t dirY, Double_t dirZ)
00568 {
00569   this->CalcResiduals( vtxX, vtxY, vtxZ, vtxTime,
00570                        dirX, dirY, dirZ );
00571 
00572   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00573     fDelta[idigit] = fPointResidual[idigit];
00574   }
00575 
00576   return;
00577 }
00578 
00579 void WCSimVertexGeometry::CalcExtendedResiduals(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t vtxTime, Double_t dirX, Double_t dirY, Double_t dirZ )
00580 {
00581   this->CalcResiduals( vtxX, vtxY, vtxZ, vtxTime,
00582                        dirX, dirY, dirZ );
00583 
00584   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00585     fDelta[idigit] = fExtendedResidual[idigit];
00586   }
00587 
00588   return;
00589 }
00590 
00591 void WCSimVertexGeometry::CalcResiduals(Double_t vtxX, Double_t vtxY, Double_t vtxZ, Double_t vtxTime, Double_t dirX, Double_t dirY, Double_t dirZ )
00592 {
00593   // reset arrays
00594   // ============
00595   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00596     fConeAngle[idigit] = 0.0;
00597     fZenith[idigit] = 0.0;
00598     fAzimuth[idigit] = 0.0;
00599     fSolidAngle[idigit] = 0.0;
00600     fDistPoint[idigit] = 0.0;
00601     fDistTrack[idigit] = 0.0;
00602     fDistPhoton[idigit] = 0.0;  
00603     fDistScatter[idigit] = 0.0; 
00604     fDeltaTime[idigit] = 0.0; 
00605     fDeltaSigma[idigit] = 0.0;
00606     fDeltaAngle[idigit] = 0.0;
00607     fDeltaPoint[idigit] = 0.0;
00608     fDeltaTrack[idigit] = 0.0;
00609     fDeltaPhoton[idigit] = 0.0;
00610     fDeltaScatter[idigit] = 0.0;
00611     fPointPath[idigit] = 0.0;
00612     fExtendedPath[idigit] = 0.0;
00613     fPointResidual[idigit] = 0.0;
00614     fExtendedResidual[idigit] = 0.0;
00615     fDelta[idigit] = 0.0;
00616   }
00617 
00618   // cone angle
00619   // ==========
00620   Double_t thetadeg = WCSimParameters::CherenkovAngle(); // degrees
00621   Double_t theta = thetadeg*(TMath::Pi()/180.0); // degrees->radians
00622 
00623   // loop over digits
00624   // ================
00625   for( Int_t idigit=0; idigit<fNDigits; idigit++ ){
00626     Double_t dx = fDigitX[idigit]-vtxX;
00627     Double_t dy = fDigitY[idigit]-vtxY;
00628     Double_t dz = fDigitZ[idigit]-vtxZ;
00629     Double_t ds = sqrt(dx*dx+dy*dy+dz*dz);
00630 
00631     Double_t px = dx/ds;
00632     Double_t py = dy/ds;
00633     Double_t pz = dz/ds;
00634 
00635     Double_t cosphi = 1.0;
00636     Double_t sinphi = 1.0;
00637     Double_t phi = 0.0; 
00638     Double_t phideg = 0.0;
00639 
00640     Double_t ax = 0.0;
00641     Double_t ay = 0.0;
00642     Double_t az = 0.0;
00643     Double_t azideg = 0.0;
00644 
00645     // calculate angles if direction is known
00646     if( dirX*dirX + dirY*dirY + dirZ*dirZ>0.0 ){
00647 
00648       // zenith angle
00649       cosphi = px*dirX+py*dirY+pz*dirZ;
00650       phi = acos(cosphi); // radians
00651       phideg = phi/(TMath::Pi()/180.0); // radians->degrees
00652       sinphi = sqrt(1.0-cosphi*cosphi);
00653       sinphi += 0.24*exp(-sinphi/0.24);
00654       sinphi /= 0.684;  // sin(phideg)/sin(thetadeg)
00655 
00656       // azimuthal angle
00657       if( dirX*dirX+dirY*dirY>0.0 ){
00658         ax = (px*dirZ-pz*dirX) - (py*dirX-px*dirY)*(1.0-dirZ)*dirY/sqrt(dirX*dirX+dirY*dirY);
00659         ay = (py*dirZ-pz*dirY) - (px*dirY-py*dirX)*(1.0-dirZ)*dirX/sqrt(dirX*dirX+dirY*dirY);
00660         az = pz*dirZ + py*dirY + px*dirX;
00661       }
00662       else{
00663         ax = px;
00664         ay = py;
00665         az = pz;
00666       }
00667 
00668       azideg = atan2(ay,ax)/(TMath::Pi()/180.0); // radians->degrees
00669     }
00670 
00671     Double_t Lpoint = ds;
00672     Double_t Ltrack = 0.0;
00673     Double_t Lphoton = 0.0;
00674     Double_t Lscatter = 0.0;
00675  
00676     if( phi<theta ){
00677       Ltrack = Lpoint*sin(theta-phi)/sin(theta);
00678       Lphoton = Lpoint*sin(phi)/sin(theta);
00679       Lscatter = 0.0;
00680     }
00681     else{
00682       Ltrack = 0.0;
00683       Lphoton = Lpoint;
00684       Lscatter = Lpoint*(phi-theta);
00685     }
00686 
00687     Double_t fC = WCSimParameters::SpeedOfLight();
00688     Double_t fN = WCSimParameters::RefractiveIndex(Lphoton);
00689 
00690     Double_t dt = fDigitT[idigit] - vtxTime;
00691     Double_t qpes = fDigitQ[idigit];
00692     Double_t res = WCSimParameters::TimeResolution(qpes);
00693 
00694     fConeAngle[idigit] = thetadeg; // degrees
00695     fZenith[idigit] = phideg;      // degrees
00696     fAzimuth[idigit] = azideg;     // degrees
00697     fSolidAngle[idigit] = sinphi;
00698 
00699     fDistPoint[idigit] = Lpoint;       
00700     fDistTrack[idigit] = Ltrack;
00701     fDistPhoton[idigit] = Lphoton;
00702     fDistScatter[idigit] = Lscatter;
00703 
00704     fDeltaTime[idigit] = dt;
00705     fDeltaSigma[idigit] = res;   
00706     
00707     fDeltaAngle[idigit] = phideg-thetadeg; // degrees
00708     fDeltaPoint[idigit] = Lpoint/(fC/fN);
00709     fDeltaTrack[idigit] = Ltrack/fC;
00710     fDeltaPhoton[idigit] = Lphoton/(fC/fN);
00711     fDeltaScatter[idigit] = Lscatter/(fC/fN);
00712  
00713     fPointPath[idigit] = fN*Lpoint;
00714     fExtendedPath[idigit] = Ltrack + fN*Lphoton;
00715 
00716     fPointResidual[idigit] = dt - Lpoint/(fC/fN);
00717     fExtendedResidual[idigit] = dt - Ltrack/fC - Lphoton/(fC/fN);
00718 
00719     fDelta[idigit] = fExtendedResidual[idigit]; // default
00720   }
00721 
00722   return;
00723 }   
00724 
00725 Double_t WCSimVertexGeometry::GetDeltaCorrection(Int_t idigit, Double_t Length)
00726 {
00727   if( Length<=0.0
00728    || GetDistTrack(idigit)<Length ){
00729     return 0.0;
00730   }
00731 
00732   else{
00733     Double_t Lpoint = GetDistPoint(idigit);
00734     Double_t Ltrack = GetDistTrack(idigit);
00735     Double_t Lphoton = GetDistPhoton(idigit);
00736     Double_t AngleRad = (TMath::Pi()/180.0)*GetAngle(idigit);    
00737     Double_t ConeAngleRad = (TMath::Pi()/180.0)*GetConeAngle(idigit);  
00738 
00739     Double_t LtrackNew = Length;
00740     Double_t LphotonNew = sqrt( Lpoint*Lpoint + Length*Length
00741                                 -2.0*Lpoint*Length*cos(AngleRad) );
00742 
00743     Double_t theta = ConeAngleRad;
00744     Double_t sinphi = (Lpoint/LphotonNew)*sin(AngleRad);
00745     Double_t phi = asin(sinphi);
00746     Double_t alpha = theta-phi;
00747     Double_t LphotonNewCorrected = LphotonNew*alpha/sin(alpha);
00748 
00749     Double_t fC = WCSimParameters::SpeedOfLight();
00750     Double_t fN = WCSimParameters::RefractiveIndex(Lpoint);
00751 
00752     return ( Ltrack/fC + Lphoton/(fC/fN) )
00753          - ( LtrackNew/fC + LphotonNewCorrected/(fC/fN) );
00754   }
00755 }
00756 
00757 void WCSimVertexGeometry::CalcVertexSeeds(WCSimRecoEvent* myEvent, Int_t NSeeds)
00758 {
00759   // load event
00760   // ==========
00761   this->LoadEvent(myEvent);
00762 
00763   // calculate vertex seeds
00764   // ======================
00765   return this->CalcVertexSeeds(NSeeds);
00766 }
00767   
00768 void WCSimVertexGeometry::CalcVertexSeeds(Int_t NSeeds)
00769 {
00770   // reset list of seeds
00771   // ===================
00772   vSeedVtxX.clear();
00773   vSeedVtxY.clear();
00774   vSeedVtxZ.clear();
00775   vSeedVtxTime.clear();
00776 
00777   // always calculate the simple vertex
00778   // ==================================
00779   this->CalcSimpleVertex(fVtxX1,fVtxY1,fVtxZ1,fVtxTime1);
00780 
00781   // add this vertex
00782   vSeedVtxX.push_back(fVtxX1); 
00783   vSeedVtxY.push_back(fVtxY1);
00784   vSeedVtxZ.push_back(fVtxZ1);
00785   vSeedVtxTime.push_back(fVtxTime1);
00786 
00787   // check limit
00788   if( NSeeds<=1 ) return;
00789 
00790 
00791   // form list of golden digits
00792   // ==========================
00793   vSeedDigitList.clear();
00794   
00795   for( fThisDigit=0; fThisDigit<fNDigits; fThisDigit++ ){
00796     if( fIsFiltered[fThisDigit] ){
00797       if( fDigitT[fThisDigit] - fMinTime>=0
00798        && fDigitT[fThisDigit] - fMinTime<300
00799        && fDigitQ[fThisDigit]>=5.0 ){
00800         vSeedDigitList.push_back(fThisDigit);
00801       }
00802     }
00803   }
00804 
00805   // check for enough digits
00806   if( vSeedDigitList.size()<=4 ) return;
00807 
00808 
00809   // generate new list of seeds
00810   // ==========================
00811  
00812   Double_t x0 = 0.0;
00813   Double_t y0 = 0.0;
00814   Double_t z0 = 0.0;
00815   Double_t t0 = 0.0;
00816 
00817   Double_t x1 = 0.0;
00818   Double_t y1 = 0.0;
00819   Double_t z1 = 0.0;
00820   Double_t t1 = 0.0; 
00821 
00822   Double_t x2 = 0.0;
00823   Double_t y2 = 0.0;
00824   Double_t z2 = 0.0;
00825   Double_t t2 = 0.0;
00826 
00827   Double_t x3 = 0.0;
00828   Double_t y3 = 0.0;
00829   Double_t z3 = 0.0;
00830   Double_t t3 = 0.0;
00831 
00832   UInt_t counter = 0;
00833   UInt_t NSeedsTarget = NSeeds;
00834 
00835   while( GetNSeeds()<NSeedsTarget && counter<100*NSeedsTarget ){
00836     counter++;
00837 
00838     // choose next four digits
00839     this->ChooseNextQuadruple(x0,y0,z0,t0,
00840                               x1,y1,z1,t1,
00841                               x2,y2,z2,t2,
00842                               x3,y3,z3,t3);
00843        
00844     //
00845     //std::cout << "   digit0: (x,y,z,t)=(" << x0 << "," << y0 << "," << z0 << "," << t0 << ") " << std::endl;
00846     //std::cout << "   digit1: (x,y,z,t)=(" << x1 << "," << y1 << "," << z1 << "," << t1 << ") " << std::endl;
00847     //std::cout << "   digit2: (x,y,z,t)=(" << x2 << "," << y2 << "," << z2 << "," << t2 << ") " << std::endl;
00848     //std::cout << "   digit3: (x,y,z,t)=(" << x3 << "," << y3 << "," << z3 << "," << t3 << ") " << std::endl;
00849     //
00850 
00851     // find common vertex
00852     WCSimGeometry::FindVertex(x0,y0,z0,t0,
00853                               x1,y1,z1,t1,
00854                               x2,y2,z2,t2,
00855                               x3,y3,z3,t3,
00856                               fVtxX1,fVtxY1,fVtxZ1,fVtxTime1,
00857                               fVtxX2,fVtxY2,fVtxZ2,fVtxTime2);
00858 
00859     //
00860     // std::cout << "   result: (x,y,z,t)=(" << fVtxX1 << "," << fVtxY1 << "," << fVtxZ1 << "," << fVtxTime1 << ") " << std::endl
00861     //           << "   result: (x,y,z,t)=(" << fVtxX2 << "," << fVtxY2 << "," << fVtxZ2 << "," << fVtxTime2 << ") " << std::endl;
00862     //
00863 
00864     // add first digit
00865     if( WCSimGeometry::Instance()->InsideDetector(fVtxX1,fVtxY1,fVtxZ1) ){
00866       vSeedVtxX.push_back(fVtxX1); 
00867       vSeedVtxY.push_back(fVtxY1);
00868       vSeedVtxZ.push_back(fVtxZ1);
00869       vSeedVtxTime.push_back(fVtxTime1);
00870     }
00871 
00872     // add second digit
00873     if( WCSimGeometry::Instance()->InsideDetector(fVtxX2,fVtxY2,fVtxZ2) ){
00874       vSeedVtxX.push_back(fVtxX2); 
00875       vSeedVtxY.push_back(fVtxY2);
00876       vSeedVtxZ.push_back(fVtxZ2);
00877       vSeedVtxTime.push_back(fVtxTime2);
00878     }
00879 
00880   }
00881 
00882   return;
00883 }
00884 
00885 void WCSimVertexGeometry::ChooseNextQuadruple(Double_t& x0, Double_t& y0, Double_t& z0, Double_t& t0, Double_t& x1, Double_t& y1, Double_t& z1, Double_t& t1, Double_t& x2, Double_t& y2, Double_t& z2, Double_t& t2, Double_t& x3, Double_t& y3, Double_t& z3, Double_t& t3)
00886 {
00887   this->ChooseNextDigit(x0,y0,z0,t0);
00888   this->ChooseNextDigit(x1,y1,z1,t1);
00889   this->ChooseNextDigit(x2,y2,z2,t2);
00890   this->ChooseNextDigit(x3,y3,z3,t3);
00891 
00892   return;
00893 }
00894 
00895 void WCSimVertexGeometry::ChooseNextDigit(Double_t& xpos, Double_t& ypos, Double_t& zpos, Double_t& time)
00896 {
00897   // default
00898   xpos=0; ypos=0; zpos=0; time=0;
00899 
00900   // ROOT random number generator
00901   // Double_t r = gRandom->Rndm();
00902 
00903   // pseudo-random number generator
00904   Int_t numEntries = vSeedDigitList.size();
00905 
00906   fCounter++;
00907   if( fCounter>=fNDigits ) fCounter = 0;
00908   fThisDigit = vSeedDigitList.at(fLastEntry);
00909 
00910   Double_t t0 = 0.5 + fDigitT[fCounter] - fMinTime;
00911   Double_t q0 = 0.5 + fDigitQ[fCounter];
00912 
00913   Double_t t1 = 0.5 + fDigitT[fThisDigit] - fMinTime;
00914   Double_t q1 = 0.5 + fDigitQ[fThisDigit];
00915 
00916   Double_t tq = 100.0*(t0*q0+t1*q1);
00917   Double_t r = tq - TMath::Floor(tq);
00918 
00919   fLastEntry = (Int_t)(r*numEntries);
00920 
00921   // return the new digit
00922   fThisDigit = vSeedDigitList.at(fLastEntry);
00923   xpos = fDigitX[fThisDigit];
00924   ypos = fDigitY[fThisDigit];
00925   zpos = fDigitZ[fThisDigit];
00926   time = fDigitT[fThisDigit];
00927 
00928   return;
00929 }
00930 
00931 
00932