loon: CandShowerHandle.cxx Source File

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