WCSimIsoChronTransform.cc

Go to the documentation of this file.
00001 #include "WCSimIsoChronTransform.hh"
00002 #include "TObject.h"
00003 #include "TMath.h"
00004 #include "math.h"
00005 #include <iostream>
00006 #include <fstream>
00007 #include "TGraph.h"
00008 #include "WCSimRecoCluster.hh"
00009 #include "WCSimWaterModel.hh"
00010 #include "TVector3.h"
00011 #include <vector>
00012 
00013 using namespace std;
00014 
00015 ClassImp(WCSimIsoChronTransform);
00016 
00017 WCSimIsoChronTransform::WCSimIsoChronTransform(Double_t nBinsV, Double_t binscale){
00018 
00019   _n=1.34689;
00020   _c=299.792458;
00021 
00022   if((int)nBinsV%2==0){
00023 
00024     nbx=nBinsV;
00025     nby=nBinsV;
00026     nbz=nBinsV;
00027 
00028     double halfbin = binscale/2.;
00029     double nbinshalf = nBinsV/2.;
00030     
00031     lx = -(nbinshalf*binscale)-halfbin;
00032     hx = (nbinshalf*binscale)-halfbin;
00033     ly = -(nbinshalf*binscale)-halfbin;
00034     hy = (nbinshalf*binscale)-halfbin;
00035     lz = -(nbinshalf*binscale)-halfbin;
00036     hz = (nbinshalf*binscale)-halfbin;
00037 
00038     cout<<lx<<" "<<hx<<endl;
00039   }
00040   else{
00041 
00042     nbx=nBinsV;
00043     nby=nBinsV;
00044     nbz=nBinsV;
00045 
00046     double halfbin = binscale/2.;
00047     double nbinshalf = (nBinsV-1)/2.;
00048     
00049     lx = -(nbinshalf*binscale)-halfbin;
00050     hx = (nbinshalf*binscale)+halfbin;
00051     ly = -(nbinshalf*binscale)-halfbin;
00052     hy = (nbinshalf*binscale)+halfbin;
00053     lz = -(nbinshalf*binscale)-halfbin;
00054     hz = (nbinshalf*binscale)+halfbin;
00055 
00056     cout<<lx<<" "<<hx<<endl;
00057   }
00058   
00059   _htrackreco = new TH3D("soh","soh",nbx,lx,hx,nby,ly,hy,nbz,lz,hz);
00060   _thexb = new TH1D("xb","xb",nbx,lx,hx);
00061   _theyb = new TH1D("yb","yb",nby,ly,hy);
00062   _thezb = new TH1D("zb","zb",nbz,lz,hz);
00063 
00064   _maxAlphadist = new TH1D("maxalpha","maxalpha",200,-3.14159,3.14159);
00065   _maxThetaCdist = new TH1D("maxthetac","maxthetac",12000,-3.14159,3.14159);
00066   _bcontentsdist = new TH1D("bch","bch",50000,-0.5,50000.5);
00067   _minS1dist= new TH1D("ms1","ms1",50000,0.,5000.);
00068   _nDist = new TH1D("ndist","ndist",400,1.2,1.4);
00069   _S1vsAlphadist = new TH2D("s1va","s1va",5000,0.,5000.,5000,0,5000.);
00070 
00071   theprojxz = new TH2D("sohpxz","sohpxz",nbx,lx,hx,nbz,lz,hz);;
00072   theprojxy = new TH2D("sohpxy","sohpxy",nbx,lx,hx,nby,ly,hy);;
00073   theprojyz = new TH2D("sohpyz","sohpyz",nby,ly,hy,nbz,lz,hz);;
00074 
00075 }
00076 
00077 
00078 WCSimIsoChronTransform::WCSimIsoChronTransform(){
00079 
00080   _n=1.34689;
00081   _c=299.792458;
00082   
00083   nbx=110;
00084   nby=200;
00085   nbz=200;
00086   lx=-95.;
00087   hx=1005.;
00088   ly=-995.;
00089   hy=1005.;
00090   lz=-995.;
00091   hz=1005.;
00092   
00093   _htrackreco = new TH3D("soh","soh",nbx,lx,hx,nby,ly,hy,nbz,lz,hz);
00094   _thexb = new TH1D("xb","xb",nbx,lx,hx);
00095   _theyb = new TH1D("yb","yb",nby,ly,hy);
00096   _thezb = new TH1D("zb","zb",nbz,lz,hz);
00097 
00098   _maxAlphadist = new TH1D("maxalpha","maxalpha",200,-3.14159,3.14159);
00099   _maxThetaCdist = new TH1D("maxthetac","maxthetac",12000,-3.14159,3.14159);
00100   _bcontentsdist = new TH1D("bch","bch",50000,-0.5,50000.5);
00101   _minS1dist= new TH1D("ms1","ms1",50000,0.,5000.);
00102   _nDist = new TH1D("ndist","ndist",400,1.2,1.4);
00103   _S1vsAlphadist = new TH2D("s1va","s1va",5000,0.,5000.,5000,0,5000.);
00104 
00105   theprojxz = new TH2D("sohpxz","sohpxz",nbx,lx,hx,nbz,lz,hz);;
00106   theprojxy = new TH2D("sohpxy","sohpxy",nbx,lx,hx,nby,ly,hy);;
00107   theprojyz = new TH2D("sohpyz","sohpyz",nby,ly,hy,nbz,lz,hz);;
00108 
00109 }
00110 
00111   WCSimIsoChronTransform::~WCSimIsoChronTransform()
00112 {
00113 
00114   delete _htrackreco;
00115   delete _maxAlphadist;
00116   delete _maxThetaCdist;
00117   delete _bcontentsdist;
00118   delete _minS1dist;
00119   delete _S1vsAlphadist;
00120   delete _nDist;
00121   delete _thexb;
00122   delete _theyb;
00123   delete _thezb;
00124   delete theprojxz;
00125   delete theprojxy;
00126   delete theprojyz;
00127 
00128 }
00129 
00130 void Reset(){
00131 
00132 }
00133 
00134 
00135 Double_t WCSimIsoChronTransform::CalcMaxAlpha(Double_t dD, Double_t dT, Double_t ni){
00136 
00137    double qA=4*(ni*ni*ni*ni)*dD*dD;
00138    double qB=-8*(ni*ni)*(_c*dT)*dD;
00139    double qC= (4*dT*dT*_c*_c) - (4*((ni*ni) - 1)*( (ni*ni*dD*dD) - (dT*dT*_c*_c) ));
00140 
00141 
00142    double maxa=-555555;
00143    if( (qB*qB) > (4*qA*qC) ){
00144      double cosa = (-qB + sqrt(qB*qB - 4*qA*qC))/(2*qA);
00145      maxa = acos(cosa);
00146 
00147      double cosalphap = cos(maxa + 0.01);
00148      double qA1=((ni*ni)-1);
00149      double qB1=2*((_c*dT)-(ni*ni*dD*cosalphap));
00150      double qC1=(ni*ni*dD*dD)-(_c*_c*dT*dT);
00151 
00152      double cosalpham = cos(maxa - 0.01);
00153      qA1=((ni*ni)-1);
00154      qB1=2*((_c*dT)-(ni*ni*dD*cosalpham));
00155      qC1=(ni*ni*dD*dD)-(_c*_c*dT*dT);
00156 
00157    } else{
00158      //   cout<<"CalcMaxAlpha::Causally impossible "<<qB*qB<<" "<<(4*qA*qC)<<endl;
00159    }
00160 
00161    return maxa;
00162 }
00163 
00164 
00165 Double_t WCSimIsoChronTransform::CalcMaxEmissionAngle(Double_t dD, Double_t dT, Double_t ni){
00166 
00167   double maxalpha=this->CalcMaxAlpha(dD,dT,ni);
00168   double mins1 = this->CalcS1(maxalpha,dD,dT,ni);
00169   double ds2x = dD - mins1*cos(maxalpha);
00170   double ds2y = mins1*sin(maxalpha);
00171   double ds2 = sqrt( ds2x*ds2x + ds2y*ds2y );
00172   double sinwangle;
00173   if(ds2>0.) sinwangle = dD*sin(maxalpha)/ds2;
00174   double wangle = asin(sinwangle); 
00175    //  double maxemissangle = (3.14159265358979 - wangle);
00176   double maxemissangle=-55555;
00177   if(mins1>0) maxemissangle= wangle;
00178 
00179   //  cout<<sinwangle<<" "<<wangle<<" "<<maxemissangle<<" "<<mins1<<" "<<ds2<<endl;
00180  
00181   return maxemissangle;
00182 }
00183 
00184 Double_t WCSimIsoChronTransform::CalcAlpha(Double_t s1, Double_t dD, Double_t dT, Double_t ni){
00185 
00186   double s2 = (_c*dT -s1)/ni;
00187   double cosalpha;
00188   if(s1!=0) cosalpha = (s1*s1 + dD*dD - s2*s2)/(2*s1*dD);
00189   double alpha = acos(cosalpha);
00190 
00191   return alpha;
00192 }
00193 
00194 
00195 Double_t WCSimIsoChronTransform::CalcS1(Double_t alpha, Double_t dD, Double_t dT, Double_t ni)
00196 {
00197    double theS1=-55555.;
00198 
00199    if(alpha>=0){
00200      double cosalpha = cos(alpha);
00201      
00202      double qA=((ni*ni)-1);
00203      double qB=2*((_c*dT)-(ni*ni*dD*cosalpha));
00204      double qC=(ni*ni*dD*dD)-(_c*_c*dT*dT);
00205      double sqrtterm=qB*qB - 4*qA*qC;
00206      if( fabs(sqrtterm)<0.0001 ) sqrtterm=0;     
00207      if( sqrtterm>=0 ){    
00208 
00209        theS1 = (-qB + sqrt(sqrtterm))/(2*qA);
00210      }
00211      else{ 
00212        //       cout<<"CalcS1::Causally Impossible "<<(qB*qB)<<" "<<(4*qA*qC)<<" "<<(qB*qB - 4*qA*qC)<<endl; 
00213      }
00214    }
00215    return theS1;
00216 }
00217 
00218 Double_t WCSimIsoChronTransform::CalcS2(Double_t alpha, Double_t dD, Double_t dT, Double_t ni)
00219 {
00220   double mins1 = this->CalcS1(alpha,dD,dT,ni);
00221   double ds2x = dD - mins1*cos(alpha);
00222   double ds2y = mins1*sin(alpha);
00223   double theS2 = sqrt( ds2x*ds2x + ds2y*ds2y );
00224 
00225   return theS2;
00226 }
00227 
00228 
00229 //Double_t WCSimIsoChronTransform::FindIndex(Double_t dD, Double_t dT)
00230 TGraph* WCSimIsoChronTransform::IndexVsS1(Double_t dD, Double_t dT)
00231 {
00232 
00233   double minimalS1=1000000.;
00234   double minimalIndex=0.;
00235   double theindex_i=1.30;
00236 
00237   TGraph* s1vsindex = new TGraph(100);
00238 
00239   for(int i=0; i<100; i++)
00240     {
00241       double maxalpha_i=this->CalcMaxAlpha(dD,dT, theindex_i);
00242       double mins1_i=this->CalcS1(maxalpha_i,dD,dT, theindex_i);
00243       // double eang_i=this->CalcMaxEmissionAngle(dD,dT);
00244       bool past_thresh=false;
00245       s1vsindex->SetPoint(i,theindex_i,mins1_i);
00246       //      cout<<"fitting..."<<theindex_i<<" "<<mins1_i<<" "<<maxalpha_i<<" "<<eang_i<<" "<<acos(1/theindex_i)<<endl;
00247 
00248       theindex_i+=0.001;
00249     }
00250  
00251   return s1vsindex;
00252 }
00253 
00254 TH3D* WCSimIsoChronTransform::IsoChronTransform(WCSimRecoCluster* theHits)
00255 {
00256 
00257   return _htrackreco;
00258 }
00259 
00260 
00261 void WCSimIsoChronTransform::ApplyTransform(vector< vector<double> >  thehits, vector<double> hypVtx)
00262 {
00263   for(int i=0; i<(int)(thehits.size()); i++){
00264 
00265     if(i%1000==0) cout<<"reconstructing from hits...hitcount: "<<i<<endl;
00266 
00267     vector<double> hcoo = thehits.at(i);
00268     this->ProcessHit(hcoo,hypVtx);
00269   }
00270 }
00271 
00272 void WCSimIsoChronTransform::ApplyTransformChromCorrected(vector< vector<double> >  thehits, vector<double> hypVtx, Int_t ncolorbins, Double_t lowcolor, Double_t highcolor)
00273 {
00274   if( (!iswatermodel) || (ncolorbins==0) ){
00275     if(!iswatermodel){
00276       cout<<"need to build WCSimWaterModel before calling this routine!!!"<<endl;;
00277     }
00278     if(ncolorbins==0){
00279       cout<<"need more than 0 color bins to run this routine!!!"<<endl;
00280     }
00281   }
00282   else{
00283 
00284     for(int i=0; i<(int)(thehits.size()); i++){
00285       
00286       if(i%1000==0) cout<<"reconstructing from hits...hitcount: "<<i<<endl;
00287       vector<double> hcoo = thehits.at(i);
00288       this->ProcessHitChromCorr(hcoo,hypVtx,ncolorbins,lowcolor,highcolor);      
00289     }
00290   }
00291 }
00292 
00293 
00294   void WCSimIsoChronTransform::ProcessHit(vector<double> hitcoordinates, vector<double> hypVtx)
00295 {
00296 
00297   TVector3 hitcoor(hitcoordinates.at(0),hitcoordinates.at(1),hitcoordinates.at(2));
00298   TVector3 hvect((hitcoordinates.at(0)-hypVtx.at(0)),(hitcoordinates.at(1)-hypVtx.at(1)),(hitcoordinates.at(2)-hypVtx.at(2)));
00299 
00300   double dD=hvect.Mag();
00301   double dT=hitcoordinates.at(3)-hypVtx.at(3);
00302   double theR; 
00303   double minalpha=0;
00304   double maxalpha = this->CalcMaxAlpha(dD,dT,_n);
00305   double maxemmang = this->CalcMaxEmissionAngle(dD,dT,_n);  
00306   double minS1val = this->CalcS1(maxalpha,dD,dT,_n); 
00307 
00308   _maxAlphadist->Fill(maxalpha);
00309   _maxThetaCdist->Fill(maxemmang);
00310   _minS1dist->Fill(minS1val);
00311   _S1vsAlphadist->Fill((minS1val*cos(maxalpha)),dD);
00312 
00313   vector<double> priorpixel;
00314   priorpixel.push_back(0.); priorpixel.push_back(0.); priorpixel.push_back(0.);
00315   priorpixel.push_back(0.); priorpixel.push_back(0.); priorpixel.push_back(0.);
00316   
00317   double alphastep=(maxalpha-minalpha)/1000.;
00318   theR=0;
00319 
00320   if(maxalpha>0){
00321     for(int i=999; i<1000; i++){
00322  
00323       double thealpha=minalpha+(i*alphastep);
00324       double theR = this->CalcS1(thealpha,dD,dT,_n);
00325       
00326       vector<double> newpp = this->FillIsoChronRing_pixel(hypVtx, hvect, priorpixel, thealpha, theR, 1.);
00327       priorpixel = newpp;
00328     }
00329   }
00330 }
00331 
00332 
00333 void WCSimIsoChronTransform::ProcessHitChromCorr(vector<double> hitcoordinates, vector<double> hypVtx, Int_t ncolorbins, Double_t lowcolor, Double_t highcolor)
00334 {
00335 
00336   TVector3 hitcoor(hitcoordinates.at(0),hitcoordinates.at(1),hitcoordinates.at(2));
00337   TVector3 hvect((hitcoordinates.at(0)-hypVtx.at(0)),(hitcoordinates.at(1)-hypVtx.at(1)),(hitcoordinates.at(2)-hypVtx.at(2)));
00338 
00339   double dD=hvect.Mag();
00340   double dT=hitcoordinates.at(3)-hypVtx.at(3);
00341   double theR; 
00342   double minalpha=0;
00343 
00344   double maxalpha = this->CalcMaxAlpha(dD,dT,_n);
00345   double maxemmang = this->CalcMaxEmissionAngle(dD,dT,_n);  
00346   double minS1val = this->CalcS1(maxalpha,dD,dT,_n); 
00347 
00348   _maxAlphadist->Fill(maxalpha);
00349   _maxThetaCdist->Fill(maxemmang);
00350   _minS1dist->Fill(minS1val);
00351   _S1vsAlphadist->Fill((minS1val*cos(maxalpha)),dD);
00352 
00353   vector<double> priorpixel;
00354   priorpixel.push_back(0.); priorpixel.push_back(0.); priorpixel.push_back(0.);
00355   priorpixel.push_back(0.); priorpixel.push_back(0.); priorpixel.push_back(0.);
00356   
00357   double alphastep=(maxalpha-minalpha)/1000.;
00358   theR=0;
00359 
00360   Double_t thecolor_i = lowcolor;
00361   Double_t colorincrement = (highcolor-lowcolor)/((Double_t)ncolorbins);
00362 
00363   vector<double> s2vect;
00364   vector<double> s1vect;
00365   vector<double> alphavect;
00366   double colorNorm=0;
00367 
00368   for(int i=0; i<ncolorbins; i++){
00369     double theindex_i=_theWM->N_Index(thecolor_i);
00370     double malpha = this->CalcMaxAlpha(dD,dT,theindex_i);
00371     double s_2 = this->CalcS2(malpha,dD,dT,theindex_i);
00372     s2vect.push_back(s_2);
00373     s1vect.push_back(this->CalcS1(malpha,dD,dT,theindex_i));
00374     alphavect.push_back(malpha);
00375     colorNorm+=(_theWM->FinlSpect(thecolor_i,s_2));
00376     thecolor_i+=colorincrement;
00377   }
00378   
00379   if(colorNorm>0){
00380     thecolor_i = lowcolor;
00381     for(int i=0; i<ncolorbins; i++){
00382       double theindex_i=_theWM->N_Index(thecolor_i);
00383       double theweight_i=(_theWM->FinlSpect(thecolor_i,(s2vect.at(i))))/colorNorm;
00384       _nDist->Fill(theindex_i,theweight_i);
00385       
00386       vector<double> newpp = this->FillIsoChronRing_pixel(hypVtx, hvect, priorpixel, (alphavect.at(i)), (s1vect.at(i)), theweight_i);
00387       priorpixel = newpp;
00388       
00389       thecolor_i+=colorincrement;
00390     }
00391   }
00392   else{
00393     cout<<"Color Norm is Zero!!!"<<endl;
00394   }  
00395 }
00396 
00397 
00398 vector<double> WCSimIsoChronTransform::FillIsoChronRing_pixel(vector<double> VxHypothesis, TVector3 hvect, vector<double> priorpixel, Double_t thealpha, Double_t theR, Double_t theweight)
00399 {
00400   double thebeta=0;
00401   double xp=theR*cos(thealpha);
00402   double yp=theR*sin(thealpha);
00403   
00404   double pmag = hvect.Mag();
00405   TVector3 paradir = (1/pmag)*hvect;
00406   TVector3 yvect(0.,1.,0.);
00407   TVector3 pdu = (paradir.Cross(yvect));
00408   double perpmag = pdu.Mag();
00409   TVector3 perpdir = (1/perpmag)*pdu;
00410   TVector3 threevtx(VxHypothesis.at(0),VxHypothesis.at(1),VxHypothesis.at(2));
00411   
00412   vector<double> newpriorpixel;
00413   newpriorpixel.push_back(0.); newpriorpixel.push_back(0.); newpriorpixel.push_back(0.);
00414   newpriorpixel.push_back(0.); newpriorpixel.push_back(0.); newpriorpixel.push_back(0.);
00415 
00416 
00417   for(int j=0; j<200; j++){
00418         
00419     thebeta = (j*(3.14159265358979))/100.;
00420     
00421     TVector3 f1pos = (xp*paradir - yp*cos(thebeta)*perpdir + yp*sin(thebeta)*yvect);
00422     TVector3 f1neg = (xp*paradir + yp*cos(thebeta)*perpdir + yp*sin(thebeta)*yvect);
00423     
00424     TVector3 f2pos = (hvect-f1pos);
00425     TVector3 f2neg = (hvect-f1neg);
00426 
00427     f1pos+=threevtx;
00428     f1neg+=threevtx;
00429     f2pos+=threevtx;
00430     f2neg+=threevtx;
00431         
00432     int ixbp = _thexb->FindBin(f1pos.X());
00433     int iybp = _theyb->FindBin(f1pos.Y());
00434     int izbp = _thezb->FindBin(f1pos.Z());
00435     int ixbn = _thexb->FindBin(f1neg.X());
00436     int iybn = _theyb->FindBin(f1neg.Y());
00437     int izbn = _thezb->FindBin(f1neg.Z());
00438     
00439     if( (ixbp!=priorpixel.at(0)) || (iybp!=priorpixel.at(1)) || (izbp!=priorpixel.at(2)) || (ixbn!=priorpixel.at(3)) || (iybn!=priorpixel.at(4)) || (izbn!=priorpixel.at(5))){
00440       
00441       double bc =_htrackreco->GetBinContent(ixbp,iybp,izbp);
00442       _htrackreco->SetBinContent(ixbp,iybp,izbp,(bc+theweight));
00443       
00444       //      bc = _htrackreco->GetBinContent(ixbn,iybn,izbn);
00445       //      _htrackreco->SetBinContent(ixbn,iybn,izbn,(bc+theweight));
00446       newpriorpixel.at(0)=ixbp;
00447       newpriorpixel.at(1)=iybp;
00448       newpriorpixel.at(2)=izbp;
00449       
00450       newpriorpixel.at(3)=ixbn;
00451       newpriorpixel.at(4)=iybn;
00452       newpriorpixel.at(5)=izbn;
00453       
00454       priorpixel = newpriorpixel;
00455     }
00456   }
00457 
00458   return newpriorpixel;
00459 }
00460 
00461 
00462 TH2D* WCSimIsoChronTransform::XYProjection()
00463 {
00464   return theprojxy;
00465 }
00466 
00467 TH2D* WCSimIsoChronTransform::XZProjection(double bclimit)
00468 {
00469   cout<<"starting XZProjection..."<<endl;
00470 
00471    for(int i=0; i<nbx; i++){
00472      //     cout<<i<<endl;
00473      for(int j=0; j<nbz; j++){
00474 
00475        double zsum=0;
00476        double zmax=0;
00477 
00478        for(int k=0; k<nby; k++){
00479 
00480          double nbc = _htrackreco->GetBinContent((i+1),(k+1),(j+1));
00481          //      cout<<"here now..."<<i<<" "<<j<<" "<<k<<" "<<nbc<<endl;
00482          if(nbc>=0) _bcontentsdist->Fill(nbc);
00483          //      if(nbc>0) cout<<i<<" "<<j<<" "<<k<<" "<<nbc<<endl;
00484          if(nbc>bclimit) zsum+=nbc;
00485          if(nbc>zmax) zmax=nbc;
00486        }
00487 
00488        //       if(zsum==0) zsum=1;
00489        theprojxz->SetBinContent((i+1),(j+1),zsum);
00490        //      if(zmax>400) sohp->SetBinContent((i+1),(j+1),zmax);
00491        //       else sohp->SetBinContent((i+1),(j+1),2);
00492        //       if(zsum==0) sohp->SetBinContent((i+1),(j+1),2);
00493      }
00494    }
00495 
00496   return theprojxz;
00497 }
00498 
00499 TH2D* WCSimIsoChronTransform::YZProjection()
00500 {
00501   return theprojyz;
00502 }
00503 
00504 TH2D* WCSimIsoChronTransform::XYSlice()
00505 {
00506   TH2D* theslice;
00507   return theslice;
00508 }
00509 
00510 TH2D* WCSimIsoChronTransform::XZSlice()
00511 {
00512   TH2D* theslice;
00513   return theslice;
00514 }
00515 
00516 TH2D* WCSimIsoChronTransform::YZSlice()
00517 {
00518   TH2D* theslice;
00519   return theslice;
00520 }
00521 
00522 TH1D* WCSimIsoChronTransform::getBCDist()
00523 {
00524   return _bcontentsdist;
00525 }
00526 
00527 TH1D* WCSimIsoChronTransform::getMaxAlphaDist()
00528 {
00529   return _maxAlphadist;
00530 }
00531 
00532 TH1D* WCSimIsoChronTransform::getMaxThetaCDist()
00533 {
00534   return _maxThetaCdist;
00535 }
00536 
00537 
00538 TH1D* WCSimIsoChronTransform::getMinS1Dist()
00539 {
00540   return _minS1dist;
00541 }
00542 
00543 TH1D* WCSimIsoChronTransform::getnDist()
00544 {
00545   return _nDist;
00546 }
00547 
00548 TH2D* WCSimIsoChronTransform::getS1vsAlpha()
00549 {
00550   return _S1vsAlphadist;
00551 }
00552 
00553 void WCSimIsoChronTransform::SetConstantIndexofRefraction(double newn)
00554 {
00555   _n=newn;
00556 }
00557 
00558 void WCSimIsoChronTransform::SetConstantSpeedofParticle(double newc)
00559 {
00560   _c=newc;
00561 }
00562 
00563 void WCSimIsoChronTransform::AddWaterModel(WCSimWaterModel* wm)
00564 {
00565   iswatermodel=true;
00566   _theWM = wm;
00567 }