WCSimGeometry.cc

Go to the documentation of this file.
00001 #include "WCSimGeometry.hh"
00002 #include "WCSimParameters.hh"
00003 
00004 #include "TFile.h"
00005 #include "TTree.h"
00006 #include "TDirectory.h"
00007 #include "TMath.h"
00008 #include "TVector3.h"
00009 #include "TMatrixD.h"
00010 #include "TError.h"
00011 
00012 #include <iostream>
00013 #include <cmath>
00014 #include <cassert>
00015 
00016 ClassImp(WCSimGeometry)
00017 
00018 static WCSimGeometry* fgGeometryHandle = 0;
00019 
00020 Bool_t WCSimGeometry::TouchGeometry()
00021 {
00022   if( WCSimGeometry::Instance()->GetNumPMTs()>0 ){
00023     return true;
00024   }
00025   else{
00026     std::cout << " ***  WCSimGeometry::TouchData() *** " << std::endl;
00027     std::cout << "  <error> need to build geometry... " << std::endl
00028               << "    Call: WCSimGeometry::BuildGeometry(WCSimRootGeom*) " << std::endl; 
00029     return false;
00030   }
00031 }
00032 
00033 WCSimGeometry* WCSimGeometry::Instance()
00034 {
00035   if( !fgGeometryHandle ){
00036     fgGeometryHandle = new WCSimGeometry();
00037   }
00038   
00039   return fgGeometryHandle;
00040 }
00041 
00042 void WCSimGeometry::BuildGeometry(WCSimRootGeom* geom)
00043 {
00044   return WCSimGeometry::Instance()->SetGeometry(geom);
00045 }
00046 
00047 void WCSimGeometry::WriteGeometry(const char* filename)
00048 {
00049   WCSimGeometry::Instance()->WriteToFile(filename);
00050 }
00051 
00052 void WCSimGeometry::PrintGeometry()
00053 {
00054   std::cout << " *** WCSimGeometry::PrintGeometry() *** " << std::endl;
00055 
00056   if( WCSimGeometry::TouchGeometry() ){
00057 
00058     // geometry type
00059     if( fgGeometryHandle->GetGeoType()==WCSimGeometry::kCylinder ){
00060       std::cout << "  Detector Geometry: CYLINDER " << std::endl
00061                 << "   radius[m]=" <<  1.0e-2*fgGeometryHandle->GetCylRadius() << std::endl
00062                 << "   length[m]=" << 1.0e-2*fgGeometryHandle->GetCylLength() << std::endl              
00063                 << "   fiducial_radius[m]=" <<  1.0e-2*fgGeometryHandle->GetCylFiducialRadius() << std::endl
00064                 << "   fiducial_length[m]=" << 1.0e-2*fgGeometryHandle->GetCylFiducialLength() << std::endl;
00065       std::cout << "   area[m^2]=" << 1.0e-4*fgGeometryHandle->GetArea() << std::endl
00066                 << "   volume[m^3]=" << 1.0e-6*fgGeometryHandle->GetVolume() << std::endl
00067                 << "   fiducial_volume[m^3]=" << 1.0e-6*fgGeometryHandle->GetFiducialVolume() << std::endl;
00068     }
00069     if( fgGeometryHandle->GetGeoType()==WCSimGeometry::kMailBox ){
00070       std::cout << "  Detector Geometry: MAILBOX " << std::endl
00071                 << "   dx[m]=" << 1.0e-2*fgGeometryHandle->GetMailBoxX() << std::endl
00072                 << "   dy[m]=" << 1.0e-2*fgGeometryHandle->GetMailBoxY() << std::endl
00073                 << "   dz[m]=" << 1.0e-2*fgGeometryHandle->GetMailBoxZ() << std::endl
00074                 << "   fiducial_dx[m]=" << 1.0e-2*fgGeometryHandle->GetMailBoxFiducialX() << std::endl
00075                 << "   fiducial_dy[m]=" << 1.0e-2*fgGeometryHandle->GetMailBoxFiducialY() << std::endl
00076                 << "   fiducial_dz[m]=" << 1.0e-2*fgGeometryHandle->GetMailBoxFiducialZ() << std::endl;
00077       std::cout << "   area[m^2]=" << 1.0e-4*fgGeometryHandle->GetArea() << std::endl
00078                 << "   volume[m^3]=" << 1.0e-6*fgGeometryHandle->GetVolume() << std::endl
00079                 << "   fiducial_volume[m^3]=" << 1.0e-6*fgGeometryHandle->GetFiducialVolume() << std::endl;
00080     }
00081     if( fgGeometryHandle->GetGeoType()==WCSimGeometry::kUnknown ){
00082       std::cout << "  Detector Geometry: UNKNOWN " << std::endl;
00083       std::cout << "   <warning> geometry unknown! " << std::endl;
00084     }  
00085 
00086     // Dteector PMTs
00087     std::cout << "  Detector PMTs: " << std::endl;
00088     std::cout << "   number=" << fgGeometryHandle->GetNumPMTs() << std::endl;
00089     std::cout << "   radius[cm]=" << fgGeometryHandle->GetPMTRadius() << std::endl;
00090     std::cout << "   separation[cm]=" << fgGeometryHandle->GetPMTSeparation() << std::endl;
00091     std::cout << "   coverage=" << fgGeometryHandle->GetPMTCoverage() << std::endl;
00092 
00093     // attempt to figure out configuration...
00094     std::cout << "  Detector Configuration: " << std::endl; 
00095     GeoConfiguration_t config = (GeoConfiguration_t)(fgGeometryHandle->GetGeoConfig());
00096     std::cout << "   " << fgGeometryHandle->AsString(config) << std::endl;
00097   }
00098   
00099   else {
00100     std::cout << " <warning> no geometry loaded " << std::endl;
00101   }
00102 
00103   return;
00104 }
00105 
00106 const char* WCSimGeometry::AsString(GeoConfiguration_t config)
00107 {
00108   switch( config ){
00109     case WCSimGeometry::kDUSEL_Unknown:
00110       return "DUSEL_Unknown";
00111     case WCSimGeometry::kDUSEL_100kton_10inch_12percent:
00112       return "DUSEL_100kton_10inch_12percent";
00113     case WCSimGeometry::kDUSEL_100kton_10inch_15percent:
00114       return "DUSEL_100kton_10inch_15percent";
00115     case WCSimGeometry::kDUSEL_100kton_10inch_30percent:
00116       return "DUSEL_100kton_10inch_30percent";
00117     case WCSimGeometry::kDUSEL_100kton_10inch_40percent:
00118       return "DUSEL_100kton_10inch_40percent";
00119     case WCSimGeometry::kDUSEL_200kton_10inch_12percent:
00120       return "DUSEL_200kton_10inch_12percent";
00121     case WCSimGeometry::kDUSEL_200kton_12inch_10percent:
00122       return "DUSEL_200kton_12inch_10percent";
00123     case WCSimGeometry::kDUSEL_200kton_12inch_14percent:
00124       return "DUSEL_200kton_12inch_14percent";
00125     case WCSimGeometry::kDUSEL_150ktonMailbox_10inch_15percent:
00126       return "DUSEL_150ktonMailbox_10inch_15percent";
00127     case WCSimGeometry::kDUSEL_150ktonMailbox_10inch_30percent:
00128       return "DUSEL_150ktonMailbox_10inch_30percent";
00129     case WCSimGeometry::kDUSEL_150ktonMailbox_10inch_40percent:
00130       return "DUSEL_150ktonMailbox_10inch_40percent";
00131     default:
00132       return "DUSEL_Unknown";
00133   }
00134 }
00135 
00136 WCSimGeometry::WCSimGeometry() :
00137   fPMTs(0),
00138   fPmtX(0),
00139   fPmtY(0),
00140   fPmtZ(0),
00141   fPmtRegion(0),
00142   fWCSimRootGeom(0)
00143 {
00144   this->SetGeometry(0);
00145 }
00146 
00147 WCSimGeometry::WCSimGeometry(WCSimRootGeom* geom) :
00148   fPMTs(0),
00149   fPmtX(0),
00150   fPmtY(0),
00151   fPmtZ(0),
00152   fPmtRegion(0),
00153   fWCSimRootGeom(0)
00154 {
00155   this->SetGeometry(geom);
00156 }
00157 
00158 WCSimGeometry::~WCSimGeometry()
00159 {
00160   if( fPmtX )         delete [] fPmtX;       
00161   if( fPmtY )         delete [] fPmtY;       
00162   if( fPmtZ )         delete [] fPmtZ;       
00163   if( fPmtRegion )    delete [] fPmtRegion;  
00164 }
00165  
00166 Int_t WCSimGeometry::GetRegion( Int_t tube )
00167 {
00168   if( tube>=0 && tube<fPMTs ){
00169     return fPmtRegion[tube];
00170   }
00171 
00172   return -1;
00173 }
00174 
00175 Double_t WCSimGeometry::GetX( Int_t tube )  
00176 {
00177   if( tube>=0 && tube<fPMTs ){
00178     return fPmtX[tube];
00179   }
00180 
00181   return -99999.9;
00182 }
00183 
00184 Double_t WCSimGeometry::GetY( Int_t tube )
00185 {
00186   if( tube>=0 && tube<fPMTs ){
00187     return fPmtY[tube];
00188   } 
00189 
00190   return -99999.9;
00191 }
00192   
00193 Double_t WCSimGeometry::GetZ( Int_t tube )
00194 {
00195   if( tube>=0 && tube<fPMTs ){
00196     return fPmtZ[tube];
00197   } 
00198 
00199   return -99999.9;
00200 }
00201   
00202 void WCSimGeometry::Reset()
00203 {
00204   if( WCSimGeometry::Instance()->GetNumPMTs()>0 ){
00205     std::cout << " *** WCSimGeometry::Reset() *** " << std::endl;
00206   }
00207 
00208   WCSimGeometry::Instance()->SetGeometry(0);
00209 }
00210 
00211 void WCSimGeometry::SetGeometry(WCSimRootGeom* myGeometry)
00212 {  
00213   // Store WCSim Geometry
00214   // ====================
00215   fWCSimRootGeom = myGeometry;
00216 
00217   // Reset Geometry
00218   // ==============
00219   fPMTs = 0;
00220 
00221   if( fPmtX ) delete [] fPmtX;             fPmtX = 0;
00222   if( fPmtY ) delete [] fPmtY;             fPmtY = 0;
00223   if( fPmtZ ) delete [] fPmtZ;             fPmtZ = 0;
00224   if( fPmtRegion ) delete [] fPmtRegion;   fPmtRegion = 0;
00225 
00226   // Reset Geometry
00227   // ==============
00228   fGeoType = WCSimGeometry::kUnknown;
00229   fGeoConfig = WCSimGeometry::kDUSEL_Unknown;
00230 
00231   fCylinder = 0;
00232   fCylRadius = 0.0;
00233   fCylLength = 0.0;  
00234   fCylFiducialRadius = 0.0;
00235   fCylFiducialLength = 0.0;
00236     
00237   fMailBox = 0;
00238   fMailBoxX = 0.0;
00239   fMailBoxY = 0.0;
00240   fMailBoxZ = 0.0;  
00241   fMailBoxFiducialX = 0.0;
00242   fMailBoxFiducialY = 0.0;
00243   fMailBoxFiducialZ = 0.0;
00244 
00245   fDetectorArea = 0.0;
00246   fDetectorVolume = 0.0;
00247   fDetectorFiducialVolume = 0.0;
00248 
00249   // Sanity Check
00250   // ============
00251   if( myGeometry == 0 ){
00252     return;
00253   }
00254 
00255   // Building New Geometry
00256   // =====================
00257   std::cout << " *** WCSimGeometry::BuildGeometry(...) *** " << std::endl;
00258 
00259   // Get Geometry Type
00260   // =================
00261   fGeoTypeInput = myGeometry->GetGeo_Type();
00262 
00263   std::cout << "   reading geometry: type=" << fGeoTypeInput << std::endl;
00264 
00265   // --- hack ---
00266   if( fGeoTypeInput<0 || fGeoTypeInput>1 ){
00267     std::cout << "   <warning> Detected Unknown Geometry: " << fGeoTypeInput << std::endl;
00268     std::cout << "             Horrible hack! Assume it's a cylinder..." << std::endl;
00269     std::cout << "             Setting GeoType = 0 " << std::endl;
00270     fGeoTypeInput = 0;
00271   }
00272 
00273   if( fGeoTypeInput==0 ){  // geotype=0: cylinder
00274     std::cout << "   building cylindrical geometry " << std::endl;
00275     fGeoType = WCSimGeometry::kCylinder;
00276     fCylinder = 1;
00277     fCylRadius = myGeometry->GetWCCylRadius();
00278     fCylLength = myGeometry->GetWCCylLength();
00279     if( fCylRadius<3000.0 ){                      // hack for 100kton -> 200kton transition
00280       fCylFiducialRadius = fCylRadius - 250.0;    // note: fiducial volume formerly defined
00281     }                                             //       as 2.5m from cylinder edge, but
00282     else{                                         //       is now reduced to 2.0m from edge.
00283       fCylFiducialRadius = fCylRadius - 200.0;    //                 
00284     }                                             //
00285     fCylFiducialLength = fCylLength - 400.0;      //
00286     fDetectorArea = 2.0*TMath::Pi()*fCylRadius*(fCylRadius
00287                                               + fCylLength);
00288     fDetectorVolume = TMath::Pi()*fCylRadius*fCylRadius*fCylLength;
00289     fDetectorFiducialVolume = TMath::Pi()*fCylFiducialRadius*fCylFiducialRadius*fCylFiducialLength;
00290   }
00291   
00292   if( fGeoTypeInput==1 ){  // geotype=1: mailbox
00293     std::cout << "   building mailbox geometry " << std::endl;
00294     fGeoType = WCSimGeometry::kMailBox;
00295     fMailBox = 1;
00296     fMailBoxX = myGeometry->GetMailBox_x();
00297     fMailBoxY = myGeometry->GetMailBox_y();
00298     fMailBoxZ = myGeometry->GetMailBox_z();
00299     fMailBoxFiducialX = fMailBoxX - 400.0;          // note: fiducial volume defined
00300     fMailBoxFiducialY = fMailBoxY - 400.0;          //       2.0m from mailbox edge
00301     fMailBoxFiducialZ = fMailBoxZ - 400.0;
00302     fDetectorArea = 2.0*(fMailBoxX*fMailBoxY 
00303                        + fMailBoxY*fMailBoxZ 
00304                        + fMailBoxZ*fMailBoxX);
00305     fDetectorVolume = fMailBoxX*fMailBoxY*fMailBoxZ;
00306     fDetectorFiducialVolume = fMailBoxFiducialX*fMailBoxFiducialY*fMailBoxFiducialZ;
00307   }
00308 
00309   // Look up PMT Locations (NB: +1 to stop overflow)
00310   // ===============================================
00311   std::cout << "   reading PMTs: " << myGeometry->GetWCNumPMT() << std::endl;
00312   fPMTs = myGeometry->GetWCNumPMT() + 1;
00313 
00314   fPmtX = new Double_t[fPMTs];
00315   fPmtY = new Double_t[fPMTs];
00316   fPmtZ = new Double_t[fPMTs];
00317   fPmtRegion = new Int_t[fPMTs];
00318 
00319   for( Int_t n=0; n<fPMTs; n++ ){
00320     fPmtX[n] = -99999.9;
00321     fPmtY[n] = -99999.9;
00322     fPmtZ[n] = -99999.9;
00323     fPmtRegion[n] = -1;
00324   }
00325 
00326   for( Int_t n=0; n<myGeometry->GetWCNumPMT(); n++ ){
00327 
00328     // get next PMT
00329     WCSimRootPMT myPmt = myGeometry->GetPMT(n);
00330 
00331     Int_t tube = myPmt.GetTubeNo();
00332     Int_t region = -1; // myPmt.GetCylLoc();
00333     Double_t x = myPmt.GetPosition(0);
00334     Double_t y = myPmt.GetPosition(1);
00335     Double_t z = myPmt.GetPosition(2);
00336 
00337     // hack the region (for now)
00338     Double_t delta = 10.0; // cm
00339 
00340     if( fCylinder ){
00341       if( z>=-0.5*fCylLength+delta 
00342        && z<=+0.5*fCylLength-delta ) region = WCSimGeometry::kSide;
00343       if( z<-0.5*fCylLength+delta )  region = WCSimGeometry::kBottom;
00344       if( z>+0.5*fCylLength-delta )  region = WCSimGeometry::kTop;
00345     }
00346 
00347     if( fMailBox ){
00348       //if( z>=-0.5*fMailBoxZ+delta 
00349       // && z<=+0.5*fMailBoxZ-delta ) region = WCSimGeometry::kSide;
00350       if( z<-0.5*fMailBoxZ+delta )  region = WCSimGeometry::kBottom;
00351       if( z>+0.5*fMailBoxZ-delta )  region = WCSimGeometry::kTop;
00352       if( x<-0.5*fMailBoxX+delta )  region = WCSimGeometry::kBack;
00353       if( x>+0.5*fMailBoxX-delta )  region = WCSimGeometry::kFront;
00354       if( y<-0.5*fMailBoxY+delta )  region = WCSimGeometry::kRight;
00355       if( y>+0.5*fMailBoxY-delta )  region = WCSimGeometry::kLeft;
00356     }
00357 
00358     // add to array
00359     if( tube>=0 && tube<fPMTs ){
00360       fPmtX[tube] = x;
00361       fPmtY[tube] = y;
00362       fPmtZ[tube] = z;
00363       fPmtRegion[tube] = region;
00364     }
00365   }
00366 
00367   // Calculate PMT Coverage
00368   // ======================
00369 
00370   // horrible hack (WCSim switched units on me!)   
00371   if( myGeometry->GetWCPMTRadius()<100.0 ){
00372     fPMTRadius = 1.0*myGeometry->GetWCPMTRadius(); // [must be cm already]
00373   }
00374   else{
00375     fPMTRadius = 0.1*myGeometry->GetWCPMTRadius(); // [must be mm, so convert to cm]
00376   }
00377 
00378   fPMTSurfaceArea = 0.0;
00379   fPMTCoverage = 0.0;
00380   fPMTSeparation = 0.0;
00381 
00382   if( fPMTs>0 && fDetectorArea>0.0 ){
00383     fPMTSurfaceArea = fPMTs*TMath::Pi()*fPMTRadius*fPMTRadius;
00384     fPMTCoverage = fPMTSurfaceArea/fDetectorArea;
00385     fPMTSeparation = sqrt(fDetectorArea/fPMTs);
00386   }
00387 
00388   // Determine Detector Configuration
00389   // ================================
00390   std::cout << "   determing detector configuration " << std::endl;
00391   
00392   Double_t kton = 0.5+1.0e-9*fDetectorFiducialVolume;
00393   Int_t kton_rounded = (Int_t)(10.0*(Int_t)(0.1*kton));
00394   Double_t pmtdiam = 0.5+(2.0/2.54)*fPMTRadius;
00395   Int_t pmtdiam_rounded = (Int_t)(2.0*(Int_t)(0.5*pmtdiam));
00396   Double_t coverage = 0.5 + 100.0*fPMTCoverage;
00397   Int_t coverage_rounded = (Int_t)(coverage);
00398 
00399   std::cout << "    DUSEL_" << kton_rounded << "kton"; 
00400   if( fgGeometryHandle->GetGeoType()==WCSimGeometry::kMailBox ) std::cout << "Mailbox";
00401   std::cout << "_";
00402   std::cout << pmtdiam_rounded
00403             << "inch_" << coverage_rounded
00404             << "perCent " << std::endl;
00405 
00406   // cylindrical geometry
00407   if( fGeoType==WCSimGeometry::kCylinder ){
00408     if( kton_rounded==100 && pmtdiam_rounded==10 && coverage_rounded==12 ) fGeoConfig = WCSimGeometry::kDUSEL_100kton_10inch_12percent;    
00409     if( kton_rounded==100 && pmtdiam_rounded==10 && coverage_rounded==15 ) fGeoConfig = WCSimGeometry::kDUSEL_100kton_10inch_15percent;    
00410     if( kton_rounded==100 && pmtdiam_rounded==10 && coverage_rounded==30 ) fGeoConfig = WCSimGeometry::kDUSEL_100kton_10inch_30percent;  
00411     if( kton_rounded==100 && pmtdiam_rounded==10 && coverage_rounded==40 ) fGeoConfig = WCSimGeometry::kDUSEL_100kton_10inch_40percent; 
00412     if( kton_rounded==200 && pmtdiam_rounded==10 && coverage_rounded==12 ) fGeoConfig = WCSimGeometry::kDUSEL_200kton_10inch_12percent; 
00413     if( kton_rounded==200 && pmtdiam_rounded==12 && coverage_rounded==10 ) fGeoConfig = WCSimGeometry::kDUSEL_200kton_12inch_10percent;
00414     if( kton_rounded==200 && pmtdiam_rounded==12 && coverage_rounded==14 ) fGeoConfig = WCSimGeometry::kDUSEL_200kton_12inch_14percent;
00415   }
00416 
00417   // mailbox geometry
00418   else if( fGeoType==WCSimGeometry::kMailBox ){
00419     if( kton_rounded==150 && pmtdiam_rounded==10 && coverage_rounded==15 ) fGeoConfig = WCSimGeometry::kDUSEL_150ktonMailbox_10inch_15percent;
00420     if( kton_rounded==150 && pmtdiam_rounded==10 && coverage_rounded==30 ) fGeoConfig = WCSimGeometry::kDUSEL_150ktonMailbox_10inch_30percent;
00421     if( kton_rounded==150 && pmtdiam_rounded==10 && coverage_rounded==40 ) fGeoConfig = WCSimGeometry::kDUSEL_150ktonMailbox_10inch_40percent;
00422   }
00423 
00424   // --- sanity check ---
00425   else{
00426     std::cout << " <error> unknown geometry" << std::endl;
00427     assert(0);
00428   }
00429 
00430   // Print Geometry
00431   // ==============
00432   PrintGeometry();
00433 
00434   return;
00435 }
00436 
00437 void WCSimGeometry::WriteToFile(const char* filename)
00438 {
00439   std::cout << " *** WCSimGeometry::WriteGeometry() *** " << std::endl;
00440 
00441   if( fPMTs<=0 ){
00442     std::cout << "  warning: no geometry loaded " << std::endl;
00443     return;
00444   }
00445 
00446   std::cout << "  writing geometry to: "  << filename << std::endl;
00447 
00448   TDirectory* tmpd = 0;
00449   TFile* fFileGeometry = 0;
00450   TTree* fTreeGeometry = 0;
00451 
00452   if( !fFileGeometry ){
00453     tmpd = gDirectory;
00454     fFileGeometry = new TFile(filename,"recreate");
00455     fTreeGeometry = new TTree("ntuple","detector geometry");
00456     fTreeGeometry->Branch("tube",&fTube,"tube/I");
00457     fTreeGeometry->Branch("region",&fRegion,"region/I");
00458     fTreeGeometry->Branch("x",&fXpos,"x/D");
00459     fTreeGeometry->Branch("y",&fYpos,"y/D");
00460     fTreeGeometry->Branch("z",&fZpos,"z/D");
00461     fTreeGeometry->SetAutoSave(100);
00462     gDirectory = tmpd;  
00463   }
00464   
00465   for( Int_t n =0; n<fPMTs; n++ ){    
00466     fTube = n;
00467     fRegion = fPmtRegion[n];
00468     fXpos = fPmtX[n];
00469     fYpos = fPmtY[n];
00470     fZpos = fPmtZ[n];
00471 
00472     if( fRegion>=0 ){
00473       tmpd = gDirectory;
00474       fFileGeometry->cd();
00475       fTreeGeometry->Fill();
00476       gDirectory = tmpd;
00477     }
00478   }
00479  
00480   tmpd = gDirectory;
00481   fFileGeometry->cd();
00482   fTreeGeometry->Write();
00483   fFileGeometry->Close();
00484   fFileGeometry = 0;
00485   fTreeGeometry = 0;
00486   gDirectory = tmpd;
00487 
00488   return;
00489 }
00490  
00491 Bool_t WCSimGeometry::InsideDetector(Double_t x, Double_t y, Double_t z)
00492 {
00493   if( fGeoType==WCSimGeometry::kCylinder ){
00494     if( z>=-0.5*fCylLength && z<=+0.5*fCylLength
00495      && x*x+y*y<=fCylRadius*fCylRadius ){
00496       return 1;
00497     }
00498   }
00499   if( fGeoType==WCSimGeometry::kMailBox ){
00500     if( x>=-0.5*fMailBoxX && x<=+0.5*fMailBoxX
00501      && y>=-0.5*fMailBoxY && y<=+0.5*fMailBoxY
00502      && z>=-0.5*fMailBoxZ && z<=+0.5*fMailBoxZ ){
00503       return 1;
00504     }
00505   }
00506 
00507   return 0;
00508 }
00509 
00510 Bool_t WCSimGeometry::InsideFiducialVolume(Double_t x, Double_t y, Double_t z)
00511 {
00512   if( fGeoType==WCSimGeometry::kCylinder ){
00513     if( z>=-0.5*fCylFiducialLength && z<=+0.5*fCylFiducialLength
00514      && x*x+y*y<=fCylFiducialRadius*fCylFiducialRadius ){
00515       return 1;
00516     }
00517   }
00518   if( fGeoType==WCSimGeometry::kMailBox ){
00519     if( x>=-0.5*fMailBoxFiducialX && x<=+0.5*fMailBoxFiducialX
00520      && y>=-0.5*fMailBoxFiducialY && y<=+0.5*fMailBoxFiducialY
00521      && z>=-0.5*fMailBoxFiducialZ && z<=+0.5*fMailBoxFiducialZ ){
00522       return 1;
00523     }
00524   }
00525 
00526   return 0;
00527 }
00528 
00529 Bool_t WCSimGeometry::InsideDetector(Double_t vx, Double_t vy, Double_t vz, Double_t ex, Double_t ey, Double_t ez)
00530 {
00531   Double_t dx = ex-vx;
00532   Double_t dy = ey-vy;
00533   Double_t dz = ez-vz;
00534   Double_t ds = sqrt(dx*dx+dy*dy+dz*dz);
00535 
00536   Double_t px = dx/ds;
00537   Double_t py = dy/ds;
00538   Double_t pz = dz/ds;
00539 
00540   if( this->ForwardProjectionToEdge(vx,vy,vz,px,py,pz) >= 0.0
00541    && this->BackwardProjectionToEdge(ex,ey,ez,px,py,pz) >= 0.0 ){
00542     return 1;
00543   }
00544   else{
00545     return 0;
00546   }
00547 }
00548 
00549 Double_t WCSimGeometry::DistanceToEdge(Double_t x, Double_t y, Double_t z)
00550 {
00551   // Cylindrical Geometry
00552   // ====================
00553   if( fGeoType==WCSimGeometry::kCylinder ){
00554 
00555     // inside detector (convention: +ve dr)
00556     if( this->InsideDetector(x,y,z) ){
00557       Double_t dr = 0.0;
00558       if( fCylRadius>dr ) dr = fCylRadius;
00559       if( 0.5*fCylLength>dr ) dr = 0.5*fCylLength;
00560       if( -sqrt(x*x+y*y)+fCylRadius<dr ) dr = -sqrt(x*x+y*y)+fCylRadius;
00561       if( -z+0.5*fCylLength<dr ) dr = -z+0.5*fCylLength;
00562       if( +z+0.5*fCylLength<dr ) dr = +z+0.5*fCylLength;
00563       return dr;
00564     }
00565 
00566     // outside detector (convention: -ve dr)
00567     else{
00568 
00569       // side region
00570       if( z>=-0.5*fCylLength && z<=+0.5*fCylLength ){
00571         return -sqrt(x*x+y*y)+fCylRadius;
00572       }
00573 
00574       // top region
00575       if( z<=-0.5*fCylLength
00576        && x*x+y*y<fCylRadius*fCylRadius ){
00577         return +z+0.5*fCylLength;
00578       }
00579       if( z>=+0.5*fCylLength
00580        && x*x+y*y<fCylRadius*fCylRadius ){
00581         return -z+0.5*fCylLength;
00582       }
00583 
00584       // corner regions
00585       if( z>=+0.5*fCylLength
00586        && x*x+y*y>=fCylRadius ){
00587         Double_t dr = sqrt(x*x+y*y)-fCylRadius;
00588         Double_t dz = -z+0.5*fCylLength;
00589         return -sqrt(dr*dr+dz*dz);
00590       }
00591       if( z<=-0.5*fCylLength
00592        && x*x+y*y>=fCylRadius ){
00593         Double_t dr = sqrt(x*x+y*y)-fCylRadius;
00594         Double_t dz = +z+0.5*fCylLength;
00595         return -sqrt(dr*dr+dz*dz);
00596       }
00597     }
00598   }
00599 
00600   // Mailbox Geometry
00601   // ================
00602   if( fGeoType==WCSimGeometry::kMailBox ){
00603 
00604     // inside detector (convention: +ve dr)
00605     if( this->InsideDetector(x,y,z) ){
00606       Double_t dr = 0.0;
00607       if( 0.5*fMailBoxX>dr ) dr = 0.5*fMailBoxX;
00608       if( 0.5*fMailBoxY>dr ) dr = 0.5*fMailBoxY;
00609       if( 0.5*fMailBoxZ>dr ) dr = 0.5*fMailBoxZ;
00610       if( -x+0.5*fMailBoxX<dr ) dr = -x+0.5*fMailBoxX;
00611       if( +x+0.5*fMailBoxX<dr ) dr = +x+0.5*fMailBoxX;
00612       if( -y+0.5*fMailBoxY<dr ) dr = -y+0.5*fMailBoxY;
00613       if( +y+0.5*fMailBoxY<dr ) dr = +y+0.5*fMailBoxY;
00614       if( -z+0.5*fMailBoxZ<dr ) dr = -z+0.5*fMailBoxZ;
00615       if( +z+0.5*fMailBoxZ<dr ) dr = +z+0.5*fMailBoxZ;
00616       return dr;
00617     }
00618 
00619     // outside detector (convention: -ve dr)
00620     else{
00621 
00622       // six edge regions
00623       if( x<=-0.5*fMailBoxX 
00624        && y>=-0.5*fMailBoxY && y<=+0.5*fMailBoxY 
00625        && z>=-0.5*fMailBoxZ && z<=+0.5*fMailBoxZ ){
00626         return +x+0.5*fMailBoxX;
00627       }
00628       if( x>=+0.5*fMailBoxX 
00629        && y>=-0.5*fMailBoxY && y<=+0.5*fMailBoxY  
00630        && z>=-0.5*fMailBoxZ && z<=+0.5*fMailBoxZ ){
00631         return -x+0.5*fMailBoxX;
00632       }
00633       if( x>=-0.5*fMailBoxX && x<=+0.5*fMailBoxX
00634        && y<=-0.5*fMailBoxY 
00635        && z>=-0.5*fMailBoxZ && z<=+0.5*fMailBoxZ ){
00636         return +y+0.5*fMailBoxY;
00637       }
00638       if( x>=-0.5*fMailBoxX && x<=+0.5*fMailBoxX
00639        && y>=+0.5*fMailBoxY 
00640        && z>=-0.5*fMailBoxZ && z<=+0.5*fMailBoxZ ){
00641         return -y+0.5*fMailBoxY;
00642       }
00643       if( x>=-0.5*fMailBoxX && x<=+0.5*fMailBoxX 
00644        && y>=-0.5*fMailBoxY && y<=+0.5*fMailBoxY  
00645        && z<=-0.5*fMailBoxZ ){
00646         return +z+0.5*fMailBoxZ;
00647       }
00648       if( x>=-0.5*fMailBoxX && x<=+0.5*fMailBoxX 
00649        && y>=-0.5*fMailBoxY && y<=+0.5*fMailBoxY 
00650        && z>=+0.5*fMailBoxZ ){
00651         return -z+0.5*fMailBoxZ;
00652       }
00653 
00654       // eight corner regions
00655       if( x<=-0.5*fMailBoxX 
00656        && y<=-0.5*fMailBoxY 
00657        && z<=-0.5*fMailBoxZ ){
00658         Double_t dx = +x+0.5*fMailBoxX;
00659         Double_t dy = +y+0.5*fMailBoxY;
00660         Double_t dz = +z+0.5*fMailBoxZ;
00661         return -sqrt(dx*dx+dy*dy+dz*dz);
00662       }
00663       if( x>=+0.5*fMailBoxX 
00664        && y<=-0.5*fMailBoxY 
00665        && z<=-0.5*fMailBoxZ ){
00666         Double_t dx = -x+0.5*fMailBoxX;
00667         Double_t dy = +y+0.5*fMailBoxY;
00668         Double_t dz = +z+0.5*fMailBoxZ;
00669         return -sqrt(dx*dx+dy*dy+dz*dz);
00670       }
00671       if( x<=-0.5*fMailBoxX 
00672        && y>=+0.5*fMailBoxY 
00673        && z<=-0.5*fMailBoxZ ){
00674         Double_t dx = +x+0.5*fMailBoxX;
00675         Double_t dy = -y+0.5*fMailBoxY;
00676         Double_t dz = +z+0.5*fMailBoxZ;
00677         return -sqrt(dx*dx+dy*dy+dz*dz);
00678       }
00679       if( x<=-0.5*fMailBoxX 
00680        && y<=-0.5*fMailBoxY 
00681        && z>=+0.5*fMailBoxZ ){
00682         Double_t dx = +x+0.5*fMailBoxX;
00683         Double_t dy = +y+0.5*fMailBoxY;
00684         Double_t dz = -z+0.5*fMailBoxZ;
00685         return -sqrt(dx*dx+dy*dy+dz*dz);
00686       }
00687       if( x>=+0.5*fMailBoxX 
00688        && y>=+0.5*fMailBoxY 
00689        && z<=-0.5*fMailBoxZ ){
00690         Double_t dx = -x+0.5*fMailBoxX;
00691         Double_t dy = -y+0.5*fMailBoxY;
00692         Double_t dz = +z+0.5*fMailBoxZ;
00693         return -sqrt(dx*dx+dy*dy+dz*dz);
00694       }
00695       if( x>=+0.5*fMailBoxX 
00696        && y<=-0.5*fMailBoxY 
00697        && z>=+0.5*fMailBoxZ ){
00698         Double_t dx = -x+0.5*fMailBoxX;
00699         Double_t dy = +y+0.5*fMailBoxY;
00700         Double_t dz = -z+0.5*fMailBoxZ;
00701         return -sqrt(dx*dx+dy*dy+dz*dz);
00702       }
00703       if( x<=-0.5*fMailBoxX 
00704        && y>=+0.5*fMailBoxY 
00705        && z>=+0.5*fMailBoxZ ){
00706         Double_t dx = +x+0.5*fMailBoxX;
00707         Double_t dy = -y+0.5*fMailBoxY;
00708         Double_t dz = -z+0.5*fMailBoxZ;
00709         return -sqrt(dx*dx+dy*dy+dz*dz);
00710       }
00711       if( x>=+0.5*fMailBoxX 
00712        && y>=+0.5*fMailBoxY 
00713        && z>=+0.5*fMailBoxZ ){
00714         Double_t dx = -x+0.5*fMailBoxX;
00715         Double_t dy = -y+0.5*fMailBoxY;
00716         Double_t dz = -z+0.5*fMailBoxZ;
00717         return -sqrt(dx*dx+dy*dy+dz*dz);
00718       }
00719     }
00720   }
00721 
00722   return -99999.9;
00723 }
00724 
00725 Double_t WCSimGeometry::ForwardProjectionToEdge(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz)
00726 {
00727   Double_t xproj = 0.0;
00728   Double_t yproj = 0.0;
00729   Double_t zproj = 0.0;
00730   Int_t regionproj = 0;
00731 
00732   this->ProjectToNearEdge(x,y,z,
00733                           px,py,pz,
00734                           xproj,yproj,zproj,
00735                           regionproj);
00736 
00737   if( regionproj > WCSimGeometry::kUnknown ){
00738     return sqrt( (xproj-x)*(xproj-x)
00739                + (yproj-y)*(yproj-y)
00740                + (zproj-z)*(zproj-z) );
00741   }
00742 
00743   this->ProjectToNearEdge(x,y,z,
00744                           -px,-py,-pz,
00745                           xproj,yproj,zproj,
00746                           regionproj);
00747 
00748   if( regionproj > WCSimGeometry::kUnknown ){
00749     return -sqrt( (xproj-x)*(xproj-x)
00750                 + (yproj-y)*(yproj-y)
00751                 + (zproj-z)*(zproj-z) );
00752   }
00753 
00754   return -99999.9;
00755 }
00756 
00757 Double_t WCSimGeometry::BackwardProjectionToEdge(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz)
00758 {
00759   Double_t xproj = 0.0;
00760   Double_t yproj = 0.0;
00761   Double_t zproj = 0.0;
00762   Int_t regionproj = 0;
00763 
00764   this->ProjectToNearEdge(x,y,z,
00765                           -px,-py,-pz,
00766                           xproj,yproj,zproj,
00767                           regionproj);
00768 
00769   if( regionproj > WCSimGeometry::kUnknown ){
00770     return sqrt( (xproj-x)*(xproj-x)
00771                + (yproj-y)*(yproj-y)
00772                + (zproj-z)*(zproj-z) );
00773   }
00774 
00775   this->ProjectToNearEdge(x,y,z,
00776                           px,py,pz,
00777                           xproj,yproj,zproj,
00778                           regionproj);
00779 
00780   if( regionproj > WCSimGeometry::kUnknown ){
00781     return -sqrt( (xproj-x)*(xproj-x)
00782                 + (yproj-y)*(yproj-y)
00783                 + (zproj-z)*(zproj-z) );
00784   }
00785 
00786   return -99999.9;
00787 }
00788 
00789 void WCSimGeometry::ProjectToNearEdge(Double_t x0, Double_t y0, Double_t z0, Double_t px, Double_t py, Double_t pz, Double_t& xproj, Double_t& yproj, Double_t&zproj, Int_t& regionproj)
00790 {
00791   this->ProjectToEdge(0,
00792                       x0,y0,z0,
00793                       px,py,pz,
00794                       xproj,yproj,zproj,
00795                       regionproj);
00796 
00797   //
00798   // this->ProjectToEdgeOld(0,
00799   //                        x0,y0,z0,
00800   //                        px,py,pz,
00801   //                        xproj,yproj,zproj,
00802   //                        regionproj);
00803   //
00804 
00805   return;
00806 }
00807 
00808 void WCSimGeometry::ProjectToFarEdge(Double_t x0, Double_t y0, Double_t z0, Double_t px, Double_t py, Double_t pz, Double_t& xproj, Double_t& yproj, Double_t&zproj, Int_t& regionproj)
00809 {  
00810   this->ProjectToEdge(1,
00811                       x0,y0,z0,
00812                       px,py,pz,
00813                       xproj,yproj,zproj,
00814                       regionproj);
00815 
00816   //
00817   // this->ProjectToEdgeOld(1,
00818   //                        x0,y0,z0,
00819   //                        px,py,pz,
00820   //                        xproj,yproj,zproj,
00821   //                        regionproj);
00822   //
00823   
00824   return;
00825 }
00826 
00827 void WCSimGeometry::ProjectToEdgeOld(Bool_t useFarEdge, Double_t x0, Double_t y0, Double_t z0, Double_t px, Double_t py, Double_t pz, Double_t& xproj, Double_t& yproj, Double_t&zproj, Int_t& regionproj)
00828 {
00829   // default locations
00830   // =================
00831   xproj = -99999.9;
00832   yproj = -99999.9;
00833   zproj = -99999.9;
00834   regionproj = WCSimGeometry::kUnknown;
00835 
00836   // CYLINDRICAL GEOMETRY
00837   // ====================
00838   if( fGeoType==WCSimGeometry::kCylinder ){
00839 
00840     //
00841     // std::cout << " --- ProjectToEdge --- " << std::endl
00842     //           << "  Vtx=(" << x0 << "," << y0 << "," << z0 << ") " << std::endl
00843     //           << "  Dir=(" << px << "," << py << "," << pz << ") --- " << std::endl; 
00844     //
00845 
00846     // projections
00847     // =========== 
00848     Double_t r = fCylRadius;
00849     Double_t L = fCylLength;
00850 
00851     Double_t x = 0.0;
00852     Double_t y = 0.0;
00853     Double_t z = 0.0;
00854     Int_t region = WCSimGeometry::kUnknown;
00855 
00856     Double_t mxy = 0.0;
00857     Double_t mz = 0.0;    
00858     Double_t xlow = 0.0;
00859     Double_t ylow = 0.0;
00860     Double_t zlow = 0.0;
00861     Double_t slow = 0.0;
00862     Double_t xhigh = 0.0;
00863     Double_t yhigh = 0.0;
00864     Double_t zhigh = 0.0;
00865     Double_t shigh = 0.0;
00866 
00867     Bool_t foundProjectionXY = 0;
00868     Bool_t foundProjectionX = 0;
00869     Bool_t foundProjectionY = 0;
00870     Bool_t foundProjectionZ = 0;
00871     Bool_t foundProjection = 0;
00872 
00873     // find radial intersections
00874     // =========================
00875     if( px!=0.0 ){
00876       mxy = py/px;
00877 
00878       // sanity check: is there an intersection?
00879       if( (1.0+mxy*mxy)*r*r - (mxy*x0-y0)*(mxy*x0-y0) > 0.0 ){
00880 
00881         // intersection at low x
00882         xlow = -(x0+mxy*y0) - sqrt( (1.0+mxy*mxy)*r*r - (mxy*x0-y0)*(mxy*x0-y0) );
00883         xlow /= (1.0+mxy*mxy);
00884         ylow = mxy*xlow;
00885 
00886         // intersection at high x
00887         xhigh = -(x0+mxy*y0) + sqrt( (1.0+mxy*mxy)*r*r - (mxy*x0-y0)*(mxy*x0-y0) );
00888         xhigh /= (1.0+mxy*mxy);
00889         yhigh = mxy*xhigh;
00890 
00891         // forward or backward
00892         if( px>=0 ){
00893           if( xlow>=0.0 ) slow = +1.0; else slow = -1.0;
00894           if( xhigh>=0.0 ) shigh = +1.0; else shigh = -1.0;
00895         }
00896         else{
00897           if( xlow<=0.0 ) slow = +1.0; else slow = -1.0;
00898           if( xhigh<=0.0 ) shigh = +1.0; else shigh = -1.0;
00899         }
00900 
00901         // intersection z values
00902         mz = pz/sqrt(px*px+py*py);
00903         zlow = z0 + mz*slow*sqrt(xlow*xlow+ylow*ylow);
00904         zhigh = z0 + mz*shigh*sqrt(xhigh*xhigh+yhigh*yhigh);
00905 
00906         // found radial projection
00907         foundProjectionXY = 1;
00908       }
00909     }
00910 
00911     // travelling along y-z plane
00912     else if( py!=0.0 ) { 
00913 
00914       // sanity check: is there an intersection?
00915       if( r*r-x0*x0 > 0.0 ){ 
00916 
00917         // intersection at low y
00918         xlow = 0.0;
00919         ylow = -sqrt( r*r-x0*x0 ) - y0;
00920 
00921         // intersection at high y
00922         xhigh = 0.0;
00923         yhigh = +sqrt( r*r-x0*x0 ) - y0;
00924 
00925         // forward or backward
00926         if( py>=0 ){
00927           if( ylow>=0.0 ) slow = +1.0; else slow = -1.0;
00928           if( yhigh>=0.0 ) shigh = +1.0; else shigh = -1.0;
00929         }
00930         else{
00931           if( ylow<=0.0 ) slow = +1.0; else slow = -1.0;
00932           if( yhigh<=0.0 ) shigh = +1.0; else shigh = -1.0;
00933         }
00934 
00935         // intersection z values
00936         mz = pz/sqrt(px*px+py*py);
00937         zlow = z0 + mz*slow*sqrt(xlow*xlow+ylow*ylow);
00938         zhigh = z0 + mz*shigh*sqrt(xhigh*xhigh+yhigh*yhigh);
00939 
00940         // found radial projection
00941         foundProjectionXY = 1;
00942       }
00943     }
00944 
00945     // travelling along z-axis
00946     else if( pz!=0.0 ) {     
00947 
00948       // sanity check: is there an intersection?
00949       if( x0*x0+y0*y0 < r*r ){ 
00950         xlow = -x0;
00951         ylow = -y0;
00952         zlow = -L/2.0;
00953 
00954         xhigh = -x0;
00955         yhigh = -y0;
00956         zhigh = +L/2.0;
00957         
00958         // found radial projection
00959         foundProjectionXY = 1;
00960       }
00961     }
00962 
00963     // check projection with cylinder
00964     // ==============================
00965     if( foundProjectionXY ){
00966 
00967       // project backwards along x-axis
00968       if( px<0.0 ){
00969         if( xlow<0 || xhigh<=0 ){
00970           foundProjectionX = 1;
00971         }
00972       }
00973 
00974       // project forwards along x-axis
00975       if( px>0.0 ){
00976         if( xlow>=0 || xhigh>0 ){
00977           foundProjectionX = 1;
00978         }
00979       }
00980 
00981       // project backwards along y-axis
00982       if( py<0.0 ){
00983         if( ylow<0 || yhigh<=0 ){
00984           foundProjectionY = 1;
00985         }
00986       }
00987 
00988       // project forwards along y-axis
00989       if( py>0.0 ){
00990         if( ylow>=0 || yhigh>0 ){
00991           foundProjectionY = 1;
00992         }
00993       }
00994 
00995       // project backwards along z-axis
00996       if( pz<0.0 ){
00997         if( z0>-L/2.0 
00998           && ( zlow<=+L/2.0 || zhigh<=+L/2.0 ) ){
00999           foundProjectionZ = 1;
01000         }
01001       }
01002 
01003       // project forwards along z-axis
01004       if( pz>0 ){
01005         if( z0<+L/2.0 
01006          && ( zlow>=-L/2.0 || zhigh>=-L/2.0 ) ){
01007           foundProjectionZ = 1;
01008         }
01009       }
01010 
01011       // overall projection
01012       if( foundProjectionX 
01013        && foundProjectionY
01014        && foundProjectionZ ){
01015         foundProjection = 1;
01016       }
01017     }
01018 
01019     // apply top/bottom correction
01020     // ===========================
01021     if( foundProjection ){
01022 
01023       if( zlow>+L/2.0 ){
01024         xlow *= (L/2.0-z0)/(zlow-z0);
01025         ylow *= (L/2.0-z0)/(zlow-z0);  
01026         zlow = +L/2.0; 
01027       }
01028 
01029       if( zlow<-L/2.0 ){
01030         xlow *= (-L/2.0-z0)/(zlow-z0);
01031         ylow *= (-L/2.0-z0)/(zlow-z0);
01032         zlow = -L/2.0;
01033       }
01034  
01035       if( zhigh>+L/2.0 ){
01036         xhigh *= (L/2.0-z0)/(zhigh-z0);
01037         yhigh *= (L/2.0-z0)/(zhigh-z0);  
01038         zhigh = +L/2.0; 
01039       }
01040 
01041       if( zhigh<-L/2.0 ){
01042         xhigh *= (-L/2.0-z0)/(zhigh-z0);
01043         yhigh *= (-L/2.0-z0)/(zhigh-z0);
01044         zhigh = -L/2.0;
01045       }
01046 
01047     }
01048 
01049     //
01050     // std::cout << "  FoundProjection = " << foundProjection << std::endl;
01051     // std::cout << "   low(x,y,z) = (" << xlow+x0 << "," << ylow+y0 << "," << zlow << ") " << std::endl; 
01052     // std::cout << "   high(x,y,z) = (" << xhigh+x0 << "," << yhigh+y0 << "," << zhigh << ") " << std::endl;
01053     //
01054 
01055     // choose intersection point
01056     // =========================
01057     if( foundProjection ){
01058 
01059       if( px!=0.0 ){
01060 
01061         // project backwards along x-axis
01062         if( px<0 ){
01063           if( useFarEdge ){
01064             if( xlow<=0.0 ){
01065               x = xlow+x0;
01066               y = ylow+y0;
01067               z = zlow;
01068             }
01069           }
01070           else{
01071             if( xhigh<0.0 ){
01072               x = xhigh+x0;
01073               y = yhigh+y0;
01074               z = zhigh;
01075             }
01076             else if( xlow<=0.0 ){
01077               x = xlow+x0;
01078               y = ylow+y0;
01079               z = zlow;
01080             }
01081           }
01082         }        
01083 
01084         // project forwards along x-axis
01085         if( px>0.0 ){
01086           if( useFarEdge ){
01087             if( xhigh>=0.0 ){
01088               x = xhigh+x0;
01089               y = yhigh+y0;
01090               z = zhigh;
01091             }
01092           }
01093           else{
01094             if( xlow>0.0 ){
01095               x = xlow+x0;
01096               y = ylow+y0;
01097               z = zlow;
01098             }
01099             else if( xhigh>=0.0 ){
01100               x = xhigh+x0;
01101               y = yhigh+y0;
01102               z = zhigh;
01103             }
01104           }
01105         }
01106       }
01107 
01108       else if( py!=0.0 ){
01109 
01110         // project backwards along y-axis
01111         if( py<0.0 ){
01112           if( useFarEdge ){
01113             if( ylow<=0.0 ){
01114               x = xlow+x0;
01115               y = ylow+y0;
01116               z = zlow;
01117             }
01118           }
01119           else{
01120             if( yhigh<0.0 ){
01121               x = xhigh+x0;
01122               y = yhigh+y0;
01123               z = zhigh;
01124             }
01125             else if( ylow<=0.0 ){
01126               x = xlow+x0;
01127               y = ylow+y0;
01128               z = zlow;
01129             }
01130           }
01131         }
01132 
01133         // project forwards along y-axis
01134         if( py>0.0 ){
01135           if( useFarEdge ){
01136             if( yhigh>=0.0 ){
01137               x = xhigh+x0;
01138               y = yhigh+y0;
01139               z = zhigh;
01140             }
01141           }
01142           else{
01143             if( ylow>0.0 ){
01144               x = xlow+x0;
01145               y = ylow+y0;
01146               z = zlow;
01147             }
01148             else if( yhigh>=0.0 ){
01149               x = xhigh+x0;
01150               y = yhigh+y0;
01151               z = zhigh;
01152             }
01153           }
01154         }
01155 
01156       }
01157 
01158       else{
01159 
01160         // project backwards along z-axis
01161         if( pz<0.0 ){
01162           if( useFarEdge ){
01163             if( z0>=-L/2.0 ){
01164               x = xlow+x0;
01165               y = ylow+y0;
01166               z = zlow;
01167             }
01168           }
01169           else{
01170             if( z0>+L/2 ){
01171               x = xhigh+x0;
01172               y = yhigh+y0;
01173               z = zhigh;
01174             }
01175             else if( z0>=-L/2.0 ){
01176               x = xlow+x0;
01177               y = ylow+y0;
01178               z = zlow;
01179             }
01180           }
01181         }
01182 
01183         // project forwards along z-axis
01184         if( pz>0.0 ){
01185           if( useFarEdge ){
01186             if( z0<=+L/2.0 ){
01187               x = xhigh+x0;
01188               y = yhigh+y0;
01189               z = zhigh;
01190             }
01191           }
01192           else{
01193             if( z0<-L/2 ){
01194               x = xlow+x0;
01195               y = ylow+y0;
01196               z = zlow;
01197             }
01198             else if( z0<=+L/2.0 ){
01199               x = xhigh+x0;
01200               y = yhigh+y0;
01201               z = zhigh;
01202             }
01203           }
01204         }
01205 
01206       }
01207 
01208       // assign region
01209       if( -L/2.0<z && z<+L/2.0 ){
01210         region = WCSimGeometry::kSide; // region 1
01211       }
01212 
01213       if( z>=+L/2.0 ){ 
01214         region = WCSimGeometry::kTop; // region 0 
01215       }
01216 
01217       if( z<=-L/2.0 ){ 
01218         region = WCSimGeometry::kBottom; // region 2
01219       }
01220 
01221     }
01222     
01223     //   
01224     // std::cout << "  Return Projection: x=" << x << " y=" << y << " z=" << z << " region=" << region << std::endl;
01225     //
01226 
01227     // return projection
01228     // =================
01229     if( region > WCSimGeometry::kUnknown ){
01230 
01231       xproj = x;
01232       yproj = y;
01233       zproj = z;
01234       regionproj = region;
01235 
01236       return;
01237     }
01238   }
01239 
01240   // MAILBOX GEOMETRY
01241   // ================
01242   if( fGeoType==WCSimGeometry::kMailBox ){
01243 
01244     Double_t dX = fMailBoxX;
01245     Double_t dY = fMailBoxY;
01246     Double_t dZ = fMailBoxZ;
01247 
01248     Double_t x_top    = -99999.9;
01249     Double_t y_top    = -99999.9;
01250     Double_t z_top    = -99999.9;
01251 
01252     Double_t x_bottom = -99999.9;
01253     Double_t y_bottom = -99999.9;
01254     Double_t z_bottom = -99999.9;
01255 
01256     Double_t x_front  = -99999.9;
01257     Double_t y_front  = -99999.9;
01258     Double_t z_front  = -99999.9;    
01259 
01260     Double_t x_back   = -99999.9;
01261     Double_t y_back   = -99999.9;
01262     Double_t z_back   = -99999.9;   
01263 
01264     Double_t x_right  = -99999.9;
01265     Double_t y_right  = -99999.9;
01266     Double_t z_right  = -99999.9;
01267 
01268     Double_t x_left   = -99999.9;
01269     Double_t y_left   = -99999.9;
01270     Double_t z_left   = -99999.9;  
01271 
01272     Double_t x = 0.0;
01273     Double_t y = 0.0;
01274     Double_t z = 0.0;
01275     Int_t region = WCSimGeometry::kUnknown;
01276   
01277     // calculate projections in each plane
01278     // ===================================
01279     if( pz!=0.0 ){
01280       x_top    = x0 + (px/pz)*(+0.5*dZ-z0);
01281       y_top    = y0 + (py/pz)*(+0.5*dZ-z0);
01282       z_top    = +0.5*dZ;
01283 
01284       x_bottom = x0 + (px/pz)*(-0.5*dZ-z0);
01285       y_bottom = y0 + (py/pz)*(-0.5*dZ-z0);
01286       z_bottom = -0.5*dZ;
01287     }
01288 
01289     if( py!=0.0 ){
01290       z_left  = z0 + (pz/py)*(+0.5*dY-y0);
01291       x_left  = x0 + (px/py)*(+0.5*dY-y0);   
01292       y_left  = +0.5*dY;
01293     
01294       z_right   = z0 + (pz/py)*(-0.5*dY-y0);
01295       x_right   = x0 + (px/py)*(-0.5*dY-y0);
01296       y_right   = -0.5*dY;
01297     }
01298 
01299     if( px!=0.0 ){
01300       y_front  = y0 + (py/px)*(+0.5*dX-x0);
01301       z_front  = z0 + (pz/px)*(+0.5*dX-x0);
01302       x_front  = +0.5*dX;
01303 
01304       y_back   = y0 + (py/px)*(-0.5*dX-x0);
01305       z_back   = z0 + (pz/px)*(-0.5*dX-x0);
01306       x_back   = -0.5*dX;
01307     }
01308 
01309     // calculate intersections with box
01310     // ================================
01311     if( x_top>=-0.5*dX && x_top<=+0.5*dX
01312      && y_top>=-0.5*dY && y_top<=+0.5*dY ){
01313       if( pz<0.0 && z0>z_top ){
01314         if( !useFarEdge ){
01315           x = x_top;  
01316           y = y_top;  
01317           z = z_top;
01318           region = WCSimGeometry::kTop;
01319         }
01320       }
01321       if( pz>0.0 && z0<=z_top ){
01322         if( z0>z_bottom
01323          || useFarEdge ){
01324           x = x_top;  
01325           y = y_top;  
01326           z = z_top;
01327           region = WCSimGeometry::kTop;
01328         }
01329       }
01330     }   
01331 
01332     if( x_bottom>=-0.5*dX && x_bottom<=+0.5*dX
01333      && y_bottom>=-0.5*dY && y_bottom<=+0.5*dY ){
01334       if( pz>0.0 && z0<z_bottom ){
01335         if( !useFarEdge ){
01336           x = x_bottom;  
01337           y = y_bottom;  
01338           z = z_bottom;
01339           region = WCSimGeometry::kBottom;
01340         }
01341       }
01342       if( pz<0.0 && z0>=z_bottom ){
01343         if( z0<z_top
01344          || useFarEdge ){
01345           x = x_bottom;  
01346           y = y_bottom;  
01347           z = z_bottom;
01348           region = WCSimGeometry::kBottom;
01349         }
01350       }
01351     }
01352 
01353     if( y_front>=-0.5*dY && y_front<=+0.5*dY
01354      && z_front>=-0.5*dZ && z_front<=+0.5*dZ ){
01355       if( px<0.0 && x0>x_front ){
01356         if( !useFarEdge ){
01357           x = x_front;  
01358           y = y_front;  
01359           z = z_front;
01360           region = WCSimGeometry::kFront;
01361         }
01362       }
01363       if( px>0.0 && x0<=x_front ){
01364         if( x0>x_back
01365          || useFarEdge ){
01366           x = x_front;  
01367           y = y_front;  
01368           z = z_front;
01369           region = WCSimGeometry::kFront;
01370         }
01371       }
01372 
01373     }       
01374 
01375     if( y_back>=-0.5*dY && y_back<=+0.5*dY
01376      && z_back>=-0.5*dZ && z_back<=+0.5*dZ ){
01377       if( px>0.0 && x0<x_back ){
01378         if( !useFarEdge ){
01379           x = x_back;  
01380           y = y_back;  
01381           z = z_back;
01382           region = WCSimGeometry::kBack;
01383         }
01384       }
01385       if( px<0.0 && x0>=x_back ){
01386         if( x0<x_front
01387          || useFarEdge ){
01388           x = x_back;  
01389           y = y_back;  
01390           z = z_back;
01391           region = WCSimGeometry::kBack;
01392         }
01393       }
01394     }     
01395 
01396     if( z_left>=-0.5*dZ && z_left<=+0.5*dZ
01397      && x_left>=-0.5*dX && x_left<=+0.5*dX ){
01398       if( py<0.0 && y0>y_left ){
01399         if( !useFarEdge ){
01400           x = x_left;  
01401           y = y_left;  
01402           z = z_left;
01403           region = WCSimGeometry::kLeft;
01404         }
01405       }
01406       if( py>0.0 && y0<=y_left ){
01407         if( y0>y_right
01408          || useFarEdge ){
01409           x = x_left;  
01410           y = y_left;  
01411           z = z_left;
01412           region = WCSimGeometry::kLeft;
01413         }
01414       }    
01415     }          
01416 
01417     if( z_right>=-0.5*dZ && z_right<=+0.5*dZ
01418      && x_right>=-0.5*dX && x_right<=+0.5*dX ){
01419       if( py>0.0 && y0<y_right ){
01420         if( !useFarEdge ){
01421           x = x_right;  
01422           y = y_right;  
01423           z = z_right;
01424           region = WCSimGeometry::kRight;
01425         }
01426       }
01427       if( py<0.0 && y0>=y_right ){
01428         if( y0<y_left
01429          || useFarEdge ){
01430           x = x_right;  
01431           y = y_right;  
01432           z = z_right;
01433           region = WCSimGeometry::kRight;
01434         }
01435       }
01436     }      
01437 
01438     // return projection
01439     // =================
01440     if( region > WCSimGeometry::kUnknown ){
01441  
01442       xproj = x;
01443       yproj = y;
01444       zproj = z;
01445       regionproj = region;
01446 
01447       return;
01448     }
01449   }
01450 
01451   return;
01452 }
01453  
01454 void WCSimGeometry::ProjectToEdge(Bool_t useFarEdge, Double_t x0, Double_t y0, Double_t z0, Double_t px, Double_t py, Double_t pz, Double_t& xproj, Double_t& yproj, Double_t&zproj, Int_t& regionproj)
01455 {
01456   // default locations
01457   // =================
01458   xproj = -99999.9;
01459   yproj = -99999.9;
01460   zproj = -99999.9;
01461   regionproj = WCSimGeometry::kUnknown;
01462 
01463   Double_t xNear = -99999.9;
01464   Double_t yNear = -99999.9;
01465   Double_t zNear = -99999.9;
01466   Int_t regionNear = WCSimGeometry::kUnknown;
01467 
01468   Double_t xFar = -99999.9;
01469   Double_t yFar = -99999.9;
01470   Double_t zFar = -99999.9;
01471   Int_t regionFar = WCSimGeometry::kUnknown;
01472 
01473 
01474   // CYLINDRICAL GEOMETRY
01475   // ====================
01476   if( fGeoType==WCSimGeometry::kCylinder ){  
01477 
01478     Double_t r = fCylRadius;
01479     Double_t L = fCylLength;
01480 
01481     Bool_t foundProjectionXY = 0;
01482     Bool_t foundProjectionZ = 0;
01483 
01484     Double_t t1 = 0.0;
01485     Double_t x1 = 0.0;
01486     Double_t y1 = 0.0;
01487     Double_t z1 = 0.0;
01488     Int_t region1 = -1;
01489 
01490     Double_t t2 = 0.0;  
01491     Double_t x2 = 0.0;
01492     Double_t y2 = 0.0;
01493     Double_t z2 = 0.0;
01494     Int_t region2 = -1;
01495 
01496     Double_t rSq = r*r;
01497     Double_t r0r0 = x0*x0 + y0*y0;
01498     Double_t r0p = x0*px + y0*py;
01499     Double_t pSq = px*px+py*py;
01500     
01501     // calculate intersection in XY
01502     if( pSq>0.0 ){
01503       if( r0p*r0p - pSq*(r0r0-rSq)>0.0 ){
01504         t1 = ( -r0p - sqrt(r0p*r0p-pSq*(r0r0-rSq)) ) / pSq;
01505         t2 = ( -r0p + sqrt(r0p*r0p-pSq*(r0r0-rSq)) ) / pSq;
01506         foundProjectionXY = 1;
01507       }
01508     }
01509 
01510     // propagation along z-axis
01511     else if( r0r0<=rSq ){
01512 
01513       if( pz>0 ){
01514         t1 = -L/2.0 - z0;
01515         t2 = +L/2.0 - z0;
01516       }
01517       else{
01518         t1 = -L/2.0 + z0;
01519         t2 = +L/2.0 + z0;
01520       }
01521       foundProjectionXY = 1;
01522     }
01523     
01524     // found intersection in XY
01525     if( foundProjectionXY ){
01526 
01527       z1 = z0 + t1*pz;
01528       z2 = z0 + t2*pz;
01529 
01530       if( ( z1>=-L/2.0 && z2<=+L/2.0 )
01531        || ( z2>=-L/2.0 && z1<=+L/2.0 ) ){
01532         foundProjectionZ = 1;
01533       }
01534     }
01535 
01536     // found intersection in Z
01537     if( foundProjectionZ ){
01538 
01539       // first intersection
01540       if( z1>-L/2.0 && z1<+L/2.0 ){
01541         region1 = WCSimGeometry::kSide;
01542       }
01543       if( z1>=+L/2.0 ){
01544         region1 = WCSimGeometry::kTop;
01545         if( z1>+L/2.0 ){
01546           z1 = +L/2.0; 
01547           t1 = (+L/2.0-z0)/pz;
01548         }
01549       }
01550       if( z1<=-L/2.0 ){
01551         region1 = WCSimGeometry::kBottom;
01552         if( z1<-L/2.0 ){
01553           z1 = -L/2.0; 
01554           t1 = (-L/2.0-z0)/pz;
01555         }
01556       }
01557 
01558       x1 = x0 + t1*px;
01559       y1 = y0 + t1*py;
01560 
01561       // second intersection
01562       if( z2>-L/2.0 && z2<+L/2.0 ){
01563         region2 = WCSimGeometry::kSide;
01564       }
01565       if( z2>=+L/2.0 ){
01566         region2 = WCSimGeometry::kTop;
01567         if( z2>+L/2.0 ){
01568           z2 = +L/2.0; 
01569           t2 = (+L/2.0-z0)/pz;
01570         }
01571       }
01572       if( z2<=-L/2.0 ){
01573         region2 = WCSimGeometry::kBottom;
01574         if( z2<-L/2.0 ){
01575           z2 = -L/2.0; 
01576           t2 = (-L/2.0-z0)/pz;
01577         }
01578       }
01579 
01580       x2 = x0 + t2*px;
01581       y2 = y0 + t2*py;
01582 
01583       // near/far projection
01584       if( t1>=0 ){
01585         xNear = x1;
01586         yNear = y1;
01587         zNear = z1;
01588         regionNear = region1;
01589 
01590         xFar = x2;
01591         yFar = y2;
01592         zFar = z2;
01593         regionFar = region2;
01594       }
01595       else if( t2>0 ){
01596         xNear = x2;
01597         yNear = y2;
01598         zNear = z2;
01599         regionNear = region2;
01600 
01601         xFar = x2;
01602         yFar = y2;
01603         zFar = z2;
01604         regionFar = region2;
01605       }
01606     }
01607  
01608   }
01609 
01610   // MAILBOX GEOMETRY
01611   // ================
01612   if( fGeoType==WCSimGeometry::kMailBox ){
01613 
01614     Double_t dX = fMailBoxX;
01615     Double_t dY = fMailBoxY;
01616     Double_t dZ = fMailBoxZ;
01617 
01618     Bool_t foundProjection = 0;
01619 
01620     Double_t t = 0.0;
01621     Double_t x = 0.0;
01622     Double_t y = 0.0;
01623     Double_t z = 0.0;
01624 
01625     Double_t t1 = 0.0;
01626     Double_t x1 = 0.0;
01627     Double_t y1 = 0.0;
01628     Double_t z1 = 0.0;
01629     Int_t region1 = -1;
01630 
01631     Double_t t2 = 0.0;  
01632     Double_t x2 = 0.0;
01633     Double_t y2 = 0.0;
01634     Double_t z2 = 0.0;
01635     Int_t region2 = -1;
01636 
01637     if( pz!=0 ){
01638 
01639       // top face (+dZ/2.0)
01640       t = (+dZ/2.0-z0)/pz;
01641       x = x0 + t*px;
01642       y = y0 + t*py;
01643       if( x>=-dX/2.0 && x<=+dX/2.0 
01644        && y>=-dY/2.0 && y<=+dY/2.0 ){
01645         if( pz>0.0 ){
01646           t2 = t;
01647           region2 = WCSimGeometry::kTop;
01648           foundProjection = 1;
01649         }
01650         else{
01651           t1 = t;
01652           region1 = WCSimGeometry::kTop;
01653           foundProjection = 1;
01654         }
01655       }
01656 
01657       // bottom face (-dZ/2.0)
01658       t = (-dZ/2.0-z0)/pz;
01659       x = x0 + t*px;
01660       y = y0 + t*py;
01661       if( x>=-dX/2.0 && x<=+dX/2.0 
01662        && y>=-dY/2.0 && y<=+dY/2.0 ){
01663         if( pz<0.0 ){
01664           t2 = t;
01665           region2 = WCSimGeometry::kBottom;
01666           foundProjection = 1;
01667         }
01668         else{
01669           t1 = t;
01670           region1 = WCSimGeometry::kBottom;
01671           foundProjection = 1;
01672         }
01673       }
01674     }
01675 
01676     if( px!=0 ){
01677 
01678       // front face (+dX/2.0)
01679       t = (+dX/2.0-x0)/px;
01680       y = y0 + t*py;
01681       z = z0 + t*pz;
01682       if( y>=-dY/2.0 && y<=+dY/2.0 
01683        && z>=-dZ/2.0 && z<=+dZ/2.0 ){
01684         if( px>0.0 ){
01685           t2 = t;
01686           region2 = WCSimGeometry::kFront;
01687           foundProjection = 1;
01688         }
01689         else{
01690           t1 = t;
01691           region1 = WCSimGeometry::kFront;
01692           foundProjection = 1;
01693         }
01694       }
01695 
01696       // back face (-dX/2.0)
01697       t = (-dX/2.0-x0)/px;
01698       y = y0 + t*py;
01699       z = z0 + t*pz;
01700       if( y>=-dY/2.0 && y<=+dY/2.0 
01701        && z>=-dZ/2.0 && z<=+dZ/2.0 ){
01702         if( px<0.0 ){
01703           t2 = t;
01704           region2 = WCSimGeometry::kBack;
01705           foundProjection = 1;
01706         }
01707         else{
01708           t1 = t;
01709           region1 = WCSimGeometry::kBack;
01710           foundProjection = 1;
01711         }
01712       }
01713     }
01714 
01715     if( py!=0 ){
01716 
01717       // left face (+dY/2.0)
01718       t = (+dY/2.0-y0)/py;
01719       z = z0 + t*pz;
01720       x = x0 + t*px;
01721       if( z>=-dZ/2.0 && z<=+dZ/2.0 
01722        && x>=-dX/2.0 && x<=+dX/2.0 ){
01723         if( py>0.0 ){
01724           t2 = t;
01725           region2 = WCSimGeometry::kLeft;
01726           foundProjection = 1;
01727         }
01728         else{
01729           t1 = t;
01730           region1 = WCSimGeometry::kLeft;
01731           foundProjection = 1;
01732         }
01733       }
01734     
01735       // right face (-dY/2.0)
01736       t = (-dY/2.0-y0)/py;
01737       z = z0 + t*pz;
01738       x = x0 + t*px;
01739       if( z>=-dZ/2.0 && z<=+dZ/2.0 
01740        && x>=-dX/2.0 && x<=+dX/2.0 ){
01741         if( py<0.0 ){
01742           t2 = t;
01743           region2 = WCSimGeometry::kRight;
01744           foundProjection = 1;
01745         }
01746         else{
01747           t1 = t;
01748           region1 = WCSimGeometry::kRight;
01749           foundProjection = 1;
01750         }
01751       }
01752     }
01753 
01754     // found projection
01755     if( foundProjection = 1 ){
01756 
01757       x1 = x0 + t1*px;
01758       y1 = y0 + t1*py;
01759       z1 = z0 + t1*pz;
01760 
01761       x2 = x0 + t2*px;
01762       y2 = y0 + t2*py;
01763       z2 = z0 + t2*pz;
01764 
01765       // near/far projection
01766       if( t1>=0 ){
01767         xNear = x1;
01768         yNear = y1;
01769         zNear = z1;
01770         regionNear = region1;
01771 
01772         xFar = x2;
01773         yFar = y2;
01774         zFar = z2;
01775         regionFar = region2;
01776       }
01777       else if( t2>0 ){
01778         xNear = x2;
01779         yNear = y2;
01780         zNear = z2;
01781         regionNear = region2;
01782 
01783         xFar = x2;
01784         yFar = y2;
01785         zFar = z2;
01786         regionFar = region2;
01787       }
01788     }
01789   }
01790 
01791   // set projections
01792   // ===============
01793   if( useFarEdge ){
01794     xproj = xFar;
01795     yproj = yFar;
01796     zproj = zFar;
01797     regionproj = regionFar;
01798   }
01799   else{
01800     xproj = xNear;
01801     yproj = yNear;
01802     zproj = zNear;
01803     regionproj = regionNear;
01804   }
01805 
01806   return;
01807 }
01808 
01809 void WCSimGeometry::XYZtoUV(Int_t region, Double_t x, Double_t y, Double_t z, Double_t& u, Double_t& v)
01810 {
01811   // default locations
01812   u = -99999.9;
01813   v = -99999.9;
01814 
01815   // Note:
01816   // ====
01817   //  assume that beam points along -ve x direction
01818 
01819   // CYLINDRICAL GEOMETRY
01820   // ====================
01821   if( fGeoType==WCSimGeometry::kCylinder ){
01822 
01823     // PMTs on top face
01824     if( region==WCSimGeometry::kTop ){
01825       u = y;
01826       v = +0.5*fCylLength+fCylRadius+x;
01827     }
01828 
01829     // PMTs on side
01830     if( region==WCSimGeometry::kSide ){
01831       Double_t theta = 0.0;
01832       if( x!=0.0 ) theta = TMath::ATan(-y/x);
01833       if( x>=0 ){
01834         if( y>0.0 ) theta += TMath::Pi();
01835         if( y<0.0 ) theta -= TMath::Pi();
01836       }
01837       u = fCylRadius*theta;
01838       v = z;
01839     }
01840 
01841     // PMTs on bottom face
01842     if( region==WCSimGeometry::kBottom ){
01843       u = y;
01844       v = -0.5*fCylLength-fCylRadius-x;
01845     }
01846   }
01847 
01848   // MAILBOX GEOMETRY
01849   // ================
01850   //
01851   // Default Coordinate System:
01852   //
01853   //               +ve z(Top)
01854   //                   ^   
01855   //                   |  / +ve y(Left)
01856   //                   | /
01857   //    -ve x(Back) -- --> +ve x(Front)          
01858   //                 / |
01859   //   -ve y(Right) /  |
01860   //                  -ve z(Bottom)
01861   //
01862   // Rolled out as follows:
01863   // 
01864   //              -----
01865   //             |  T  |
01866   //    --- ----- ----- ----- ---
01867   //   | F |  R  |  B  |  L  | F |
01868   //    --- ----- ----- ----- ---
01869   //             |  B  |
01870   //              -----
01871   //
01872 
01873   if( fGeoType==WCSimGeometry::kMailBox ){
01874 
01875     // PMTs on top face
01876     if( region==WCSimGeometry::kTop ){
01877       u = y;
01878       v = +0.5*fMailBoxZ+0.5*fMailBoxX+x;
01879     }
01880 
01881     // PMTs on front face
01882     if( region==WCSimGeometry::kFront ){      
01883       if( y>0.0 ) u = +fMailBoxX+fMailBoxY-y;
01884       if( y<0.0 ) u = -fMailBoxX-fMailBoxY-y;
01885       v = z;
01886     }
01887 
01888     // PMTs on left face
01889     if( region==WCSimGeometry::kLeft ){      
01890       u = +0.5*fMailBoxX+0.5*fMailBoxY+x;
01891       v = z;
01892     }
01893 
01894     // PMTs on right face
01895     if( region==WCSimGeometry::kRight ){
01896       u = -0.5*fMailBoxX-0.5*fMailBoxY-x;
01897       v = z;
01898     }
01899 
01900     // PMTs on back face
01901     if( region==WCSimGeometry::kBack ){
01902       u = y;
01903       v = z;
01904     }    
01905 
01906     // PMTs on bottom face
01907     if( region==WCSimGeometry::kBottom ){
01908       u = y;
01909       v = -0.5*fMailBoxZ-0.5*fMailBoxX-x;
01910     }
01911 
01912   }
01913 
01914   return;
01915 }
01916 
01917 void WCSimGeometry::FindCircle(Double_t x0, Double_t y0, Double_t z0, Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2, Double_t& rx, Double_t& ry, Double_t& rz, Double_t& nx, Double_t& ny, Double_t& nz, Double_t& r)
01918 {
01919   Double_t centre[3] = {0.0,0.0,0.0};
01920   Double_t normal[3] = {0.0,0.0,0.0};
01921   Double_t radius = 0.0;
01922 
01923   Double_t p01[3] = {x0-x1, y0-y1, z0-z1};
01924   Double_t p12[3] = {x1-x2, y1-y2, z1-z2};
01925   Double_t p20[3] = {x2-x0, y2-y0, z2-z0};
01926 
01927   Double_t n[3] = {p01[1]*p12[2]-p01[2]*p12[1],
01928                    p01[2]*p12[0]-p01[0]*p12[2],
01929                    p01[0]*p12[1]-p01[1]*p12[0]};
01930 
01931   Double_t D2 = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
01932 
01933   if( D2>0.0 ){
01934     Double_t a0 = -(p12[0]*p12[0]+p12[1]*p12[1]+p12[2]*p12[2])*(p01[0]*p20[0]+p01[1]*p20[1]+p01[2]*p20[2])/(2.0*D2);
01935     Double_t a1 = -(p20[0]*p20[0]+p20[1]*p20[1]+p20[2]*p20[2])*(p12[0]*p01[0]+p12[1]*p01[1]+p12[2]*p01[2])/(2.0*D2);
01936     Double_t a2 = -(p01[0]*p01[0]+p01[1]*p01[1]+p01[2]*p01[2])*(p20[0]*p12[0]+p20[1]*p12[1]+p20[2]*p12[2])/(2.0*D2);
01937     Double_t D = sqrt(D2);
01938 
01939     centre[0] = a0*x0 + a1*x1 + a2*x2;
01940     centre[1] = a0*y0 + a1*y1 + a2*y2;
01941     centre[2] = a0*z0 + a1*z1 + a2*z2;
01942 
01943     radius = sqrt( (p01[0]*p01[0]+p01[1]*p01[1]+p01[2]*p01[2])
01944                   *(p12[0]*p12[0]+p12[1]*p12[1]+p12[2]*p12[2])
01945                   *(p20[0]*p20[0]+p20[1]*p20[1]+p20[2]*p20[2]))/(2.0*D);
01946 
01947     if( n[0]*centre[0]
01948       + n[1]*centre[1]
01949       + n[2]*centre[2]>0.0 ){
01950       normal[0] = +n[0]/D;
01951       normal[1] = +n[1]/D;
01952       normal[2] = +n[2]/D;
01953     }
01954     else{
01955       normal[0] = -n[0]/D;
01956       normal[1] = -n[1]/D;
01957       normal[2] = -n[2]/D;
01958     }
01959   }
01960 
01961   rx = centre[0];
01962   ry = centre[1];
01963   rz = centre[2];
01964 
01965   nx = normal[0];
01966   ny = normal[1];
01967   nz = normal[2];
01968 
01969   r = radius;
01970 
01971   return;
01972 }
01973 
01974 void WCSimGeometry::FindCircleOld( Double_t xp, Double_t yp, Double_t zp, Double_t x0, Double_t y0, Double_t z0, Double_t angle_degrees, Double_t omega_degrees, Double_t& rx, Double_t& ry, Double_t& rz, Double_t& nx, Double_t& ny, Double_t& nz, Double_t& r ) 
01975 {
01976   // inputs
01977   // ======
01978   Double_t phi = 0.0;
01979   Double_t theta = 0.0;
01980   Double_t angle = (TMath::Pi()/180.0)*angle_degrees;
01981   Double_t omega = (TMath::Pi()/180.0)*omega_degrees;
01982 
01983   Double_t x = xp;
01984   Double_t y = yp;
01985   Double_t z = zp;
01986 
01987   Double_t xtemp = 0.0;
01988   Double_t ytemp = 0.0;
01989   Double_t ztemp = 0.0;
01990 
01991   Double_t dx = 0.0;
01992   Double_t dy = 0.0;
01993   Double_t dz = 0.0;
01994   Double_t ds = 0.0;
01995 
01996   Double_t radius = 0.0;
01997 
01998   // subtract vertex
01999   // ===============
02000   x -= x0;
02001   y -= y0;
02002   z -= z0;
02003 
02004   // forward rotation (x,y,z)->(x,0,z)
02005   // =================================
02006   phi = 0.0;
02007 
02008   if( x!=0.0 ){
02009     phi = atan(y/x);
02010   }
02011   if( x<=0.0 ){
02012     if( y>0.0 ) phi += TMath::Pi();
02013     if( y<0.0 ) phi -= TMath::Pi();
02014   }
02015 
02016   xtemp = x;
02017   ytemp = y;
02018   ztemp = z;  
02019 
02020   x = + xtemp*cos(phi) + ytemp*sin(phi);
02021   y = - xtemp*sin(phi) + ytemp*cos(phi);
02022   z = + ztemp;
02023 
02024   // forward rotation (x,0,z)->(0,0,z)
02025   // =================================
02026   theta = 0.0;
02027 
02028   if( z!=0.0 ){
02029     theta = atan(x/z);
02030   }
02031   if( z<=0.0 ){
02032     if( x>0.0 ) theta += TMath::Pi();
02033     if( x<0.0 ) theta -= TMath::Pi();
02034   }
02035 
02036   xtemp = x;
02037   ytemp = y;
02038   ztemp = z;  
02039  
02040   x = + xtemp*cos(theta) - ztemp*sin(theta);
02041   y = + ytemp;
02042   z = + xtemp*sin(theta) + ztemp*cos(theta);
02043 
02044   // apply rotation: cone angle + azimuthal angle
02045   // ============================================  
02046   xtemp = x;
02047   ytemp = y;
02048   ztemp = z;  
02049 
02050   x = ztemp*sin(angle)*cos(omega);
02051   y = ztemp*sin(angle)*sin(omega);
02052   z = ztemp*cos(angle);
02053 
02054   radius = ztemp*sin(angle);
02055 
02056   // backward rotation (0,0,z)->(x,0,z)
02057   // ==================================
02058   xtemp = x;
02059   ytemp = y;
02060   ztemp = z;  
02061  
02062   x = + xtemp*cos(-theta) - ztemp*sin(-theta);
02063   y = + ytemp;
02064   z = + xtemp*sin(-theta) + ztemp*cos(-theta);
02065 
02066   // backward rotation (x,0,z)->(x,y,z)
02067   // =================================
02068   xtemp = x;
02069   ytemp = y;
02070   ztemp = z;  
02071 
02072   x = + xtemp*cos(-phi) + ytemp*sin(-phi);
02073   y = - xtemp*sin(-phi) + ytemp*cos(-phi);
02074   z = + ztemp;
02075 
02076   // add on vertex
02077   // =============
02078   x += x0;
02079   y += y0;
02080   z += z0;
02081 
02082   // return coordinates
02083   // ==================
02084   rx = x;
02085   ry = y;
02086   rz = z;
02087 
02088   // return normal 
02089   // =============
02090   dx = rx - x0;
02091   dy = ry - y0;
02092   dz = rz - z0;
02093   ds = sqrt(dx*dx+dy*dy+dz*dz);
02094 
02095   nx = dx/ds;
02096   ny = dy/ds;
02097   nz = dz/ds;
02098 
02099   // return radius
02100   // =============
02101   r = radius;
02102 
02103   return;
02104 }
02105 
02106 void WCSimGeometry::FindCircle( Double_t xp, Double_t yp, Double_t zp, Double_t x0, Double_t y0, Double_t z0, Double_t angle_degrees, Double_t omega_degrees, Double_t& rx, Double_t& ry, Double_t& rz, Double_t& nx, Double_t& ny, Double_t& nz, Double_t& r ) 
02107 {
02108   // default
02109   // =======
02110   rx = xp;
02111   ry = yp;
02112   rz = zp;
02113   
02114   nx = 0.0;
02115   ny = 0.0;
02116   nz = 0.0;
02117 
02118   r = 0.0;
02119 
02120   // inputs
02121   // ======
02122   Double_t angle = (TMath::Pi()/180.0)*angle_degrees; // radians
02123   Double_t omega = (TMath::Pi()/180.0)*omega_degrees; // radians
02124 
02125   Double_t dx = xp-x0;
02126   Double_t dy = yp-y0;
02127   Double_t dz = zp-z0;
02128   Double_t ds = sqrt(dx*dx+dy*dy+dz*dz);
02129 
02130   Double_t px = 0.0;
02131   Double_t py = 0.0;
02132   Double_t pz = 0.0;
02133 
02134   if( ds>0.0 ){
02135     px = dx/ds;
02136     py = dy/ds;
02137     pz = dz/ds;
02138   }
02139   else{
02140     return;
02141   }
02142 
02143   // do the rotations
02144   // ================
02145   TVector3 myVtx(x0,y0,z0);
02146   TVector3 myDir(px,py,pz);
02147 
02148   TVector3 initDir = myDir;
02149   TVector3 orthDir = myDir.Orthogonal();
02150   myDir.Rotate(angle,orthDir);
02151   myDir.Rotate(omega,initDir);
02152  
02153   // outputs
02154   // =======
02155   rx = x0 + ds*myDir.x();
02156   ry = y0 + ds*myDir.y();
02157   rz = z0 + ds*myDir.z();
02158   
02159   nx = myDir.x();
02160   ny = myDir.y();
02161   nz = myDir.z();
02162 
02163   r = ds*sin(angle);
02164 
02165   return;
02166 }
02167 
02168 void WCSimGeometry::FindVertex(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, Double_t& vxm, Double_t& vym, Double_t& vzm, Double_t& vtm, Double_t& vxp, Double_t& vyp, Double_t& vzp, Double_t& vtp)
02169 {
02170   // switch off error messages
02171   // =========================
02172   // suppress messages like: 
02173   //  Error in <TDecompLU::DecomposeLUCrout>: matrix is singular
02174   Int_t myIgnoreLevel = gErrorIgnoreLevel;
02175   gErrorIgnoreLevel = kFatal;
02176   
02177   // default vertex
02178   // ==============
02179   vxm = -99999.9; 
02180   vym = -99999.9;
02181   vzm = -99999.9;
02182   vtm = -99999.9;
02183 
02184   vxp = -99999.9; 
02185   vyp = -99999.9;
02186   vzp = -99999.9;
02187   vtp = -99999.9;  
02188 
02189   // speed of light in water
02190   // =======================
02191   Double_t c = 29.98/1.33; // cm/ns [water]
02192 
02193   // causality checks
02194   // ================
02195   if( (x1-x0)*(x1-x0) + (y1-y0)*(y1-y0) + (z1-z0)*(z1-z0) >= c*c*(t1-t0)*(t1-t0)
02196    && (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1) >= c*c*(t2-t1)*(t2-t1)
02197    && (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) + (z3-z2)*(z3-z2) >= c*c*(t3-t2)*(t3-t2)
02198    && (x2-x0)*(x2-x0) + (y2-y0)*(y2-y0) + (z2-z0)*(z2-z0) >= c*c*(t2-t0)*(t2-t0)
02199    && (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1) + (z3-z1)*(z3-z1) >= c*c*(t3-t1)*(t3-t1)
02200    && (x3-x0)*(x3-x0) + (y3-y0)*(y3-y0) + (z3-z0)*(z3-z0) >= c*c*(t3-t0)*(t3-t0) ){
02201 
02202     // [Note: for causality, require that |x_{i}-x_{j}| >= c*|t_{i}-t_{j}|
02203     //        for each pair of points]
02204 
02205     Double_t dx1 = x1-x0;  Double_t dy1 = y1-y0;  Double_t dz1 = z1-z0;  Double_t dt1 = c*(t1-t0);
02206     Double_t dx2 = x2-x0;  Double_t dy2 = y2-y0;  Double_t dz2 = z2-z0;  Double_t dt2 = c*(t2-t0);
02207     Double_t dx3 = x3-x0;  Double_t dy3 = y3-y0;  Double_t dz3 = z3-z0;  Double_t dt3 = c*(t3-t0);
02208 
02209     Double_t epsilon = 1.0e-7;
02210 
02211     // check that points don't all lie in a plane
02212     if( !( fabs(dx1)<epsilon && fabs(dx2)<epsilon && fabs(dx3)<epsilon )
02213      && !( fabs(dy1)<epsilon && fabs(dy2)<epsilon && fabs(dy3)<epsilon )
02214      && !( fabs(dz1)<epsilon && fabs(dz2)<epsilon && fabs(dz3)<epsilon ) 
02215      && !( fabs(dx1)<epsilon && fabs(dy1)<epsilon && fabs(dz1)<epsilon )
02216      && !( fabs(dx2)<epsilon && fabs(dy2)<epsilon && fabs(dz2)<epsilon )
02217      && !( fabs(dx3)<epsilon && fabs(dy3)<epsilon && fabs(dz3)<epsilon ) ){
02218 
02219       // [Note: this is a problem for detectors with flat faces!]
02220 
02221       Double_t Mdata[9] = { dx1, dy1, dz1,
02222                             dx2, dy2, dz2,
02223                             dx3, dy3, dz3 };
02224 
02225       Double_t Qdata[3] = { 0.5*( dx1*dx1 + dy1*dy1 + dz1*dz1 - dt1*dt1 ),
02226                             0.5*( dx2*dx2 + dy2*dy2 + dz2*dz2 - dt2*dt2 ),
02227                             0.5*( dx3*dx3 + dy3*dy3 + dz3*dz3 - dt3*dt3 ) };
02228                       
02229       Double_t Tdata[3] = { dt1,
02230                             dt2,
02231                             dt3 };
02232      
02233       TMatrixD M(3,3,Mdata);
02234       TMatrixD Q(3,1,Qdata);
02235       TMatrixD T(3,1,Tdata);
02236 
02237       if( M.Determinant() != 0.0 ){
02238 
02239         TMatrixD A(3,1);
02240         TMatrixD B(3,1);
02241 
02242         M.Invert();
02243         A.Mult(M,T);
02244         B.Mult(M,Q);
02245     
02246         Double_t ax = A(0,0);
02247         Double_t ay = A(1,0);
02248         Double_t az = A(2,0);
02249 
02250         Double_t bx = B(0,0);
02251         Double_t by = B(1,0);
02252         Double_t bz = B(2,0);
02253 
02254         Double_t ab = ax*bx + ay*by + az*bz;
02255         Double_t a2 = ax*ax + ay*ay + az*az;
02256         Double_t b2 = bx*bx + by*by + bz*bz;
02257 
02258         Double_t qa = a2-1.0;
02259         Double_t qb = 2.0*ab;
02260         Double_t qc = b2;
02261 
02262         // check for solutions
02263         if( qb*qb-4.0*qa*qc>0.0 ){
02264 
02265           // The common vertex is given by a quadratic equation, which has two solutions.
02266           // Typically, one solution corresponds to photons travelling forwards in time,
02267           // and the other solution corresponds to photons travelling backwards in time.
02268           // However, sometimes there appear to be two valid solutions.
02269 
02270           Double_t ctm = ( -qb - sqrt(qb*qb-4.0*qa*qc) ) / ( 2.0*qa );
02271           Double_t ctp = ( -qb + sqrt(qb*qb-4.0*qa*qc) ) / ( 2.0*qa );
02272 
02273           Double_t tm = t0 + ctm/c;
02274           Double_t xm = x0 + ctm*ax + bx;
02275           Double_t ym = y0 + ctm*ay + by;
02276           Double_t zm = z0 + ctm*az + bz;
02277           Bool_t foundVertexM = 0;
02278 
02279           if( tm<t0 && tm<t1 
02280            && tm<t2 && tm<t3 ){
02281             vxm = xm;
02282             vym = ym;
02283             vzm = zm;
02284             vtm = tm;
02285             foundVertexM = 1;
02286           }
02287 
02288           Double_t tp = t0 + ctp/c;
02289           Double_t xp = x0 + ctp*ax + bx;
02290           Double_t yp = y0 + ctp*ay + by;
02291           Double_t zp = z0 + ctp*az + bz;
02292           Bool_t foundVertexP = 0;
02293 
02294           if( tp<t0 && tp<t1 
02295            && tp<t2 && tp<t3 ){
02296             vxp = xp;
02297             vyp = yp;
02298             vzp = zp;
02299             vtp = tp;
02300             foundVertexP = 1;
02301           }
02302 
02303           // std::cout << "  Vertex: [a=" << qa << ",b=" << qb << ",c=" << qc << "] [ctm=" << ctm << ",ctp=" << ctp << "] " << std::endl
02304           //           << "   Low T:  (x,y,z,t)=(" << xm << "," << ym << "," << zm << "," << tm << ") [foundVtx=" << foundVertexM << "] " << std::endl
02305           //           << "   High T: (x,y,z,t)=(" << xp << "," << yp << "," << zp << "," << tp << ") [foundVtx=" << foundVertexP << "] " << std::endl;
02306           
02307         }
02308       }
02309     }
02310   }
02311 
02312   // switch on error messages
02313   // ========================
02314   gErrorIgnoreLevel = myIgnoreLevel;
02315 
02316   return;
02317 }
02318 
02319 
02320 
02321 void WCSimGeometry::DistanceToIntersectLine(Double_t x0, Double_t y0, Double_t z0, Double_t vx, Double_t vy, Double_t vz, Double_t ex, Double_t ey, Double_t ez, Double_t& x, Double_t& y, Double_t& z, Double_t& L)
02322 {
02323   Double_t dx = 0.0;
02324   Double_t dy = 0.0;
02325   Double_t dz = 0.0;
02326   Double_t ds = 0.0;
02327   
02328   Double_t px = 0.0;
02329   Double_t py = 0.0;
02330   Double_t pz = 0.0;
02331   Double_t dsTrack = 0.0;
02332 
02333   Double_t qx = 0.0;
02334   Double_t qy = 0.0;
02335   Double_t qz = 0.0;
02336   Double_t dsPmt = 0.0;
02337 
02338   // vector from vertex to end
02339   dx = ex - vx;
02340   dy = ey - vy;
02341   dz = ez - vz;
02342   ds = sqrt( dx*dx+dy*dy+dz*dz );
02343 
02344   if( ds>0.0 ){
02345     px = dx/ds;
02346     py = dy/ds;
02347     pz = dz/ds;
02348     dsTrack = ds;
02349   }
02350 
02351   // vector from vertex to pmt
02352   dx = x0 - vx;
02353   dy = y0 - vy;
02354   dz = z0 - vz;
02355   ds = sqrt( dx*dx+dy*dy+dz*dz );
02356 
02357   if( ds>0.0 ){
02358     qx = dx/ds;
02359     qy = dy/ds;
02360     qz = dz/ds;
02361     dsPmt = ds;
02362   }
02363      
02364   // Cherenkov geometry
02365   // ==================
02366   Double_t cosphi = px*qx+py*qy+pz*qz;
02367   Double_t phi = acos(cosphi); 
02368   Double_t theta = WCSimParameters::CherenkovAngle()*TMath::Pi()/180.0;    
02369 
02370   Double_t Ltrack = 0.0;
02371   Double_t Lphoton = 0.0;
02372     
02373   if( phi<theta ){
02374     Ltrack  = +dsPmt*sin(theta-phi)/sin(theta);
02375     Lphoton = +dsPmt*sin(phi)/sin(theta);
02376   }
02377   else{
02378     Ltrack  = -dsPmt*sin(phi-theta)/sin(phi);
02379     Lphoton = -dsPmt*sin(phi)/sin(theta);
02380   }
02381  
02382   // return intersection
02383   // ===================
02384   x = vx + px*Ltrack;
02385   y = vy + py*Ltrack;
02386   z = vz + pz*Ltrack;
02387   L = Lphoton;
02388 
02389   return;
02390 }
02391 
02392 
02393 Double_t WCSimGeometry::DistanceToIntersectLine(Double_t x0, Double_t y0, Double_t z0, Double_t sx, Double_t sy, Double_t sz, Double_t ex, Double_t ey, Double_t ez, Double_t& x, Double_t& y, Double_t& z)
02394 {
02395   Double_t L = -1.0;
02396 
02397   WCSimGeometry::DistanceToIntersectLine(x0,y0,z0,
02398                                          sx,sy,sz,
02399                                          ex,ey,ez,
02400                                          x,y,z,L);
02401   return L;
02402 }
02403 
02404 Double_t WCSimGeometry::DistanceToIntersectLine(Double_t* pos, Double_t* start, Double_t* end, Double_t* intersection)
02405 {
02406   Double_t x = -999.9;
02407   Double_t y = -999.9;
02408   Double_t z = -999.9; 
02409   Double_t L = -1.0;
02410 
02411   WCSimGeometry::DistanceToIntersectLine(pos[0],   pos[1],   pos[2],
02412                                          start[0], start[1], start[2],
02413                                          end[0],   end[1],   end[2],
02414                                          x,y,z,L);
02415 
02416   intersection[0]=x; intersection[1]=y; intersection[2]=z;  
02417   
02418   return L;
02419 }
02420 
02421 //
02422 //
02423 // // Note: have replaced the below code with an implementation
02424 // //       using the sine rule, which I think is quicker
02425 //
02426 // Double_t WCSimGeometry::DistanceToIntersectLine(Double_t x0, Double_t y0, Double_t z0, 
02427 //                                    Double_t sx, Double_t sy, Double_t sz,
02428 //                                    Double_t ex, Double_t ey, Double_t ez, 
02429 //                                    Double_t x, Double_t y, Double_t z)
02430 // {
02431 //   Double_t vstart[3], dist[3], pos[3], start[3], end[3], intersection[3]; /// dist = vector(end - start), vstart = vector(pos - start)
02432 //   pos[0] = x0; pos[1] = y0; pos[2] = z0;
02433 //   start[0] = sx; start[1] = sy; start[2] = sz;
02434 //   end[0] = ex; end[1] = ey; end[2] = ez; 
02435 //   intersection[0] = x; intersection[1] = y; intersection[2] = z; 
02436 //   Double_t ctc    = WCSimParameters::CosThetaC;
02437 //   Double_t length = WCSimMathUtil::Magnitude(start,end);
02438 //   Double_t arg_a  = pow(length,4)*(1-WCSimMathUtil::sqr(ctc));
02439 //
02440 //  
02441 //   WCSimMathUtil::SubtractVec(pos,start,vstart);
02442 //   WCSimMathUtil::SubtractVec(end,start,dist);
02443 //  
02444 //   Double_t arg_b  = -2*WCSimMathUtil::sqr(length)*WCSimMathUtil::DotProduct(vstart,dist)*(1-WCSimMathUtil::sqr(ctc));
02445 //   Double_t arg_c = WCSimMathUtil::sqr(WCSimMathUtil::DotProduct(vstart,dist))-WCSimMathUtil::sqr(WCSimMathUtil::Magnitude(pos,start))*WCSimMathUtil::sqr(length)*WCSimMathUtil::sqr(ctc);
02446 //
02447 //   Double_t firstSol;
02448 //   Double_t secondSol;
02449 //   Double_t det = WCSimMathUtil::Determinant(arg_a, arg_b, arg_c);
02450 //
02451 //   if (det >= 0) {
02452 //     WCSimMathUtil::RealRoots(arg_a, arg_b, sqrt(det), firstSol, secondSol);  // roots are real
02453 //     Double_t newPos[3], newvstart[3], newdist[3];
02454 //     Double_t theta1, theta2;
02455 //     Double_t lim = 0.001;
02456 //     if((firstSol  >= 0.0) && (firstSol  <= 1.0)) { 
02457 //       for(Int_t i=0;i<3;i++) newPos[i] = start[i] + firstSol*(end[i]-start[i]);
02458 //       WCSimMathUtil::SubtractVec(pos,newPos,newvstart);
02459 //       WCSimMathUtil::SubtractVec(end,newPos,newdist);
02460 //       theta1 = WCSimMathUtil::DotProduct(newvstart,newdist)/(WCSimMathUtil::Magnitude(pos,newPos)*WCSimMathUtil::Magnitude(end,newPos));
02461 //
02462 //     }
02463 //     if((secondSol >= 0.0) && (secondSol <= 1.0)) {
02464 //       for(Int_t i=0;i<3;i++) newPos[i] = start[i] + secondSol*(end[i]-start[i]);
02465 //       WCSimMathUtil::SubtractVec(pos,newPos,newvstart);
02466 //       WCSimMathUtil::SubtractVec(end,newPos,newdist);
02467 //       theta2 = WCSimMathUtil::DotProduct(newvstart,newdist)/(WCSimMathUtil::Magnitude(pos,newPos)*WCSimMathUtil::Magnitude(end,newPos));
02468 //     }
02469 //     if((fabs(theta1-ctc) < lim) && (fabs(theta2-ctc) < lim)) return -1.0; // two good solutions
02470 //     if((fabs(theta1+ctc) < lim) && (fabs(theta2+ctc) < lim)) return -1.0; // two bad solutions
02471 //     if((fabs(theta1-ctc) < lim) && (theta2 <= 0)) {
02472 //       for(Int_t i=0;i<3;i++) intersection[i] = start[i] + firstSol*(end[i]-start[i]); 
02473 //       return WCSimMathUtil::Magnitude(pos,intersection);
02474 //     }  
02475 //     if(fabs(theta2-ctc) < lim && theta1 <= 0){
02476 //       for(Int_t i=0;i<3;i++) intersection[i] = start[i] + secondSol*(end[i]-start[i]); 
02477 //       return WCSimMathUtil::Magnitude(pos,intersection);
02478 //     }
02479 //     return -1.0; // no classification          
02480 //   } // end if
02481 //   else { 
02482 //     return -1.0;}
02483 //   return -1.0;    
02484 // }
02485 //
02486 // Double_t WCSimGeometry::DistanceToIntersectLine(Double_t *pos, Double_t *start, Double_t *end, Double_t *intersection)
02487 // {
02488 //   Double_t ctc    = WCSimParameters::CosThetaC;
02489 //   Double_t length = WCSimMathUtil::Magnitude(start,end);
02490 //   Double_t arg_a  = pow(length,4)*(1-WCSimMathUtil::sqr(ctc));
02491 //
02492 //   Double_t vstart[3], dist[3];    /// dist = vector(end - start), vstart = vector(pos - start) 
02493 //   WCSimMathUtil::SubtractVec(pos,start,vstart);
02494 //   WCSimMathUtil::SubtractVec(end,start,dist);
02495 //  
02496 //   Double_t arg_b  = -2*WCSimMathUtil::sqr(length)*WCSimMathUtil::DotProduct(vstart,dist)*(1-WCSimMathUtil::sqr(ctc));
02497 //  
02498 //   Double_t arg_c = WCSimMathUtil::sqr(WCSimMathUtil::DotProduct(vstart,dist))-WCSimMathUtil::sqr(WCSimMathUtil::Magnitude(pos,start))*WCSimMathUtil::sqr(length)*WCSimMathUtil::sqr(ctc);
02499 //  
02500 //   Double_t firstSol;
02501 //   Double_t secondSol;
02502 //   Double_t det = WCSimMathUtil::Determinant(arg_a, arg_b, arg_c);
02503 //  
02504 //   if (det >= 0) {
02505 //     WCSimMathUtil::RealRoots(arg_a, arg_b, sqrt(det), firstSol, secondSol);  // roots are real
02506 //  
02507 //     Double_t newPos[3], newvstart[3], newdist[3];
02508 //     Double_t costheta1 = 0.0;
02509 //     Double_t costheta2 = 0.0;
02510 //     Double_t lim = 0.001;
02511 //
02512 //     if((firstSol  >= 0.0) && (firstSol  <= 1.0)) { 
02513 //       for(Int_t i=0;i<3;i++) newPos[i] = start[i] + firstSol*(end[i]-start[i]);
02514 //       WCSimMathUtil::SubtractVec(pos,newPos,newvstart);
02515 //       WCSimMathUtil::SubtractVec(end,newPos,newdist);
02516 //       costheta1 = WCSimMathUtil::DotProduct(newvstart,newdist)/(WCSimMathUtil::Magnitude(pos,newPos)*WCSimMathUtil::Magnitude(end,newPos));
02517 //     }
02518 //     if((secondSol >= 0.0) && (secondSol <= 1.0)) {
02519 //       for(Int_t i=0;i<3;i++) newPos[i] = start[i] + secondSol*(end[i]-start[i]);
02520 //       WCSimMathUtil::SubtractVec(pos,newPos,newvstart);
02521 //       WCSimMathUtil::SubtractVec(end,newPos,newdist);
02522 //       costheta2 = WCSimMathUtil::DotProduct(newvstart,newdist)/(WCSimMathUtil::Magnitude(pos,newPos)*WCSimMathUtil::Magnitude(end,newPos));
02523 //     }
02524 //
02525 //     if((fabs(costheta1-ctc) < lim) && (fabs(costheta2-ctc) < lim)) { return -1.0; }// two good solutions
02526 //     if((fabs(costheta1+ctc) < lim) && (fabs(costheta2+ctc) < lim)) { return -1.0; }// two bad solutions
02527 //     if((fabs(costheta1-ctc) < lim) && (costheta2 == 0)) {
02528 //       for(Int_t i=0;i<3;i++) intersection[i] = start[i] + firstSol*(end[i]-start[i]);    
02529 //       return WCSimMathUtil::Magnitude(pos,intersection);
02530 //     }  
02531 //     if((fabs(costheta2-ctc) < lim) && (costheta1 == 0)){
02532 //       for(Int_t i=0;i<3;i++) intersection[i] = start[i] + secondSol*(end[i]-start[i]); 
02533 //       return WCSimMathUtil::Magnitude(pos,intersection);
02534 //     }
02535 //     return -1.0; // no classification          
02536 //   } // end if
02537 //   else {   
02538 //     return -1.0;
02539 //   }
02540 //   return -1.0;    
02541 // }
02542 //
02543 //
02544 
02545