loon: CandShowerHandle.cxx Source File

00001 
00002 // $Id: CandShowerHandle.cxx,v 1.52 2008/11/05 17:33:15 rodriges Exp $
00003 //
00004 // CandShowerHandle.cxx
00005 //
00006 // CandShowerHandle is the specialized access handle to CandShower.
00007 //
00008 // Each concrete CandHandle must define a DupHandle function.
00009 //
00010 // Author:  R. Lee
00012 
00013 #include "TMath.h"
00014 
00015 #include <cassert>
00016 #include <cmath>
00017 #include <iostream>
00018 
00019 #include "MessageService/MsgService.h"
00020 #include "RecoBase/CandClusterHandle.h"
00021 #include "RecoBase/CandSliceHandle.h"
00022 #include "RecoBase/CandShowerHandle.h"
00023 #include "RecoBase/CandTrackHandle.h"
00024 #include "RecoBase/CandFitTrackHandle.h"
00025 #include "RecoBase/CandShower.h"
00026 #include "Navigation/NavKey.h"
00027 #include "Navigation/NavSet.h"
00028 #include "Validity/VldContext.h"
00029 #include "Algorithm/AlgConfig.h"
00030 #include "UgliGeometry/UgliGeomHandle.h"
00031 #include "UgliGeometry/UgliScintPlnHandle.h"
00032 #include "DataUtil/PlaneOutline.h"
00033 #include "DataUtil/EnergyCorrections.h"
00034 
00035 ClassImp(CandShowerHandle)
00036 
00037 //______________________________________________________________________
00038 CVSID("$Id: CandShowerHandle.cxx,v 1.52 2008/11/05 17:33:15 rodriges Exp $");
00039 
00040 //______________________________________________________________________
00041 CandShowerHandle::CandShowerHandle()
00042 {
00043 }
00044 
00045 //______________________________________________________________________
00046 CandShowerHandle::CandShowerHandle(const CandShowerHandle &cdh) :
00047   CandRecoHandle(cdh)
00048 {
00049 }
00050 
00051 //______________________________________________________________________
00052 CandShowerHandle::CandShowerHandle(CandShower *cd) :
00053   CandRecoHandle(cd)
00054 {
00055 }
00056 
00057 //______________________________________________________________________
00058 CandShowerHandle::~CandShowerHandle()
00059 {
00060 }
00061 
00062 //______________________________________________________________________
00063 CandShowerHandle *CandShowerHandle::DupHandle() const
00064 {
00065    return (new CandShowerHandle(*this));
00066 }
00067 
00068 //______________________________________________________________________
00069 void CandShowerHandle::Trace(const char *c) const
00070 {
00071   MSG("Cand", Msg::kDebug)
00072     << "**********Begin CandShowerHandle::Trace(\"" << c << "\")" << endl
00073            << "Information from CandShowerHandle's CandHandle: " << endl;
00074   CandHandle::Trace(c);
00075   MSG("Cand", Msg::kDebug)
00076      << "**********End CandShowerHandle::Trace(\"" << c << "\")" << endl;
00077 }
00078 
00079 
00080 
00081 
00082 
00083 
00084 //______________________________________________________________________
00085 void CandShowerHandle::SetU(Int_t plane, Float_t tpos) 
00086 {
00087   CandShower *shower = dynamic_cast<CandShower*>(GetOwnedCandBase());
00088   shower->fUPos[plane] = tpos;
00089 }
00090 
00091 Float_t CandShowerHandle::GetU(Int_t plane) const
00092 {
00093   const CandShower *shower = dynamic_cast<const CandShower*>(GetCandBase());
00094   if ((shower->fUPos).count(plane)) return shower->fUPos[plane];
00095   return -99999.;
00096 }
00097 
00098 //______________________________________________________________________
00099 
00100 void CandShowerHandle::SetV(Int_t plane, Float_t tpos) 
00101 {
00102   CandShower *shower = dynamic_cast<CandShower*>(GetOwnedCandBase());
00103   shower->fVPos[plane] = tpos;
00104 }
00105 
00106 Float_t CandShowerHandle::GetV(Int_t plane) const
00107 {
00108   const CandShower *shower = dynamic_cast<const CandShower*>(GetCandBase());
00109   if ((shower->fVPos).count(plane)) return shower->fVPos[plane];
00110   return -99999.;
00111 }
00112 
00113 //______________________________________________________________________
00114 
00115 
00116 void CandShowerHandle::SetMinStripPE(Double_t minstrippe) 
00117 {
00118   CandShower *shower = dynamic_cast<CandShower*>(GetOwnedCandBase());
00119   shower->fMinStripPE=minstrippe;
00120 }
00121 
00122 Double_t CandShowerHandle::GetMinStripPE() const
00123 {
00124   const CandShower *shower = dynamic_cast<const CandShower*>(GetCandBase());
00125   return shower->fMinStripPE;
00126 }
00127 
00128 //______________________________________________________________________
00129 
00130 Float_t CandShowerHandle::GetZ(Int_t plane) const
00131 {
00132   TIter stripItr(GetDaughterIterator());
00133   while (CandStripHandle *strip =
00134                            dynamic_cast<CandStripHandle*>(stripItr())) {
00135     if (strip->GetPlane()==plane) return strip->GetZPos();
00136   }
00137   return -1.;
00138 }
00139 
00140 //______________________________________________________________________
00141 
00142 void CandShowerHandle::SetT(Int_t plane, StripEnd::StripEnd_t stripend_t, Double_t time)
00143 {
00144   CandShower *shower = dynamic_cast<CandShower*>(GetOwnedCandBase());
00145   switch (stripend_t) {
00146   case StripEnd::kNegative:
00147     shower->fTime[0][plane] = time;
00148     break;
00149   case StripEnd::kPositive:
00150     shower->fTime[1][plane] = time;
00151     break;
00152   default:
00153     shower->fTime[0][plane] = time;
00154     break;
00155   }
00156 }
00157 
00158 //______________________________________________________________________
00159 
00160 void CandShowerHandle::ClearUVT()
00161 {
00162   CandShower *shower = dynamic_cast<CandShower*>(GetOwnedCandBase());
00163   shower->fUPos.clear();
00164   shower->fVPos.clear();
00165   shower->fTime[0].clear();
00166   shower->fTime[1].clear();
00167 }
00168 
00169 //_____________________________________________________________________
00170 
00171 Float_t CandShowerHandle::GetMinU(Int_t iplane, Double_t minPE) const{
00172   CandStripHandleItr sItr(GetDaughterIterator());
00173   CandStripHandle * begstrip = sItr();
00174   if(!begstrip)return -999;
00175   const VldContext * vld = begstrip->GetVldContext();
00176   if(!vld)return -999;
00177   PlexPlaneId plnid(vld->GetDetector(),iplane,false);
00178   if(plnid.GetPlaneView()!=PlaneView::kU)return -999.;
00179   CandStripHandleItr stripItr(GetDaughterIterator());
00180   CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00181   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00182   stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00183   stripKf = 0;
00184   stripItr.GetSet()->Slice(iplane);
00185 
00186   Float_t minU=999.;
00187   while (const CandStripHandle *strip =
00188          dynamic_cast<const CandStripHandle*>(stripItr())) {
00189     
00190     if(iplane<GetBegPlane() || iplane > GetEndPlane()) return 0.;
00191     
00192     if(strip->GetTPos() < minU && strip->GetCharge(CalDigitType::kPE) > minPE )
00193        minU=strip->GetTPos();
00194   }
00195   if(minU==999)minU=0;
00196   return minU;
00197 }
00198 //_____________________________________________________________________
00199 
00200 Float_t CandShowerHandle::GetMinV(Int_t iplane, Double_t minPE) const{
00201   CandStripHandleItr sItr(GetDaughterIterator());
00202   CandStripHandle * begstrip = sItr();
00203   if(!begstrip)return -999;
00204   const VldContext * vld = begstrip->GetVldContext();
00205   if(!vld)return -999;
00206   PlexPlaneId plnid(vld->GetDetector(),iplane,false);
00207   if(plnid.GetPlaneView()!=PlaneView::kV)return -999.;
00208   if(iplane<GetBegPlane() || iplane > GetEndPlane()) return 0.;
00209   CandStripHandleItr stripItr(GetDaughterIterator());
00210   CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00211   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00212   stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00213   stripKf = 0;
00214   stripItr.GetSet()->Slice(iplane);
00215  
00216   Float_t minV=999.;
00217   while (const CandStripHandle *strip =
00218          dynamic_cast<const CandStripHandle*>(stripItr())) {
00219     if(strip->GetTPos()<minV && strip->GetCharge(CalDigitType::kPE) > minPE )
00220        minV=strip->GetTPos();
00221   }
00222   if(minV==999.)minV=0;
00223   return minV;
00224 }
00225 
00226 //_____________________________________________________________________
00227 
00228 Float_t CandShowerHandle::GetMaxU(Int_t iplane, Double_t minPE) const{
00229   CandStripHandleItr sItr(GetDaughterIterator());
00230   CandStripHandle * begstrip = sItr();
00231   if(!begstrip)return 999;
00232   const VldContext * vld = begstrip->GetVldContext();
00233   if(!vld)return 999;
00234   PlexPlaneId plnid(vld->GetDetector(),iplane,false);
00235   if(plnid.GetPlaneView()!=PlaneView::kU)return 999.;
00236   if(iplane<GetBegPlane() || iplane > GetEndPlane()) return 0.;
00237 
00238   CandStripHandleItr stripItr(GetDaughterIterator());
00239   CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00240   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00241   stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00242   stripKf = 0;
00243   stripItr.GetSet()->Slice(iplane);
00244  
00245   Float_t maxU=-999.;
00246   while (const CandStripHandle *strip =
00247          dynamic_cast<const CandStripHandle*>(stripItr())) {
00248     if(strip->GetTPos() > maxU&& strip->GetCharge(CalDigitType::kPE) > minPE ) maxU=strip->GetTPos();
00249   }
00250   if(maxU==-999.)maxU=0;
00251   return maxU;
00252 }
00253 
00254 //_____________________________________________________________________
00255 
00256 Float_t CandShowerHandle::GetMaxV(Int_t iplane, Double_t minPE) const {
00257   CandStripHandleItr sItr(GetDaughterIterator());
00258   CandStripHandle * begstrip = sItr();
00259   if(!begstrip)return 999;
00260   const VldContext * vld = begstrip->GetVldContext();
00261   if(!vld)return 999;
00262   PlexPlaneId plnid(vld->GetDetector(),iplane,false);
00263   if(plnid.GetPlaneView()!=PlaneView::kV)return 999.;
00264   if(iplane<GetBegPlane() || iplane > GetEndPlane()) return 0.;
00265 
00266   CandStripHandleItr stripItr(GetDaughterIterator());
00267   CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00268   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00269   stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00270   stripKf = 0;
00271   stripItr.GetSet()->Slice(iplane);
00272   Float_t maxV=-999.;
00273   while (const CandStripHandle *strip =
00274          dynamic_cast<const CandStripHandle*>(stripItr())) {
00275     if(strip->GetTPos()>maxV && strip->GetCharge(CalDigitType::kPE)  > minPE)
00276        maxV=strip->GetTPos();
00277   }
00278   if(maxV==-999)maxV=0;
00279   return maxV;
00280 }
00281 //_____________________________________________________________________
00282 
00283 Int_t CandShowerHandle::GetNStrips(Int_t iplane){
00284   Int_t nstrips=0;
00285   CandStripHandleItr stripItr(GetDaughterIterator());
00286   while (const CandStripHandle *strip =
00287          dynamic_cast<const CandStripHandle*>(stripItr())) {
00288     if(strip->GetPlane()==iplane)nstrips++;
00289   }
00290   return nstrips;
00291 }
00292 
00293 //_____________________________________________________________________
00294 
00295 Double_t CandShowerHandle::GetT(Int_t plane,StripEnd::StripEnd_t stripend_t) const
00296 {
00297   const CandShower *shower = dynamic_cast<const CandShower*>(GetCandBase());
00298   if (stripend_t==StripEnd::kNegative) {
00299     if ((shower->fTime[0]).count(plane)) {
00300       return shower->fTime[0][plane];
00301     }
00302     return -99999.;
00303   }
00304   if (stripend_t==StripEnd::kPositive) {
00305     if ((shower->fTime[1]).count(plane)) {
00306       return shower->fTime[1][plane];
00307     }
00308     return -99999.;
00309   }
00310   if ((shower->fTime[0]).count(plane) && (shower->fTime[1]).count(plane)) {
00311     return min(shower->fTime[0][plane],shower->fTime[1][plane]);
00312   }
00313   else if ((shower->fTime[0]).count(plane)) {
00314     return shower->fTime[0][plane];
00315   }
00316   else if ((shower->fTime[1]).count(plane)) {
00317     return shower->fTime[1][plane];
00318   }
00319   return -99999.;
00320 }
00321 
00322 Double_t CandShowerHandle::GetT(StripEnd::StripEnd_t stripend_t,Int_t plane) const
00323 {
00324   return this->GetT(plane,stripend_t);
00325 }
00326 
00327 Double_t CandShowerHandle::GetT(Int_t plane) const
00328 {
00329   return this->GetT(plane,StripEnd::kWhole);
00330 }
00331 
00332 //__________________________________________________________________
00333 Bool_t CandShowerHandle::IsTPosValid(Int_t plane) const
00334 {
00335   const CandShower *shower = dynamic_cast<const CandShower*>(GetCandBase());
00336   if ((shower->fUPos).count(plane) && (shower->fVPos).count(plane)) {
00337     return kTRUE;
00338   }
00339   return kFALSE;
00340 }
00341 //_____________________________________________________________________
00342 Bool_t CandShowerHandle::BelongsWithTrack(CandTrackHandle * trk, 
00343                                           AlgConfig & ac, 
00344                                           const VldContext * vldcptr, 
00345                                           Double_t tolTPos2, Double_t tolZPos, Double_t tolTime){
00346 
00347   if(!trk)return false;
00348 
00349   Int_t pSMPlaneLast = ac.GetInt("SMPlaneLast");
00350   Int_t pSMPlaneFirst = ac.GetInt("SMPlaneFirst");
00351   VldContext vldc = *vldcptr;
00352   UgliGeomHandle ugh(vldc); 
00353   Double_t zGapSM=0;
00354   // "Fixed": MAK, Feb 8, 2005
00355   // GetScintPlnHandle was being called in the ND for planes which didn't
00356   // have any scint. This fix only stops the block below from being called
00357   // in the near detector.  It doesn't require that 
00358   // pSMPlaneLast, pSMPlaneFirst actually refer to planes with scintillator!!
00359   if(vldc.GetDetector() == Detector::kFar){  
00360      // calculate Z gap between SM for later use
00361      PlexPlaneId scintlastid(vldc.GetDetector(),pSMPlaneLast,kFALSE);
00362      PlexPlaneId scintfirstid(vldc.GetDetector(),pSMPlaneFirst,kFALSE);
00363      UgliScintPlnHandle scintlast = ugh.GetScintPlnHandle(scintlastid);
00364      UgliScintPlnHandle scintfirst = ugh.GetScintPlnHandle(scintfirstid);
00365      if (scintlast.IsValid() && scintfirst.IsValid()) {  
00366         zGapSM=scintfirst.GetZ0()-scintlast.GetZ0()-0.0594;
00367      }
00368   }
00369 
00370   MSG("RecoBase",Msg::kDebug) 
00371     << " comparing shower at " 
00372     << this->GetVtxPlane() 
00373     << " " << this->GetVtxU() 
00374     << " " << this->GetVtxV() 
00375     << " with track at " << trk->GetVtxPlane() 
00376     << " " << trk->GetVtxU() 
00377     << " " << trk->GetVtxV() << endl;
00378   Float_t tolTPos = TMath::Sqrt(tolTPos2); 
00379   Int_t fwdTrkPlane=min(trk->GetVtxPlane(),trk->GetEndPlane());
00380   Int_t bckTrkPlane=max(trk->GetVtxPlane(),trk->GetEndPlane());
00381   Int_t fwdShwPlane=min(this->GetVtxPlane(),this->GetEndPlane());
00382   Int_t bckShwPlane=max(this->GetVtxPlane(),this->GetEndPlane());
00383 
00384   Double_t dt = trk->GetVtxT()-this->GetVtxT();
00385   //  if(trk->IsTPosValid(this->GetVtxPlane())){
00386   //   dt = this->GetVtxT()-trk->GetT(this->GetVtxPlane());
00387   //   }
00388 
00389   if( bckShwPlane>=fwdTrkPlane && fwdShwPlane<=bckTrkPlane && fabs(dt)<tolTime){   // trk overlaps in Z and Time
00390     MSG("RecoBase",Msg::kDebug) << " trk and shower overlap in Z and time " << endl;
00391     Bool_t matchu=false;
00392     Bool_t matchv=false;
00393     for(Int_t iloop=0;iloop<2;iloop++){      
00394       if(iloop >0 && (matchu || matchv))tolTPos=tolTPos*2;
00395       // if on second loop and no matches so far, give up
00396       if(iloop >0 && !matchu && !matchv) break;
00397       for (Int_t iplane=fwdShwPlane; iplane<=bckShwPlane;iplane++){  // look for trk passing through shower
00398         
00399         Double_t trkU(0.),trkV(0.);
00400         PlexPlaneId plnid(GetVldContext()->GetDetector(),iplane,false);
00401         if(plnid.GetPlaneView()==PlaneView::kU){
00402           trkU=trk->GetU(iplane);
00403           if(!trk->IsTPosValid(iplane)){
00404             if(abs(iplane-trk->GetVtxPlane())<abs(iplane-trk->GetEndPlane())){
00405               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00406               Double_t dz=scintpln.GetZ0()-trk->GetVtxZ();
00407               trkU = trk->GetVtxU() + dz*trk->GetVtxDirCosU()/trk->GetVtxDirCosZ();
00408             }
00409             else{
00410               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00411               Double_t dz=scintpln.GetZ0()-trk->GetEndZ();
00412               trkU = (trk->GetEndU() + dz*trk->GetEndDirCosU()/trk->GetEndDirCosZ());
00413             }
00414           }
00415           if(trkU>=this->GetMinU(iplane)-tolTPos && trkU<=this->GetMaxU(iplane)+tolTPos  && this->GetNStrips(iplane)>0){
00416             MSG("RecoBase", Msg::kDebug)<< " u match! (details follow)" << endl;
00417             matchu=true;
00418           }
00419         }
00420         else if(plnid.GetPlaneView()==PlaneView::kV){
00421           trkV=trk->GetV(iplane);
00422           if(!trk->IsTPosValid(iplane)){
00423             if(abs(iplane-trk->GetVtxPlane())<abs(iplane-trk->GetEndPlane())){
00424               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00425               Double_t dz=scintpln.GetZ0()-trk->GetVtxZ();
00426               trkV = (trk->GetVtxV() + dz*trk->GetDirCosV()/trk->GetDirCosZ());
00427             }
00428             else{
00429               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00430               Double_t dz=scintpln.GetZ0()-trk->GetEndZ();
00431               trkV= (trk->GetEndV() + dz*trk->GetEndDirCosV()/trk->GetEndDirCosZ());
00432             }
00433           }
00434           if(trkV>=this->GetMinV(iplane)-tolTPos && trkV<=this->GetMaxV(iplane)+tolTPos && this->GetNStrips(iplane)>0){
00435             MSG("RecoBase", Msg::kDebug)<< " v match! (details follow)" << endl;
00436             matchv=true;
00437           }
00438         }
00439         MSG("RecoBase",Msg::kDebug) << " plane  " << iplane 
00440                                     << " trk u " << trkU 
00441                                     << " min/max U " << this->GetMinU(iplane) 
00442                                     << "/" << this->GetMaxU(iplane) 
00443                                     << " trk v " << trkV 
00444                                     << " min/max V " << this->GetMinV(iplane) 
00445                                     << "/" << this->GetMaxV(iplane) << endl; 
00446         if(matchu && matchv) {
00447           MSG("RecoBase", Msg::kDebug)<< " 3D  match! " << endl;
00448           return true;
00449         }
00450       }
00451     }
00452   }
00453   if( this->GetVtxZ()>min(trk->GetVtxZ(),trk->GetEndZ())-tolZPos  &&
00454       this->GetVtxZ()<max(trk->GetVtxZ(),trk->GetEndZ())+tolZPos) {
00455     Double_t du= trk->GetVtxU()-this->GetVtxU();
00456     Double_t dv= trk->GetVtxV()-this->GetVtxV();        
00457     if(trk->IsTPosValid(this->GetVtxPlane())){
00458       MSG("RecoBase",Msg::kDebug) << " TPos is valid @ shower vtx " 
00459                                   << " track u/v " << trk->GetU(this->GetVtxPlane()) 
00460                                   << " " << trk->GetV(this->GetVtxPlane()) << endl;
00461       
00462       du = this->GetVtxU()-trk->GetU(this->GetVtxPlane());
00463       dv = this->GetVtxV()-trk->GetV(this->GetVtxPlane());
00464   }
00465     MSG("RecoBase",Msg::kDebug)
00466       << "    at plane " << this->GetVtxPlane() << " dt2  " 
00467       << du*du+dv*dv << "/" << tolTPos2
00468       << " dt " << dt*1.e9 
00469       << "/" << tolTime*1e9 <<  "\n";
00470     if(dt<tolTime && 
00471        du*du+dv*dv<tolTPos2) return true;   // check whether vertex within minimum distance from trk 
00472   }
00473   // finally, check distance between vertices 
00474   Double_t trkZ=trk->GetVtxZ();
00475   Double_t trkU=trk->GetVtxU();
00476   Double_t trkV=trk->GetVtxV();
00477   Double_t trkT=trk->GetVtxT();
00478   Double_t trkDirU=trk->GetVtxDirCosU();
00479   Double_t trkDirV=trk->GetVtxDirCosV();
00480   Double_t trkDirZ=trk->GetVtxDirCosZ();
00481   Int_t trkPlane=trk->GetVtxPlane();
00482   if(fabs(this->GetVtxZ()-trk->GetEndZ())<fabs(this->GetVtxZ()-trk->GetVtxZ())){
00483     trkZ=trk->GetEndZ();
00484     trkU=trk->GetEndU();
00485     trkV=trk->GetEndV();
00486     trkT=trk->GetEndT();
00487     trkPlane=trk->GetEndPlane();
00488     trkDirU=trk->GetEndDirCosU();
00489     trkDirV=trk->GetEndDirCosV();
00490     trkDirZ=trk->GetEndDirCosZ();
00491   }
00492   Double_t dz = this->GetVtxZ()-trkZ;
00493   dt = this->GetVtxT()-trkT;
00494   //compensate dz for SM gap if necessary
00495   if (vldcptr->GetDetector()==Detector::kFar &&
00496       this->GetVtxPlane()<=pSMPlaneLast &&
00497       trkPlane>=pSMPlaneFirst) dz+=zGapSM;
00498   else if (vldcptr->GetDetector()==Detector::kFar &&
00499            this->GetVtxPlane()>=pSMPlaneLast &&
00500            trkPlane<=pSMPlaneFirst) dz-=zGapSM;
00501   
00502   Double_t du = (trkU + dz*trkDirU/trkDirZ) - this->GetVtxU();
00503   Double_t dv = (trkV + dz*trkDirV/trkDirZ) - this->GetVtxV();
00504 
00505   MSG("RecoBase",Msg::kDebug)
00506     << "    dvertex shower/track " << du
00507     << " " << dv << " " << dz <<" "  << dt*1.e9 << "\n";
00508   if (du*du+dv*dv<tolTPos2 &&
00509       fabs(dz)<tolZPos && 
00510       fabs(dt)<tolTime) return true;
00511 
00512   dz = this->GetEndZ()-trkZ;
00513   dt = this->GetVtxT()-trkT;
00514   //compensate dz for SM gap if necessary
00515   if (vldcptr->GetDetector()==Detector::kFar &&
00516       this->GetVtxPlane()<=pSMPlaneLast &&
00517       trkPlane>=pSMPlaneFirst) dz+=zGapSM;
00518   else if (vldcptr->GetDetector()==Detector::kFar &&
00519            this->GetVtxPlane()>=pSMPlaneLast &&
00520            trkPlane<=pSMPlaneFirst) dz-=zGapSM;
00521   
00522   du = (trkU + dz*trkDirU/trkDirZ) - this->GetEndU();
00523   dv = (trkV + dz*trkDirV/trkDirZ) - this->GetEndV();
00524 
00525   MSG("RecoBase",Msg::kDebug)
00526     << "    dvertex shower/track " << du
00527     << " " << dv << " " << dz <<" "  << dt*1.e9 << "\n";
00528   if (du*du+dv*dv<tolTPos2 &&
00529       fabs(dz)<tolZPos && 
00530       fabs(dt)<tolTime) return true;
00531 
00532   return false;
00533 }
00534 //_____________________________________________________________________
00535 Bool_t CandShowerHandle::BelongsWithShower(CandShowerHandle * shw, 
00536                                           AlgConfig & ac, 
00537                                           const VldContext * vldcptr, 
00538                                           Double_t tolTPos2, Double_t tolZPos, Double_t tolTime){
00539 
00540   if(!shw)return false;
00541 
00542   Int_t pSMPlaneLast = ac.GetInt("SMPlaneLast");
00543   Int_t pSMPlaneFirst = ac.GetInt("SMPlaneFirst");
00544   VldContext vldc = *vldcptr;
00545   UgliGeomHandle ugh(vldc); 
00546   Double_t zGapSM=0;
00547 
00548   // "Fixed": MAK, Feb 8, 2005
00549   // GetScintPlnHandle was being called in the ND for planes which didn't
00550   // have any scint. This fix only stops the block below from being called
00551   // in the near detector.  It doesn't require that 
00552   // pSMPlaneLast, pSMPlaneFirst actually refer to planes with scintillator!!
00553   if(vldc.GetDetector() == Detector::kFar){
00554      // calculate Z gap between SM for later use
00555      PlexPlaneId scintlastid(vldc.GetDetector(),pSMPlaneLast,kFALSE);
00556      PlexPlaneId scintfirstid(vldc.GetDetector(),pSMPlaneFirst,kFALSE);
00557      UgliScintPlnHandle scintlast = ugh.GetScintPlnHandle(scintlastid);
00558      UgliScintPlnHandle scintfirst = ugh.GetScintPlnHandle(scintfirstid);
00559      if (scintlast.IsValid() && scintfirst.IsValid()) {  
00560         zGapSM=scintfirst.GetZ0()-scintlast.GetZ0()-0.0594;
00561      } 
00562   }
00563 
00564   MSG("RecoBase",Msg::kDebug) << " shower extents " << this->GetVtxZ() << " " << this->GetEndZ() << " " << shw->GetVtxZ() << " " << shw->GetEndZ() << endl;
00565   Double_t dz = this->GetVtxZ()-shw->GetVtxZ();
00566   Double_t dzend= this->GetEndZ()-shw->GetVtxZ();
00567   Double_t dzend2= this->GetVtxZ()-shw->GetEndZ();
00568   Double_t dzend3= this->GetEndZ()-shw->GetEndZ();
00569   Double_t du= this->GetVtxU()-shw->GetVtxU();
00570   Double_t dv= this->GetVtxV()-shw->GetVtxV();  
00571   Double_t dt= this->GetVtxT()-shw->GetVtxT();
00572  
00573   //compensate dz for SM gap if necessary
00574   if (vldcptr->GetDetector()==Detector::kFar &&
00575       this->GetVtxPlane()<=pSMPlaneLast &&
00576       shw->GetVtxPlane()>=pSMPlaneFirst) dz+=zGapSM;
00577   else if (vldcptr->GetDetector()==Detector::kFar &&
00578            this->GetVtxPlane()>=pSMPlaneLast &&
00579            shw->GetVtxPlane()<=pSMPlaneFirst) dz-=zGapSM;
00580   
00581   if (vldcptr->GetDetector()==Detector::kFar &&
00582       this->GetEndPlane()<=pSMPlaneLast &&
00583       shw->GetVtxPlane()>=pSMPlaneFirst) dzend+=zGapSM;
00584   else if (vldcptr->GetDetector()==Detector::kFar &&
00585            this->GetEndPlane()>=pSMPlaneLast &&
00586            shw->GetVtxPlane()<=pSMPlaneFirst) dzend-=zGapSM;
00587   if(fabs(dz)>fabs(dzend)){
00588     du= this->GetEndU()-shw->GetVtxU();
00589     dv= this->GetEndV()-shw->GetVtxV(); 
00590     dz=dzend;
00591   }
00592   if (vldcptr->GetDetector()==Detector::kFar &&
00593       this->GetVtxPlane()<=pSMPlaneLast &&
00594       shw->GetEndPlane()>=pSMPlaneFirst) dzend2+=zGapSM;
00595   else if (vldcptr->GetDetector()==Detector::kFar &&
00596            this->GetVtxPlane()>=pSMPlaneLast &&
00597            shw->GetEndPlane()<=pSMPlaneFirst) dzend2-=zGapSM;
00598   if(fabs(dz)>fabs(dzend2)){
00599     du= this->GetVtxU()-shw->GetEndU();
00600     dv= this->GetVtxV()-shw->GetEndV(); 
00601     dz=dzend2;
00602   }
00603   if (vldcptr->GetDetector()==Detector::kFar &&
00604       this->GetEndPlane()<=pSMPlaneLast &&
00605       shw->GetEndPlane()>=pSMPlaneFirst) dzend3+=zGapSM;
00606   else if (vldcptr->GetDetector()==Detector::kFar &&
00607            this->GetEndPlane()>=pSMPlaneLast &&
00608            shw->GetEndPlane()<=pSMPlaneFirst) dzend3-=zGapSM;
00609   if(fabs(dz)>fabs(dzend3)){
00610     du= this->GetEndU()-shw->GetEndU();
00611     dv= this->GetEndV()-shw->GetEndV(); 
00612     dz=dzend3;
00613   }
00614 
00615  Double_t tolTp=tolTPos2;
00616   Double_t tolZ=tolZPos;
00617   if(this->GetEnergy()>15.0 || shw->GetEnergy()>15.0){
00618     tolTp=tolTPos2*4;
00619   }
00620   MSG("RecoBase",Msg::kDebug)
00621     << "    dvertex shower/shower " << du
00622     << " " << dv << " " << dz <<" "  << dt*1.e9 << "\n";
00623   if(( (du*du+dv*dv<tolTp &&fabs(dz)<tolZ) ||  
00624        (du*du<tolTp/4. && dv*dv<tolTp*4. && fabs(dz)<tolZ/2.) ||
00625        (du*du<tolTp*4 && dv*dv<tolTp/4. && fabs(dz)<tolZ/2.))
00626      &&  fabs(dt)<tolTime) return true;
00627 
00628   return false;
00629 }
00630 
00631 //______________________________________________________________________
00632 void CandShowerHandle::AddCluster(CandClusterHandle *cluster)
00633 {
00634   const TObjArray &clusterlist =
00635                (dynamic_cast<CandShower*>(GetOwnedCandBase()))->fClusterList;
00636   Bool_t found(0);
00637   for (Int_t i=0; i<=clusterlist.GetLast() && !found; i++) {
00638     CandClusterHandle *target =
00639                     dynamic_cast<CandClusterHandle*>(clusterlist.At(i));
00640     if (*cluster == *target) found = 1;
00641   }
00642   if (!found) {                 // don't want to duplicate object in list
00643     CandClusterHandle * cl=cluster->DupHandle();
00644     (dynamic_cast<CandShower*>(GetOwnedCandBase()))->fClusterList.
00645       Add(cl);
00646   }
00647   return;
00648 }
00649 //______________________________________________________________________
00650 void CandShowerHandle::RemoveCluster(CandClusterHandle *cluster)
00651 {
00652   const TObjArray &clusterlist =
00653                (dynamic_cast<CandShower*>(GetOwnedCandBase()))->fClusterList;
00654   Bool_t found(0);
00655   for (Int_t i=0; i<=clusterlist.GetLast() && !found; i++) {
00656     CandClusterHandle *target =
00657                     dynamic_cast<CandClusterHandle*>(clusterlist.At(i));
00658     if (*cluster == *target){
00659       (dynamic_cast<CandShower*>(GetOwnedCandBase()))->fClusterList.
00660         RemoveAt(i);
00661       delete target;
00662       (dynamic_cast<CandShower*>(GetOwnedCandBase()))->fClusterList.
00663                                               Compress();
00664      return;
00665     }
00666   }
00667   return;
00668 }
00669 
00670 //_____________________________________________________________________
00671 bool CandShowerHandle::IsContained(){
00672   // determine whether shower is contained.
00673 
00674   CandStripHandleItr sItr(GetDaughterIterator());
00675   CandStripHandle * begstrip = sItr();
00676   if(!begstrip)return true;
00677   const VldContext * vld = begstrip->GetVldContext();
00678   if(!vld)return true;
00679 
00680   Bool_t contained =true;
00681   double totcharge = 0;
00682   double containedcharge = 0;
00683   UgliGeomHandle ugh(*vld); 
00684   PlaneOutline pl;
00685   float detzmin=0;
00686   float detzmax=999;
00687   ugh.GetZExtent(detzmin,detzmax,-1);
00688   float u=0,v=0;
00689   float xedge,yedge,dist;
00690 
00691   while(CandStripHandle * strip = sItr()){
00692     float z = strip->GetZPos();
00693     if(z<detzmin+0.1*Munits::m ||
00694        z>detzmax-0.1*Munits::m) contained=false;
00695     PlexPlaneId plnid(vld->GetDetector(),strip->GetPlane(),false);
00696     if(plnid.GetPlaneView()==PlaneView::kV){
00697       v=strip->GetTPos();
00698       u=GetU(strip->GetPlane());
00699     }
00700     else if(plnid.GetPlaneView()==PlaneView::kU){
00701       u=strip->GetTPos();
00702       v=GetV(strip->GetPlane());
00703     }
00704     float x = 0.707*(u-v);
00705     float y = 0.707*(u+v);
00706     pl.DistanceToEdge(x, y,
00707                       plnid.GetPlaneView(),
00708                       plnid.GetPlaneCoverage(),
00709                       dist, xedge, yedge);
00710     bool isInside = pl.IsInside( x, y,
00711                                  plnid.GetPlaneView(),
00712                                  plnid.GetPlaneCoverage());
00713     
00714     totcharge += strip->GetCharge();
00715     isInside &= dist>0.1*Munits::m;
00716     if(isInside) {
00717       containedcharge +=strip->GetCharge();
00718     }
00719   }
00720   if(totcharge>0) contained = containedcharge/totcharge>0.9;
00721   return contained;
00722 }
00723 
00724 //______________________________________________________________________
00725 
00726 const CandClusterHandle *CandShowerHandle::GetCluster(Int_t i) const
00727 {
00728   const TObjArray *fClusterList = &(dynamic_cast<const CandShower*>
00729     (GetCandBase())->fClusterList);
00730   if (i>fClusterList->GetLast()) {
00731     return 0;
00732   }
00733   return dynamic_cast<const CandClusterHandle*>(fClusterList->At(i));
00734 }
00735 
00736 //______________________________________________________________________
00737 
00738 Int_t CandShowerHandle::GetLastCluster() const
00739 {
00740   return dynamic_cast<const CandShower*>(GetCandBase())->fClusterList.
00741                                                               GetLast();
00742 }
00743 
00744 //______________________________________________________________________
00745 
00746 Bool_t CandShowerHandle::IsUnphysical(Float_t xtalkFrac,Float_t xtalkCut)
00747 {
00748   if(this->GetNStrip()<=0) return true;
00749   Float_t nxtalk = 0;
00750   CandStripHandleItr stripItr(GetDaughterIterator());
00751   while (const CandStripHandle *strip =
00752          dynamic_cast<const CandStripHandle*>(stripItr())) {
00753     if(strip->GetCharge(CalDigitType::kPE) < xtalkCut) nxtalk+=1;      
00754   }
00755   if(nxtalk/Float_t(this->GetNStrip())>xtalkFrac) return true;
00756   return false;
00757 }
00758 
00759 //______________________________________________________________________
00760 
00761 void CandShowerHandle::SetEnergy(Double_t energy, 
00762                                  CandShowerHandle::ShowerType_t showertype)
00763 {
00764  
00765   switch (showertype){
00766   case  kWtCC:
00767     dynamic_cast<CandShower*>(GetOwnedCandBase())->fEnergy = energy;
00768     dynamic_cast<CandShower*>(GetOwnedCandBase())->fEnergy_wtCC = energy;
00769     break;
00770   case kCC:
00771     dynamic_cast<CandShower*>(GetOwnedCandBase())->fEnergy_CC = energy;
00772     break;
00773   case kWtNC:
00774     dynamic_cast<CandShower*>(GetOwnedCandBase())->fEnergy_wtNC = energy;
00775     break;
00776   case kNC:
00777     dynamic_cast<CandShower*>(GetOwnedCandBase())->fEnergy_NC = energy;
00778     break;
00779   case kEM:
00780     dynamic_cast<CandShower*>(GetOwnedCandBase())->fEnergy_EM = energy;
00781     break;
00782   }
00783 }
00784 
00785 Double_t CandShowerHandle::GetEnergy(CandShowerHandle::ShowerType_t showertype) const
00786 {
00787   switch (showertype){
00788   case kWtCC:
00789     return dynamic_cast<const CandShower*>(GetCandBase())->fEnergy_wtCC;
00790   case kCC:
00791     return dynamic_cast<const CandShower*>(GetCandBase())->fEnergy_CC;
00792   case kWtNC:
00793     return dynamic_cast<const CandShower*>(GetCandBase())->fEnergy_wtNC;
00794   case kNC:
00795     return dynamic_cast<const CandShower*>(GetCandBase())->fEnergy_NC;
00796   case kEM:
00797     return dynamic_cast<const CandShower*>(GetCandBase())->fEnergy_EM;
00798   default:  
00799     MSG("RecoBase",Msg::kWarning) << " requested shower energy of unknown type " << endl;
00800     return 0;
00801   }
00802 }
00803 
00804 void CandShowerHandle::CalibrateEnergy(CandTrackHandle * trk, 
00805                          AlgConfig & ac){
00806 
00807   CandStripHandleItr sItr(GetDaughterIterator());
00808   CandStripHandle * begstrip = sItr();
00809   const VldContext * vld = begstrip->GetVldContext();
00810 
00811   Double_t fCCWtLowScale = ac.GetDouble("CCWtLowScale");
00812   Double_t fCCWtLowC1 = ac.GetDouble("CCWtLowC1");
00813   Double_t fCCWtLowC2 = ac.GetDouble("CCWtLowC2");
00814   Double_t fCCWtLowC3 = ac.GetDouble("CCWtLowC3");
00815   Double_t fCCWtLowC4 = ac.GetDouble("CCWtLowC4");
00816   Double_t fCCWtHiC0 = ac.GetDouble("CCWtHiC0");
00817   Double_t fCCWtHiC1 = ac.GetDouble("CCWtHiC1");
00818 
00819   Double_t fCCLinLowScale = ac.GetDouble("CCLinLowScale");
00820   Double_t fCCLinLowC1 = ac.GetDouble("CCLinLowC1");
00821   Double_t fCCLinLowC2 = ac.GetDouble("CCLinLowC2");
00822   Double_t fCCLinHiC0 = ac.GetDouble("CCLinHiC0");
00823   Double_t fCCLinHiC1 = ac.GetDouble("CCLinHiC1");
00824 
00825   Double_t fNCWtLowScale = ac.GetDouble("NCWtLowScale");
00826   Double_t fNCWtMidScale = ac.GetDouble("NCWtMidScale");
00827   Double_t fNCWtLowC1 = ac.GetDouble("NCWtLowC1");
00828   Double_t fNCWtLowC2 = ac.GetDouble("NCWtLowC2");
00829   Double_t fNCWtLowC3 = ac.GetDouble("NCWtLowC3");
00830   Double_t fNCWtMidC0 = ac.GetDouble("NCWtMidC0");
00831   Double_t fNCWtMidC1 = ac.GetDouble("NCWtMidC1");
00832   Double_t fNCWtMidC2 = ac.GetDouble("NCWtMidC2");
00833   Double_t fNCWtMidC3 = ac.GetDouble("NCWtMidC3");
00834   Double_t fNCWtHiC0 = ac.GetDouble("NCWtHiC0");
00835   Double_t fNCWtHiC1 = ac.GetDouble("NCWtHiC1");
00836 
00837   Double_t fNCLinLowScale = ac.GetDouble("NCLinLowScale");
00838   Double_t fNCLinLowC1 = ac.GetDouble("NCLinLowC1");
00839   Double_t fNCLinLowC2 = ac.GetDouble("NCLinLowC2");
00840   Double_t fNCLinHiC0 = ac.GetDouble("NCLinHiC0");
00841   Double_t fNCLinHiC1 = ac.GetDouble("NCLinHiC1");
00842 
00843   Double_t fEMLowScale = ac.GetDouble("EMLowScale");
00844   Double_t fEMLowC1 = ac.GetDouble("EMLowC1");
00845   Double_t fEMLowC2 = ac.GetDouble("EMLowC2");
00846   Double_t fEMHiC0 = ac.GetDouble("EMHiC0");
00847   Double_t fEMHiC1 = ac.GetDouble("EMHiC1");
00848 
00849   Double_t shw_linmipsum=GetCharge(CalStripType::kMIP);
00850   Double_t shw_wtmipsum=0;
00851   Double_t totshw_wtmipsum=0;
00852   Double_t totshw_linmipsum = GetCharge(CalStripType::kMIP);
00853 
00854   // loop through strips once, getting linear version of shower sum, used
00855   // to obtain deweighting factors as well
00856 
00857   Int_t maxpln = 0;
00858   Int_t minpln = 999;
00859   Int_t shared_planes = 0;
00860   Int_t shared_strips = 0;
00861   Int_t nshwstp = 0;
00862   Double_t dedx_1 = 0.;
00863   Double_t dedx_2 = 0.;
00864   Double_t avg_dedx = 0.;
00865   Double_t reco_emu = 0.;
00866   Double_t reco_dircosz = 0.;
00867   Double_t reco_emu_endshw=0.;
00868 
00869   if (trk) {
00870     reco_dircosz = fabs(trk->GetDirCosZ());
00871     CandFitTrackHandle* fittrk = dynamic_cast<CandFitTrackHandle*>(trk); // if there is a fitted track then use it for improved tracking + direction in the shower
00872     reco_emu = sqrt(trk->GetMomentum()*trk->GetMomentum()+ 0.10566*0.10566);
00873     if(fittrk){
00874       reco_dircosz = fabs(fittrk->GetDirCosZ());
00875       if(fittrk->GetMomentum()>0.) {reco_emu = sqrt(fittrk->GetMomentum()*fittrk->GetMomentum()+ 0.10566*0.10566);}
00876     }
00877    
00878     CandStripHandleItr stripItr(GetDaughterIterator());
00879     // check all strips in shower,and see determine length of track through shower
00880     while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
00881       if (strip) {
00882         nshwstp++;
00883         if((!fittrk && trk->FindDaughter(strip)) || (fittrk && fittrk->FindDaughter(strip))){     
00884           shared_strips++;
00885           if(strip->GetPlane()>maxpln) maxpln = strip->GetPlane();
00886           if(strip->GetPlane()<minpln) minpln = strip->GetPlane(); 
00887         }
00888       }
00889     }
00890     shared_planes = maxpln - minpln + 1;
00891     
00892     reco_emu_endshw = reco_emu-shared_planes/(30.*reco_dircosz);
00893     
00894     if(fittrk && fittrk->GetPlaneQP(maxpln)!=0.){
00895       Double_t pln_qp = fittrk->GetPlaneQP(maxpln)*(1.013*fabs(fittrk->GetPlaneQP(maxpln)));
00896       reco_emu_endshw = sqrt(1./pln_qp*1./pln_qp + 0.10566*0.10566);
00897     }
00898 
00899     if(reco_emu_endshw>reco_emu){reco_emu_endshw = reco_emu-shared_planes/(30.*reco_dircosz);}// get best estimate for muon energy at the end of the shw using fittrk info if present
00900 
00901      dedx_1 = DeDx(reco_emu*1000)/1.985; // get dedx at vtx, normalised to MIPs
00902     if(reco_emu_endshw>0.106){ // stop non-physical situations where more energy is lost than initially present
00903       dedx_2 = DeDx(reco_emu_endshw*1000)/1.985; // get approx dedx at end of shw
00904       avg_dedx = (dedx_1 > dedx_2)? 0.5*(dedx_1+dedx_2):dedx_1; // stop muons which stop in the shw from skewing the dedx factor
00905     }
00906     else{avg_dedx = dedx_1;}
00907     shw_linmipsum -= avg_dedx*shared_strips/reco_dircosz; // remove muon chg from shw
00908   }
00909   else{
00910     shw_linmipsum=totshw_linmipsum;
00911   }
00912 
00913   Double_t CCEnergyLin = (shw_linmipsum<fCCLinLowScale) ? shw_linmipsum*fCCLinLowC1 + shw_linmipsum*shw_linmipsum*fCCLinLowC2: fCCLinHiC0 +  shw_linmipsum*fCCLinHiC1;
00914   if(shw_linmipsum<=0) CCEnergyLin=0;
00915         
00916   Double_t NCEnergyLin = (totshw_linmipsum<fNCLinLowScale) ? totshw_linmipsum*fNCLinLowC1 + totshw_linmipsum*totshw_linmipsum*fNCLinLowC2: fNCLinHiC0 +  totshw_linmipsum*fNCLinHiC1;
00917   if(totshw_linmipsum<=0) NCEnergyLin=0;
00918         
00919   Double_t EMEnergy = (totshw_linmipsum<fEMLowScale) ? totshw_linmipsum*fEMLowC1  + totshw_linmipsum*totshw_linmipsum*fEMLowC2: fEMHiC0 +  totshw_linmipsum*fEMHiC1;
00920   if(totshw_linmipsum<=0) EMEnergy=0;
00921 
00922   // for old software, exponential deweighting function is used.  For SS, use cubic
00923 
00924   Double_t deweightfactorCC = 1.0;
00925   Double_t deweightfactorNC = 1.0;
00926  
00927   Double_t fdeweightLowScale = ac.GetDouble("deweightLowScale");
00928   Double_t fdeweightC0 = ac.GetDouble("deweightC0");
00929   Double_t fdeweightC1 = ac.GetDouble("deweightC1");
00930   Double_t fdeweightC2 = ac.GetDouble("deweightC2");
00931   Double_t fdeweightC3 = ac.GetDouble("deweightC3");
00932   /*  
00933   if(CCEnergyLin < fdeweightLowScale){
00934     deweightfactorCC =  fdeweightC0 + 
00935       fdeweightC1 * CCEnergyLin +
00936       fdeweightC2 * CCEnergyLin * CCEnergyLin +
00937       fdeweightC3 * CCEnergyLin * CCEnergyLin * CCEnergyLin;
00938   }
00939   */
00940   // Find final deweight power
00941   double SigHalfPoint = 8.0;
00942   double SigGrad = 0.5;      
00943   deweightfactorCC = 0.25 + 0.75/(1. + exp(-(CCEnergyLin-SigHalfPoint)*SigGrad));
00944  
00945   if(NCEnergyLin < fdeweightLowScale){
00946     deweightfactorNC =  fdeweightC0 + 
00947       fdeweightC1 * NCEnergyLin +
00948       fdeweightC2 * NCEnergyLin * NCEnergyLin +
00949       fdeweightC3 * NCEnergyLin * NCEnergyLin * NCEnergyLin;
00950   }
00951   
00952   // now that deweighting factors are known, calculated de-weighted sum 
00953     CandStripHandleItr stripItr(GetDaughterIterator());
00954   // check all strips in shower,and see whether they are shared with the track 
00955   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
00956     if (strip) {
00957       totshw_wtmipsum+=pow(GetStripCharge(strip,StripEnd::kWhole),deweightfactorNC);
00958       Double_t trackEloss = 0;  
00959       // if track crosses strip, subtract chg dependent on dedx factor (calc earlier) and direction
00960       if(trk){
00961         CandFitTrackHandle* fittrk = dynamic_cast<CandFitTrackHandle*>(trk);
00962         if((!fittrk && trk->FindDaughter(strip)) || (fittrk && fittrk->FindDaughter(strip))){     
00963           trackEloss=avg_dedx/reco_dircosz;
00964         }
00965       }
00966       if(GetStripCharge(strip,StripEnd::kWhole) - trackEloss > 0){
00967         shw_wtmipsum+=pow(GetStripCharge(strip,StripEnd::kWhole)-trackEloss,deweightfactorCC);
00968       }
00969     }
00970   }
00971   
00972   Double_t CCEnergyWt =  fCCWtHiC0 +  shw_wtmipsum*fCCWtHiC1;
00973   if( shw_wtmipsum<fCCWtLowScale){
00974     CCEnergyWt =  shw_wtmipsum*fCCWtLowC1 + 
00975       shw_wtmipsum*shw_wtmipsum*fCCWtLowC2 + 
00976       shw_wtmipsum*shw_wtmipsum*shw_wtmipsum*fCCWtLowC3 + 
00977       shw_wtmipsum*shw_wtmipsum*shw_wtmipsum*shw_wtmipsum*fCCWtLowC4;
00978   }
00979   if(shw_wtmipsum<=0) CCEnergyWt=0;
00980   
00981   Double_t NCEnergyWt = fNCWtHiC0 + totshw_wtmipsum*fNCWtHiC1;
00982   if(totshw_wtmipsum<fNCWtLowScale){ 
00983     NCEnergyWt = totshw_wtmipsum*fNCWtLowC1 + 
00984       totshw_wtmipsum*totshw_wtmipsum*fNCWtLowC2 + 
00985       totshw_wtmipsum*totshw_wtmipsum*totshw_wtmipsum*fNCWtLowC3; 
00986   }
00987 
00988   if(totshw_wtmipsum>=fNCWtLowScale && totshw_wtmipsum<fNCWtMidScale){
00989     NCEnergyWt =  fNCWtMidC0 +  totshw_wtmipsum*fNCWtMidC1 + totshw_wtmipsum*totshw_wtmipsum*fNCWtMidC2 + totshw_wtmipsum*totshw_wtmipsum*totshw_wtmipsum*fNCWtMidC3;
00990   }
00991 
00992   if(totshw_wtmipsum<=0) NCEnergyWt=0;
00993 
00994 // cout << " CC deweightmips = " << shw_wtmipsum << " NC deweight_mips = " <<totshw_wtmipsum  << endl;
00995 //  cout << " CC deweight E = " << CCEnergyWt << " NC dwewightE = " << NCEnergyWt<< endl;
00996   
00997   // for near detector, compensate for near/far shower completeness difference
00998   // and observed difference between near/far mip scale on muon tracks
00999 
01000  if (vld && vld->GetDetector()==Detector::kNear){
01001     CCEnergyWt *= (.97 + exp(-(CCEnergyWt+11.6)/9.3)); 
01002     CCEnergyLin *= (.99 + exp(-(CCEnergyLin+28.2)/11.8)); 
01003     NCEnergyWt *= (1.01 + exp(-(NCEnergyWt+4.33)/5.)); 
01004     NCEnergyLin *= (1. + exp(-(NCEnergyLin+10.1)/5.)); 
01005     EMEnergy *= 1.06; 
01006  }
01007 
01008  // now apply energy corrections, corresponding to EnergyCorrections::MasakiMay17thCGScaled.
01009   
01010  if(vld){
01011    float etemplin = CCEnergyLin;
01012    CCEnergyLin = EnergyCorrections::ShowerEnergyConversionDogwood(etemplin,*vld);
01013    etemplin = CCEnergyWt;
01014    CCEnergyWt = EnergyCorrections::WeightedShowerEnergyConversionDogwood(etemplin,*vld);
01015  }
01016 
01017  
01018  SetEnergy(CCEnergyWt);
01019  SetEnergy(CCEnergyWt,kWtCC);
01020  SetEnergy(CCEnergyLin,kCC);
01021  SetEnergy(NCEnergyWt,kWtNC);
01022  SetEnergy(NCEnergyLin,kNC);
01023  SetEnergy(EMEnergy,kEM);
01024 }
01025 
01026 //______________________________________________________________________
01027 
01028 Double_t CandShowerHandle::DeDx(Double_t emu){
01029   
01030   Double_t dedx = 0.; 
01031   Double_t mm = 105.658389;
01032   Double_t mm_2 = mm*mm;
01033 
01034   if(emu*emu-mm_2<0.){return 0.;}
01035 
01036   Double_t a_2     = 1./(137.036*137.036);
01037   Double_t pi = 3.141;
01038   Double_t Na = 6.023;
01039   Double_t lamda_2 = 3.8616*3.8616;
01040   Double_t Z_A = 0.5377;
01041   Double_t me  = 0.51099906;
01042   Double_t me_2 = me*me;  
01043   Double_t beta = sqrt(emu*emu-mm_2)/emu;
01044   Double_t beta_2 = beta*beta;
01045   Double_t gamma = emu/105.658389;
01046   Double_t gamma_2 = gamma*gamma;
01047   Double_t I = 68.7e-6;
01048   Double_t I_2 = I*I;
01049   Double_t p = sqrt(emu*emu-mm_2);
01050   Double_t p_2 = p*p;
01051   Double_t Em = 2 * me * p_2 / ( me_2 + mm_2 + 2*me*emu );
01052   Double_t Em_2= Em*Em;
01053   Double_t emu_2 = emu*emu;
01054   Double_t X0 = 0.165;
01055   Double_t X1 = 2.503;
01056   Double_t X = log10(beta*gamma);
01057   Double_t a = 0.165;
01058   Double_t C = -3.3;
01059   Double_t m = 3.222;
01060   Double_t d = 0;
01061   if(X0 < X && X < X1){d = 4.6052 * X + a * pow(X1-X,m) + C;}
01062   if(X > X1){d = 4.6052 * X + C;}
01063 
01064   dedx =  a_2 * 2*pi * Na * lamda_2 * Z_A * (me / beta_2) *( log( 2*me*beta_2*gamma_2*Em/I_2 ) - 2*beta_2 + 0.25*(Em_2/emu_2) - d );
01065 
01066   return 10*dedx;
01067 
01068 }
01069 
01070 
01071 //______________________________________________________________________
01072 // Navigation Helpers
01073 //______________________________________________________________________
01074 
01075 NavKey CandShowerHandle::KeyFromSlice(const CandShowerHandle *reco)
01076 {
01077   if (reco->GetCandSlice()) {
01078     return static_cast<Int_t>(reco->GetCandSlice()->GetUidInt());
01079   }
01080   return 0;
01081 
01082 }
01083 
01084 
01085 
01086 
01087 
01088 XXXITRIMP(CandShowerHandle)