loon: AltDeMuxCalc.cxx Source File

00001 
00002 //$Id: AltDeMuxCalc.cxx,v 1.6 2006/12/01 20:39:57 rhatcher Exp $
00003 //
00004 //AltDeMuxCalc.cxx
00005 //
00006 // Some common calculations for AltDeMux
00007 //
00008 // Author:  M. Thosmson April 2003
00010 
00011 #include "AltDeMuxCalc.h"
00012 #include "MessageService/MsgService.h"
00013 #include "Conventions/CalTimeType.h"
00014 #include "Conventions/StripEnd.h"
00015 #include <cmath>
00016 
00017 ClassImp(AltDeMuxCalc)
00018 
00019 CVSID("$Id: AltDeMuxCalc.cxx,v 1.6 2006/12/01 20:39:57 rhatcher Exp $");
00020 
00021 //......................................................................
00022 AltDeMuxCalc::AltDeMuxCalc() 
00023 {
00024 
00025   for(Int_t ip=0;ip<MAX_NUMBER_OF_PLANES;ip++){
00026     for(Int_t is=0;is<MAX_NUMBER_OF_STRIPS;is++){
00027       fWLSFibreLengthE[ip][is] = -999.;
00028       fClearFibreLengthE[ip][is] = -999.;
00029       fWLSFibreLengthW[ip][is] = -999.;
00030       fClearFibreLengthW[ip][is] = -999.;
00031     }
00032   }
00033 
00034 
00035   // build the look up table for attenuation corrections
00036   // store in 1cm steps - to a maximum of 25m  
00037   for(Int_t iwls=0; iwls<2500;iwls++){
00038     Float_t lWls = static_cast<Float_t>(iwls);
00039     lWls = lWls/100.;
00040     fAttCorWLS[iwls] = 1.0/( 0.666*exp(-lWls/7.05)+ 0.333*exp(-(lWls)/1.05));
00041   }
00042   for(Int_t iclear=0; iclear<2500;iclear++){
00043     Float_t lClear = static_cast<Float_t>(iclear);
00044     lClear = lClear/100.;
00045     fAttCorClear[iclear] = exp(+(lClear)/10.);
00046   }
00047   return;
00048 }
00049 
00050 //----------------------------------------------------------------------
00051 AltDeMuxCalc::~AltDeMuxCalc()
00052 {
00053 
00054 }
00055 
00056 //......................................................................
00057 
00058 Float_t AltDeMuxCalc::AttCorWLS(Float_t lwls)
00059 { 
00060   
00061   Int_t iwls = static_cast<Int_t>(lwls*100.);
00062   if(iwls>=2500){
00063     MSG("AlgAltDeMux", Msg::kDebug) << "AltDeMuxCalc::AttWLS requested correction for too long WLS fibre" << endl;
00064     iwls = 2500;
00065   }
00066 
00067   return fAttCorWLS[iwls];
00068 }
00069 
00070 
00071 Float_t AltDeMuxCalc::AttCorClear(Float_t lclear)
00072 { 
00073   
00074   Int_t iclear = static_cast<Int_t>(lclear*100.);
00075   if(iclear>=2500){
00076     MSG("AlgAltDeMux", Msg::kDebug) << "AltDeMuxCalc::AttClear requested correction for too long Clear fibre" << endl;
00077     iclear = 2500;
00078   }
00079 
00080   return fAttCorClear[iclear];
00081 }
00082 
00083 
00084 void AltDeMuxCalc::SetEast(PlexSEIdAltL* pE, Int_t stripE)
00085 { 
00086   
00087   pAltListE = pE;
00088   fStripE = stripE;
00089   fEastIsSet = false;
00090 
00091   return;
00092 }
00093 
00094 void AltDeMuxCalc::SetWest(PlexSEIdAltL* pW, Int_t stripW)
00095 { 
00096   
00097   pAltListW = pW;
00098   fStripW = stripW;
00099   fWestIsSet = false;
00100 
00101   return;
00102 }
00103 
00104 
00105 void AltDeMuxCalc::SetWestToStrip(PlexSEIdAltL* pW, Int_t stripW)
00106 { 
00107 
00108   bool notFound = true;
00109   pW->SetFirst();
00110   while( pW->IsValid() && notFound ){
00111     Int_t is = pW->GetCurrentSEId().GetStrip();
00112     if(is==stripW)notFound=false;
00113     if(notFound)pW->Next();
00114   }
00115 
00116   if(notFound)MSG("AltDeMux", Msg::kError) << "AltDeMuxCalc::SetWestToStrip Could not find strip " << stripW << " in PlexSEIdAltL " << endl;
00117   
00118   pAltListW = pW;
00119   fStripW = stripW;
00120   fWestIsSet = false;
00121 
00122   return;
00123 }
00124 
00125 void AltDeMuxCalc::SetEastToStrip(PlexSEIdAltL* pE, Int_t stripE)
00126 { 
00127   bool notFound = true;
00128   pE->SetFirst();
00129   while( pE->IsValid() && notFound ){
00130     Int_t is = pE->GetCurrentSEId().GetStrip();
00131     if(is==stripE)notFound=false;
00132     if(notFound)pE->Next();
00133   }
00134 
00135   if(notFound)MSG("AltDeMux", Msg::kError) << "AltDeMuxCalc::SetEastToStrip Could not find strip " << stripE << " in PlexSEIdAltL " << endl;
00136   
00137   pAltListE = pE;
00138   fStripE = stripE;
00139   fEastIsSet = false;
00140 
00141   return;
00142 }
00143 
00144 
00145 void AltDeMuxCalc::CalcEastWest()
00146 {
00147   this->CalcEast();
00148   this->CalcWest();
00149 
00150   return;
00151 }
00152 
00153 void AltDeMuxCalc::CalcBestEastWest()
00154 {
00155   this->CalcBestEast();
00156   this->CalcBestWest();
00157   return;
00158 }
00159 
00160 
00161 Float_t AltDeMuxCalc::CurrentQ(PlexSEIdAltL* pAltList)
00162 {
00163 
00164   Float_t Q  = pAltList->GetCurrentItem().GetPE();
00165   Float_t Qc = pAltList->GetCurrentItem().GetSigCorr()/60.;
00166   if(Q/Qc<0.1 || Q/Qc > 10.0){
00167     Qc=Q;
00168   }
00169   
00170   return Qc;
00171 }
00172 
00173 void AltDeMuxCalc::CalcEast()
00174 {
00175 
00176   if(fEastIsSet)return;
00177 
00178   fQE  = pAltListE->GetCurrentItem().GetPE();
00179   fQEc = pAltListE->GetCurrentItem().GetSigCorr()/60.;
00180   if(fQE/fQEc<0.1 || fQE/fQEc > 10.0){
00181     MSG("AltDeMux", Msg::kDebug) << "AltDeMuxCalc::CalcEast   Large/small SigCor/SigPE ratio : " <<fQEc/fQE << " for Plane : " << fPlane << endl;
00182     fQEc=fQE;
00183   }
00184   
00185   fClearE = fClearFibreLengthE[fPlane][fStripE]; 
00186   fWlsE   = fWLSFibreLengthE[fPlane][fStripE];  
00187 
00188   // timewalk correction based on adc - ok so I cheat and use SigLin
00189   // which is OK since corrections neg for high adc
00190   float adc  = pAltListE->GetCurrentItem().GetSigLin();
00191   fTE     = pAltListE->GetCurrentItem().GetTime()*1.0E9 - this->TimeWalk(adc);
00192 
00193   fEastIsSet=true;
00194   fDtIsSet = false;
00195   fDQIsSet = false;
00196 
00197   return;
00198 }
00199 
00200 
00201 void AltDeMuxCalc::CalcWest()
00202 {
00203 
00204   if(fWestIsSet)return;
00205 
00206   fQW  = pAltListW->GetCurrentItem().GetPE();
00207   fQWc = pAltListW->GetCurrentItem().GetSigCorr()/60.;
00208   if(fQW/fQWc<0.1 || fQW/fQWc > 10.0){
00209 
00210     MSG("AltDeMux", Msg::kDebug) << "AltDeMuxCalc::CalcWest Large/small SigCor/SigPE ratio : " <<fQWc/fQW << " for Plane : " << fPlane << endl;
00211     fQWc=fQW;
00212   }
00213    
00214   fClearW = fClearFibreLengthW[fPlane][fStripW]; 
00215   fWlsW   = fWLSFibreLengthW[fPlane][fStripW];  
00216 
00217   // timewalk correction based on adc - ok so I cheat and use SigLin
00218   // which is OK since corrections neg for high adc
00219   float adc  = pAltListW->GetCurrentItem().GetSigLin();
00220   fTW        = pAltListW->GetCurrentItem().GetTime()*1.0E9 - this->TimeWalk(adc);
00221 
00222   fWestIsSet=true;
00223   fDtIsSet = false;
00224   fDQIsSet = false;
00225 
00226   return;
00227 }
00228 
00229 
00230 void AltDeMuxCalc::CalcBestEast()
00231 {
00232 
00233   if(fEastIsSet)return;
00234 
00235   fQE  = pAltListE->GetBestItem().GetPE();
00236   fQEc = pAltListE->GetBestItem().GetSigCorr()/60.;
00237   if(fQE/fQEc<0.1 || fQE/fQEc > 10.0){
00238     MSG("AltDeMux", Msg::kDebug) << "AltDeMuxCalc::CalcBestEast   Large/small SigCor/SigPE ratio : " <<fQEc/fQE << " for Plane : " << fPlane << endl;
00239     fQEc=fQE;
00240   }
00241 
00242   
00243   fClearE = fClearFibreLengthE[fPlane][fStripE]; 
00244   fWlsE   = fWLSFibreLengthE[fPlane][fStripE];  
00245 
00246   // timewalk correction based on adc - ok so I cheat and use SigLin
00247   // which is OK since corrections neg for high adc
00248   float adc   = pAltListE->GetBestItem().GetSigLin();
00249   fTE     = pAltListE->GetBestItem().GetTime()*1.0E9-this->TimeWalk(adc);
00250   fEastIsSet=true;
00251   fDtIsSet = false;
00252   fDQIsSet = false;
00253 
00254   return;
00255 }
00256 
00257 void AltDeMuxCalc::CalcBestWest()
00258 {
00259 
00260   if(fWestIsSet)return;
00261 
00262   fQW  = pAltListW->GetBestItem().GetPE();
00263   fQWc = pAltListW->GetBestItem().GetSigCorr()/60.;
00264   if(fQW/fQWc<0.1 || fQW/fQWc > 10.0){
00265     MSG("AltDeMux", Msg::kDebug) << "AltDeMuxCalc::CalcBestWest   Large/small SigCor/SigPE ratio : " <<fQWc/fQW << " for Plane : " << fPlane << endl;
00266     fQWc=fQW;
00267   }
00268 
00269   fClearW = fClearFibreLengthW[fPlane][fStripW]; 
00270   fWlsW   = fWLSFibreLengthW[fPlane][fStripW];  
00271 
00272   float adc = pAltListW->GetBestItem().GetSigLin();
00273   fTW       = pAltListW->GetBestItem().GetTime()*1.0E9- this->TimeWalk(adc);
00274 
00275   fWestIsSet=true;
00276   fDtIsSet = false;
00277   fDQIsSet = false;
00278 
00279   return;
00280 }
00281 
00282 
00283 
00284 void AltDeMuxCalc::CalcDt()
00285 {
00286 
00287   if(!fEastIsSet)this->CalcEast();
00288   if(!fWestIsSet)this->CalcWest();
00289   fDt = (fTW -fTE)-(fClearW-fClearE)/fClearFibreC-(fWlsW-fWlsE)/fWLSFibreC;
00290   fDtIsSet = true;
00291 
00292   return;
00293  
00294 }
00295 
00296 Int_t AltDeMuxCalc::StripAim()
00297 { 
00298 
00299   double stripAim=0.0;
00300   if(!fDtIsSet)this->CalcDt();
00301   
00302   if(fView==PlaneView::kU){
00303     stripAim = fNumberOfStrips/2.0+fDt*fDt2Dstrip;
00304   }
00305   if(fView==PlaneView::kV){
00306     stripAim = fNumberOfStrips/2.0-fDt*fDt2Dstrip;
00307   }
00308 
00309   if(stripAim<0.)fStripAim=0;
00310   if(stripAim>=fNumberOfStrips)fStripAim=fNumberOfStrips-1;
00311   if(stripAim>=0.&&stripAim<fNumberOfStrips)fStripAim = static_cast<Int_t>(stripAim);     
00312 
00313   return fStripAim;
00314 
00315 }
00316 
00317 Float_t AltDeMuxCalc::SigmaDQ()
00318 { 
00319 
00320   if(!fDtIsSet)this->CalcDt();
00321   if(fDQIsSet)return fDQ;
00322   Float_t le = 4.0-0.0823*fDt;
00323   Float_t lw = 4.0+0.0823*fDt;
00324   if(lw>8.0)lw = 8.0;
00325   if(le>8.0)le = 8.0;
00326   if(lw<0.0)lw = 0.0;
00327   if(le<0.0)le = 0.0;
00328   //Float_t ae = exp(+(fClearE)/10.)/( 0.666*exp(-(le+fWlsE)/7.05)+ 0.333*exp(-(le+fWlsE)/1.05));
00329   //Float_t aw = exp(+(fClearW)/10.)/( 0.666*exp(-(lw+fWlsW)/7.05)+ 0.333*exp(-(lw+fWlsW)/1.05));
00330 
00331   Float_t ae = this->AttCorWLS(le+fWlsE)*this->AttCorClear(fClearE);
00332   Float_t aw = this->AttCorWLS(lw+fWlsW)*this->AttCorClear(fClearW);
00333   Float_t ae2 = ae*ae;
00334   Float_t aw2 = aw*aw;
00335   fQAttCorE = ae*fQEc;
00336   fQAttCorW = aw*fQWc;
00337 
00338   Float_t sl1 = 0.5/6.6*(fQEc*this->AttCorWLS(le)+fQWc*this->AttCorWLS(lw));
00339   Float_t d1  = fQAttCorE - fQAttCorW;
00340   Float_t sd1 = sqrt(fQEc*ae2 + fQWc*aw2 + sl1*sl1 + fQEc*fQEc*ae2/50. + fQWc*fQWc*aw2/50.);
00341   sd1 = d1/sd1;
00342   fDQ = sd1;
00343   fDQIsSet = true;
00344   return sd1;
00345 
00346 }
00347 
00348 
00349 
00350 
00351 
00352 
00353 void AltDeMuxCalc::SetFibreLengthE(Int_t iplane, PlexSEIdAltL* paltlist)
00354 {
00355 
00356   Int_t istrip = paltlist->GetCurrentSEId().GetStrip();
00357   if(fClearFibreLengthE[iplane][istrip]>0)return;
00358 
00359   UgliStripHandle  ush = pUgh->GetStripHandle(paltlist->GetCurrentSEId());
00360   if(ush.IsValid()){
00361     fClearFibreLengthE[iplane][istrip] = ush.ClearFiber(StripEnd::kEast);
00362     fWLSFibreLengthE[iplane][istrip]   = ush.WlsPigtail(StripEnd::kEast);
00363   }else{
00364     MSG("AltDeMux", Msg::kFatal) << "AltDeMuxCalc::SetFibreLengthE => UgliStripHandle NOT VALID : plane " << iplane << ":" << istrip << " (Carry on regardless)" <<  endl;
00365   }
00366   
00367   return;
00368 
00369 }
00370 
00371 void AltDeMuxCalc::SetFibreLengthW(Int_t iplane, PlexSEIdAltL* paltlist)
00372 {
00373 
00374   Int_t istrip = paltlist->GetCurrentSEId().GetStrip();
00375   if(fWLSFibreLengthW[iplane][istrip]>0)return;
00376 
00377   UgliStripHandle ush = pUgh->GetStripHandle(paltlist->GetCurrentSEId());
00378   if(ush.IsValid()){
00379     fClearFibreLengthW[iplane][istrip] = ush.ClearFiber(StripEnd::kWest);
00380     fWLSFibreLengthW[iplane][istrip]   = ush.WlsPigtail(StripEnd::kWest);
00381   }else{
00382     MSG("AltDeMux", Msg::kFatal) << "AltDeMuxCalc::SetFibreLengthW => UgliStripHandle NOT VALID : plane " << iplane << ":" << istrip << " (Carry on regardless)" <<  endl;
00383   }
00384   return;
00385 
00386 }
00387 
00388 
00389 
00390 void AltDeMuxCalc::FindNearestEdge(Float_t x, Float_t y, Int_t &iedge, Float_t &d)
00391 {
00392   
00393   Float_t d1 = fabs(x);
00394   Float_t d2 = fabs(y);
00395   Float_t d3 = 0.70711*fabs(x-y);
00396   Float_t d4 = 0.70711*fabs(x+y);
00397 
00398   if(d1>d2){
00399     if(d3>d4){
00400       if(d1>d3){
00401         d = 4.0-d1;
00402         if(x>0)iedge=3;
00403         if(x<=0)iedge=7;
00404       }else{
00405         d = 4.0-d3;
00406         if(x-y>0)iedge=4;
00407         if(x-y<=0)iedge=8;
00408       }
00409     }else{
00410       if(d1>d4){
00411         d = 4.0-d1;
00412         if(x>0)iedge=3;
00413         if(x<=0)iedge=7;
00414       }else{
00415         d = 4.0-d4;
00416         if(x+y>0)iedge=2;
00417         if(x+y<=0)iedge=6;
00418       }
00419     }
00420   }else{
00421     if(d3>d4){
00422       if(d2>d3){
00423         d = 4.0-d2;
00424         if(y>0)iedge=1;
00425         if(y<=0)iedge=5;
00426 
00427       }else{
00428         d = 4.0-d3;
00429         if(x-y>0)iedge=4;
00430         if(x-y<=0)iedge=8;
00431       }
00432     }else{
00433       if(d2>d4){
00434         d = 4.0-d2;
00435         if(y>0)iedge=1;
00436         if(y<=0)iedge=5;
00437       }else{
00438         d = 4.0-d4;
00439         if(x+y>0)iedge=2;
00440         if(x+y<=0)iedge=6;
00441       }
00442     }
00443   }
00444 
00445   return;
00446 } 
00447 
00448 
00449 
00450 Float_t AltDeMuxCalc::TimeWalk(float adc){
00451 
00452   // Time walk correction based on DATA parameterisation
00453   // from Andy Blake (new fit using form suggested by Roy Lee)
00454 
00455   Float_t logQ  = log(adc/2.3);
00456   Float_t logQ2 = logQ*logQ;
00457   Float_t logQ3 = logQ2*logQ;
00458   Float_t timeW = 5.7 - 0.76*logQ-0.038*logQ2+0.00541*logQ3; 
00459   return  timeW/0.3;
00460 
00461 }