AlgShowerSR Class Reference

#include <AlgShowerSR.h>

Inheritance diagram for AlgShowerSR:
AlgBase AlgReco AlgShowerSS

List of all members.

Public Member Functions

 AlgShowerSR ()
virtual ~AlgShowerSR ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void Trace (const char *c) const

Protected Member Functions

void SetUVT (CandShowerHandle *)
void SetUV (CandShowerHandle *)
void SetT (CandShowerHandle *)
Int_t FindTimingDirection (CandShowerHandle &shwh, AlgConfig &ac, Double_t *fitparm, Double_t &timefitchi2)

Detailed Description

Definition at line 22 of file AlgShowerSR.h.


Constructor & Destructor Documentation

AlgShowerSR::AlgShowerSR (  ) 

Within AlgShowerSR::RunAlg, the vertex location and shower direction is determined. First, the most upstream U and V planes are found. Next, a loop over all shower clusters is used to obtain the total shower charge. The charge-weighted U/V positions of the CandClusters in the plane interval from the most upstream planes to VtxPlaneSpan planes downstream is identified as the shower vertex location. To obtain the shower direction, a linear fit is performed to the charge-weighted transverse position in a given plane vs the Z position of that plane. A weighted linear fit is performed in each viewed, weighted by the total cluster charge for this shower in a given plane. The du/dz and dv/dz slopes obtained in these fits are then used to obtain the direction cosines of the shower axis. Finally, the shower energy is set.

Definition at line 42 of file AlgShowerSR.cxx.

00043 {
00044 }

AlgShowerSR::~AlgShowerSR (  )  [virtual]

Definition at line 47 of file AlgShowerSR.cxx.

00048 {
00049 }


Member Function Documentation

Int_t AlgShowerSR::FindTimingDirection ( CandShowerHandle shwh,
AlgConfig ac,
Double_t *  fitparm,
Double_t &  timefitchi2 
) [protected]

Definition at line 801 of file AlgShowerSR.cxx.

References bfld::AsString(), UgliStripHandle::ClearFiber(), digit(), CandDigitHandle::GetCharge(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), PlexSEIdAltL::GetEnd(), UgliStripHandle::GetHalfLength(), Registry::GetInt(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandStripHandle::GetTime(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), CandStripHandle::GetZPos(), Msg::kDebug, Detector::kNear, StripEnd::kNegative, PlaneView::kU, StripEnd::kUnknown, PlaneView::kV, Msg::kVerbose, MSG, PropagationVelocity::Velocity(), LinearFit::Weighted(), and UgliStripHandle::WlsPigtail().

Referenced by RunAlg().

00805 {
00806    Double_t parmMaxTimeFitRes = ac.GetDouble("MaxTimeFitRes");
00807   Double_t parmMinTimeFitSize = ac.GetInt("MinTimeFitSize");
00808   Double_t parmMaxTimeFitOutlier = ac.GetDouble("MaxTimeFitOutlier");
00809   
00810   MSG("AlgShowerSR", Msg::kVerbose) << "in find timing direction" << endl;
00811   
00812   fitparm[0] = 0;
00813   fitparm[1] = 0;
00814   
00815   CandStripHandle *strip = 0;
00816   CandDigitHandle *digit = 0;
00817   
00818   StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown;
00819   
00820   //get VldContext and UgliGeomHandle
00821   CandStripHandleItr stripItr(shwh.GetDaughterIterator());
00822   CandStripHandle *firststrip = stripItr.Ptr();
00823   
00824   assert(firststrip);
00825   
00826   const VldContext *vldcontext = firststrip->GetVldContext();
00827   assert(vldcontext);
00828   
00829   UgliGeomHandle ugh(*vldcontext);
00830   
00831   //arrays for determining the timing direction - limit to
00832   //1000 entries - really shouldnt need more than that to 
00833   //be able to figure out if the shower is going north or 
00834   //south
00835   Double_t time[1000];
00836   Double_t pathLength[1000];
00837   Double_t weight[1000];
00838   Double_t adcsum[1000];
00839   
00840   for(Int_t i = 0; i<1000; ++i){
00841     time[i] = 0.;
00842     pathLength[i] = 0.;
00843     weight[i] = 0.;
00844     adcsum[i] = 0;
00845   }
00846   
00847   //need to put in a check to make sure you dont have more planes in one
00848   //view than the other
00849   Int_t uBeg = 1000;
00850   Int_t uEnd = 0;
00851   Int_t vBeg = 1000;
00852   Int_t vEnd = 0;
00853    while( (strip = stripItr())){
00854     Int_t plane = strip->GetPlane();
00855     if(strip->GetPlaneView() == PlaneView::kV && plane>vEnd) vEnd=plane;
00856     if(strip->GetPlaneView() == PlaneView::kV && plane<vBeg) vBeg=plane;
00857     if(strip->GetPlaneView() == PlaneView::kU && plane>uEnd) uEnd=plane;
00858     if(strip->GetPlaneView() == PlaneView::kU && plane<uBeg) uBeg=plane;
00859   }
00860   
00861   //get the endpoints  
00862   Int_t showerBeg = TMath::Min(uBeg,vBeg);
00863   Int_t showerEnd = TMath::Max(uEnd,vEnd);
00864   
00865   if(TMath::Abs(uBeg-vBeg)>3.){
00866     showerBeg = TMath::Max(uBeg,vBeg) - 1;
00867   } 
00868   if(TMath::Abs(uEnd-vEnd)>3.){
00869     showerEnd = TMath::Min(uEnd,vEnd) + 1;
00870   }
00871   MSG("AlgShowerSR", Msg::kDebug) << "u end points = " << uBeg
00872                                   << " " << uEnd << " v end points = "
00873                                   << vBeg << " " << vEnd << endl
00874                                   << " shower end points = " << showerBeg
00875                                   << " " << showerEnd << endl;
00876   Double_t halfLength = 0.;
00877   Double_t apos = 0.;
00878   Double_t distFromCenter = 0.;
00879   Double_t fiberDist = 0.;
00880   Int_t nctr = 0;
00881   Double_t maxPathLength = 0;
00882   for(Int_t iplane = showerBeg;iplane <= showerEnd;iplane++){
00883     adcsum[nctr] = 0;
00884     time[nctr] = 0; 
00885     pathLength[nctr] = 0;
00886     CandStripHandleItr shwstripItr(shwh.GetDaughterIterator());
00887     while( (strip = shwstripItr()) && nctr<1000){
00888       if(strip->GetPlane() == iplane){
00889         MSG("AlgShowerSR", Msg::kDebug) << "strip from plane " << iplane
00890                                         << " " << Detector::AsString(vldcontext->GetDetector()) << endl;
00891         
00892         //only look at double ended strips if in the far detector
00893         //and strips from planes that have orthogonal measurements 
00894         if( (strip->GetNDaughters() == 2 || vldcontext->GetDetector() == Detector::kNear)){
00895           
00896           MSG("AlgShowerSR", Msg::kVerbose) << "strip good for timing" << endl;
00897           
00898           //get the iterator over the digits in the strip
00899           CandDigitHandleItr digitItr(strip->GetDaughterIterator());
00900           
00901           while( (digit = digitItr())){     
00902             stripEnd = digit->GetPlexSEIdAltL().GetEnd();
00903             pathLength[nctr] = strip->GetZPos();
00904             maxPathLength = TMath::Max(maxPathLength,pathLength[nctr] );
00905             UgliStripHandle stripHandle = ugh.GetStripHandle(strip->GetStripEndId());
00906             halfLength = stripHandle.GetHalfLength();
00907             apos = 0.;
00908             fiberDist = 0.;
00909             
00910             if(strip->GetPlaneView() == PlaneView::kV) apos = shwh.GetU(iplane);
00911             else if(strip->GetPlaneView() == PlaneView::kU) apos = shwh.GetV(iplane);
00912             distFromCenter = 0.;
00913             if(strip->GetPlaneView() == PlaneView::kV) 
00914               distFromCenter = (stripEnd == StripEnd::kNegative) ? apos : -apos;
00915             else if(strip->GetPlaneView() == PlaneView::kU) 
00916               distFromCenter = (stripEnd == StripEnd::kNegative) ? -apos : apos;
00917             
00918             fiberDist = (halfLength + distFromCenter + stripHandle.ClearFiber(stripEnd) 
00919                          + stripHandle.WlsPigtail(stripEnd));
00920             
00921             time[nctr] += digit->GetCharge() * (strip->GetTime(stripEnd) - fiberDist/PropagationVelocity::Velocity());   
00922             
00923             weight[nctr] = 1.0;
00924             adcsum[nctr] += digit->GetCharge();
00925             Double_t temptime =  time[nctr];
00926             if(adcsum[nctr]!=0)temptime /= adcsum[nctr]; 
00927             MSG("AlgShowerSR", Msg::kDebug) << nctr 
00928                                             << " " << pathLength[nctr]
00929                                             << " " << temptime
00930                                             << " " << strip->GetTime(stripEnd)
00931                                             << " " << weight[nctr]
00932                                             << endl;
00933             
00934           }//end loop over digits in strip
00935         }//end if double ended strip
00936       }//end loop over strips in plane
00937     }
00938     if(adcsum[nctr] > 0){
00939       time[nctr] /= adcsum[nctr];        
00940       nctr++;
00941     } 
00942   }
00943   for(Int_t i=0;i<nctr;i++){
00944     MSG("AlgShowerSR",Msg::kDebug) << "TIMEFIT " << i << " " 
00945                                    <<pathLength[i] << " " 
00946                                    << (time[i]-time[0])*1.e9 << " " << weight[i] << endl;
00947   }
00948   
00949   //loop over all strips in the shower
00950   Double_t efitparm[2];
00951   Double_t maxresid = parmMaxTimeFitRes + 1.;
00952   Double_t resid = 0.;
00953   Int_t ctr = nctr;
00954   Int_t imaxresid = 0;
00955   Int_t nremove = 0;
00956   Bool_t dofit = false;
00957   while(maxresid > parmMaxTimeFitRes
00958         && nctr-nremove > parmMinTimeFitSize 
00959         && (1.*nremove) < parmMaxTimeFitOutlier*nctr
00960         ){
00961     dofit = true;
00962     LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm);
00963     maxresid = 0.;
00964     imaxresid = 0;
00965     for(int i=0; i<ctr; ++i){
00966       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1.e9;
00967       if(weight[i]>0. && TMath::Abs(resid)>maxresid){
00968         maxresid = TMath::Abs(resid);
00969         imaxresid = i;
00970       }
00971     }
00972     
00973     MSG("AlgShowerSR",Msg::kDebug) << "constrained fit, dt/ds slope, offset = " 
00974                                    << fitparm[1] << " " << fitparm[0] << endl;
00975     MSG("AlgShowerSR",Msg::kDebug) << "maximum residual " << maxresid << " " 
00976                                    << pathLength[imaxresid] << " " 
00977                                    << time[imaxresid] << " " 
00978                                    << weight[imaxresid] << endl;
00979     
00980     if(maxresid>parmMaxTimeFitRes){
00981       MSG("AlgShowerSR",Msg::kDebug) << "  removing" << endl;
00982       weight[imaxresid] = 0.;
00983       nremove++;
00984     }
00985     
00986     timefitchi2 = 0.;
00987     for(int i = 0; i < ctr; ++i){
00988       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1e9;
00989       timefitchi2 += weight[i]*resid*resid;
00990     }
00991     MSG("AlgShowerSR",Msg::kDebug) << "chi2 = " << timefitchi2 \
00992                                    << "   n = " << nctr-nremove << endl;
00993   }//end loop to find timing fit and remove outliers
00994   
00995   timefitchi2 = 0.;
00996   if(dofit){
00997     for(int i = 0; i<ctr; ++i){
00998       if(weight[i] > 0.){
00999         resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1.e9;
01000         timefitchi2 += weight[i]*resid*resid;
01001       }
01002     }
01003   }
01004 
01005   MSG("AlgShowerSR",Msg::kDebug) << " TIMEFIT " << fitparm[1] 
01006                                  << " chi^2/ndf = " 
01007                                  << timefitchi2/(1.*(nctr-nremove)) 
01008                                  << " " << TMath::Abs(uBeg-uEnd)
01009                                  << " " << TMath::Abs(vBeg-vEnd) << endl;
01010   
01011   //check the chi^2 for the fit - if it is really bad dont change the
01012   //initial guess at the direction for the event, ie just make sure
01013   //that the slope in time vs pathlength is positive;
01014   if(fitparm[1] < 0. && timefitchi2 >= 10.*(nctr-nremove))
01015     fitparm[1] *= -1.;
01016   return nctr-nremove;
01017 }

void AlgShowerSR::RunAlg ( AlgConfig ac,
CandHandle ch,
CandContext cx 
) [virtual]

Implements AlgBase.

Reimplemented in AlgShowerSS.

Definition at line 52 of file AlgShowerSR.cxx.

References AlgReco::Calibrate(), FindTimingDirection(), CandContext::GetDataIn(), Registry::GetInt(), CandStripHandle::GetPlane(), UgliGeomHandle::GetScintPlnHandle(), CandStripHandle::GetStripEndId(), CandHandle::GetVldContext(), UgliPlnHandle::GetZ0(), Msg::kDebug, Msg::kError, PlaneView::kU, PlaneView::kV, MSG, SetUV(), and LinearFit::Weighted().

00053 {
00054 
00055   assert(ch.InheritsFrom("CandShowerHandle"));
00056   CandShowerHandle &csh = dynamic_cast<CandShowerHandle &>(ch);
00057 
00058   assert(cx.GetDataIn());
00059   assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00060 
00061   // grab parameters
00062 
00063   //  Double_t MIPperGeV = ac.GetDouble("MIPperGeV");
00064   Int_t vtxplanespan = ac.GetInt("VtxPlaneSpan");
00065   Int_t isCosmic = ac.GetInt("IsCosmic");
00066 
00067   const TObjArray *tary =
00068                          dynamic_cast<const TObjArray*>(cx.GetDataIn());
00069 
00070   Int_t ubegplane=0;
00071   Int_t vbegplane=0;
00072   Int_t uendplane=0;
00073   Int_t vendplane=0;
00074   CandStripHandle *begstrip=0;
00075   CandStripHandle *endstrip=0;
00076 
00077   for (Int_t i=0; i<=tary->GetLast(); i++) {
00078     TObject *tobj = tary->At(i);
00079     assert(tobj->InheritsFrom("CandClusterHandle"));
00080     CandClusterHandle *clusterhandle =
00081                                  dynamic_cast<CandClusterHandle*>(tobj);
00082     csh.AddCluster(clusterhandle);
00083     switch (clusterhandle->GetPlaneView()) {
00084     case PlaneView::kU:
00085       if (!ubegplane || clusterhandle->GetBegPlane()<ubegplane) {
00086         ubegplane = clusterhandle->GetBegPlane();
00087       }
00088       if(!uendplane || clusterhandle->GetEndPlane()>uendplane) {
00089         uendplane = clusterhandle->GetEndPlane();
00090       }
00091       break;
00092     case PlaneView::kV:
00093       if (!vbegplane || clusterhandle->GetBegPlane()<vbegplane) {
00094         vbegplane = clusterhandle->GetBegPlane();
00095       }
00096       if(!vendplane || clusterhandle->GetEndPlane()>vendplane) {
00097         vendplane = clusterhandle->GetEndPlane();
00098       }
00099       break;
00100     default:
00101       break;
00102     }
00103     if (!begstrip) {
00104       CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator());
00105       begstrip = stripItr();
00106     }
00107   }
00108 
00109   // bail if these shower has an empty daughter list
00110   if(!begstrip){
00111     MSG("ShowerSR", Msg::kError)
00112                  << "Found empty daughter list in CandShowerSR" << endl;
00113     csh.SetDirCosU(0.);
00114     csh.SetDirCosV(0.);
00115     csh.SetDirCosZ(1.);
00116     csh.SetVtxPlane(-1);
00117     csh.SetVtxZ(-1.);
00118     csh.SetVtxU(0.);
00119     csh.SetVtxV(0.);
00120     csh.SetVtxT(0.);
00121     csh.SetEndZ(-1.);
00122     csh.SetEndU(0.);
00123     csh.SetEndV(0.);
00124     csh.SetEndT(0.);
00125     csh.SetEnergy(0);
00126 
00127     return;
00128   }
00129 
00130   const VldContext * vld = begstrip->GetVldContext();
00131   UgliGeomHandle ugh(*vld);
00132 
00133   begstrip = 0;
00134 
00135   Double_t uvtx=0.;
00136   Double_t vvtx=0.;
00137   Double_t uvtxph=0.;
00138   Double_t vvtxph=0.;
00139   Double_t tvtx=0.;
00140 
00141   Double_t uend=0.;
00142   Double_t vend=0.;
00143   Double_t uendph=0.;
00144   Double_t vendph=0.;
00145   Double_t tend=0.;
00146 
00147   Bool_t firstvtx(1);
00148   Bool_t firstend(1);
00149 
00150   Double_t *uvtxfit = new Double_t[vtxplanespan+1];
00151   Double_t *uvtxphfit = new Double_t[vtxplanespan+1];
00152   Double_t *uvtxzfit = new Double_t[vtxplanespan+1];
00153   Double_t *vvtxfit = new Double_t[vtxplanespan+1];
00154   Double_t *vvtxphfit = new Double_t[vtxplanespan+1];
00155   Double_t *vvtxzfit = new Double_t[vtxplanespan+1];
00156 
00157   Double_t *uendfit = new Double_t[vtxplanespan+1];
00158   Double_t *uendphfit = new Double_t[vtxplanespan+1];
00159   Double_t *uendzfit = new Double_t[vtxplanespan+1];
00160   Double_t *vendfit = new Double_t[vtxplanespan+1];
00161   Double_t *vendphfit = new Double_t[vtxplanespan+1];
00162   Double_t *vendzfit = new Double_t[vtxplanespan+1];
00163 
00164   for (Int_t i=0; i<=vtxplanespan; i++) {
00165     uvtxfit[i] = 0.;
00166     uvtxphfit[i] = 0.;
00167     uvtxzfit[i] = 0.;
00168     vvtxfit[i] = 0.;
00169     vvtxphfit[i] = 0.;
00170     vvtxzfit[i] = 0.;
00171     uendfit[i] = 0.;
00172     uendphfit[i] = 0.;
00173     uendzfit[i] = 0.;
00174     vendfit[i] = 0.;
00175     vendphfit[i] = 0.;
00176     vendzfit[i] = 0.;
00177   }
00178 
00179   // loop over cluster daughters, calculating a charge sum, and
00180   // constructing arrays of charge-weighted u vs z and v vs z for the 
00181   // planes within VtxPlaneSpan of beginning and end of shower
00182 
00183   Double_t charge=0.;
00184   for (Int_t i=0; i<=tary->GetLast(); i++) {
00185     TObject *tobj = tary->At(i);
00186     assert(tobj->InheritsFrom("CandClusterHandle"));
00187     CandClusterHandle *clusterhandle =
00188                                  dynamic_cast<CandClusterHandle*>(tobj);
00189     CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator());
00190     while (CandStripHandle *strip = stripItr()) {
00191       csh.AddDaughterLink(*strip);
00192       charge += strip->GetCharge();
00193       if (!begstrip || strip->GetPlane()<begstrip->GetPlane()) {
00194         begstrip = strip;
00195       }
00196       if (!endstrip || strip->GetPlane()>endstrip->GetPlane()) {
00197         endstrip = strip;
00198       }
00199 
00200       switch (strip->GetPlaneView()) {
00201       case PlaneView::kU:
00202         if (strip->GetPlane()<=ubegplane+vtxplanespan) {
00203           Double_t charge = strip->GetCharge();
00204           Double_t tpos = strip->GetTPos();
00205           uvtx += tpos*charge;
00206           uvtxph += charge;
00207           if (firstvtx || strip->GetBegTime()<tvtx) {
00208             firstvtx = 0;
00209             tvtx = strip->GetBegTime();
00210           }
00211           Int_t dplane = strip->GetPlane()-ubegplane;
00212           if(dplane<=vtxplanespan && dplane>+0){
00213             uvtxfit[dplane] += tpos*charge;
00214             uvtxphfit[dplane] += charge;
00215             UgliScintPlnHandle upid = 
00216               ugh.GetScintPlnHandle(strip->GetStripEndId());
00217             uvtxzfit[dplane] = upid.GetZ0();
00218           }
00219         }
00220         if(strip->GetPlane()>=uendplane-vtxplanespan) {
00221           Double_t charge = strip->GetCharge();
00222           Double_t tpos = strip->GetTPos();
00223           uend += tpos*charge;
00224           uendph += charge;
00225           if (firstend || strip->GetBegTime()<tend) {
00226             firstend = 0;
00227             tend = strip->GetBegTime();
00228           }
00229           Int_t dplane = uendplane-strip->GetPlane();
00230           if(dplane<=vtxplanespan && dplane>=0){
00231             uendfit[dplane] += tpos*charge;
00232             uendphfit[dplane] += charge;
00233             UgliScintPlnHandle upid = 
00234               ugh.GetScintPlnHandle(strip->GetStripEndId());
00235             uendzfit[dplane] = upid.GetZ0();
00236           }
00237         }
00238         break;
00239       case PlaneView::kV:
00240         if (strip->GetPlane()<=vbegplane+vtxplanespan) {
00241           Double_t charge = strip->GetCharge();
00242           Double_t tpos = strip->GetTPos();
00243           vvtx += tpos*charge;
00244           vvtxph += charge;
00245           if (firstvtx || strip->GetBegTime()<tvtx) {
00246             firstvtx = 0;
00247             tvtx = strip->GetBegTime();
00248           }
00249           Int_t dplane = strip->GetPlane()-vbegplane;
00250           if(dplane<=vtxplanespan && dplane>=0){
00251             vvtxfit[dplane] += tpos*charge;
00252             vvtxphfit[dplane] += charge;
00253             UgliScintPlnHandle upid = 
00254               ugh.GetScintPlnHandle(strip->GetStripEndId());
00255             vvtxzfit[dplane] = upid.GetZ0();
00256           }
00257         }
00258         if(strip->GetPlane()>=uendplane-vtxplanespan) {
00259           Double_t charge = strip->GetCharge();
00260           Double_t tpos = strip->GetTPos();
00261           vend += tpos*charge;
00262           vendph += charge;
00263           if (firstend || strip->GetBegTime()<tend) {
00264             firstend = 0;
00265             tend = strip->GetBegTime();
00266           }
00267           Int_t dplane = vendplane-strip->GetPlane();
00268           if(dplane<=vtxplanespan && dplane>=0){
00269             vendfit[dplane] += tpos*charge;
00270             vendphfit[dplane] += charge;
00271             UgliScintPlnHandle upid = 
00272               ugh.GetScintPlnHandle(strip->GetStripEndId());
00273             vendzfit[dplane] = upid.GetZ0();
00274           }
00275         }
00276         break;
00277       default:
00278         break;
00279       }
00280     }
00281   }
00282 
00283   // normalize charge weighted fit points, and remove array entries
00284   // with zero pulse height
00285 
00286   Int_t iuvtxfit=0;
00287   Int_t ivvtxfit=0;
00288   Int_t iuendfit=0;
00289   Int_t ivendfit=0;
00290 
00291   for (Int_t i=0; i<=vtxplanespan; i++) {
00292     if (uvtxphfit[i]>0.) {
00293       uvtxfit[iuvtxfit] = uvtxfit[i]/uvtxphfit[i];
00294       uvtxzfit[iuvtxfit]=uvtxzfit[i];
00295       uvtxphfit[iuvtxfit]=uvtxphfit[i];
00296 
00297       iuvtxfit++;
00298     }
00299     if (vvtxphfit[i]>0.) {
00300       vvtxfit[ivvtxfit] = vvtxfit[i]/vvtxphfit[i];
00301       vvtxzfit[ivvtxfit]=vvtxzfit[i];
00302       vvtxphfit[ivvtxfit]=vvtxphfit[i];
00303       ivvtxfit++;
00304     }
00305    if (uendphfit[i]>0.) {
00306       uendfit[iuendfit] = uendfit[i]/uendphfit[i];
00307       uendzfit[iuendfit]=uendzfit[i];
00308       uendphfit[iuendfit]=uendphfit[i];
00309       iuendfit++;
00310     }
00311     if (vendphfit[i]>0.) {
00312       vendfit[ivendfit] = vendfit[i]/vendphfit[i];
00313       vendzfit[ivendfit]=vendzfit[i];
00314       vendphfit[ivendfit]=vendphfit[i];
00315       ivendfit++;
00316     }
00317   }
00318 
00319   //perform weighted fit of u, v vs z positions
00320   Double_t uvtxparm[2];
00321   Double_t vvtxparm[2];
00322   Double_t dudz=0.;
00323   Double_t dvdz=0.;
00324   if (iuvtxfit>1) {
00325     LinearFit::Weighted(iuvtxfit,uvtxzfit,uvtxfit,uvtxphfit,uvtxparm);
00326     dudz = uvtxparm[1];
00327   }
00328   if (ivvtxfit>1) {
00329     LinearFit::Weighted(ivvtxfit,vvtxzfit,vvtxfit,vvtxphfit,vvtxparm);
00330     dvdz = vvtxparm[1];
00331   }
00332   Double_t uendparm[2];
00333   Double_t vendparm[2];
00334   if (iuendfit>1) {
00335     LinearFit::Weighted(iuendfit,uendzfit,uendfit,uendphfit,uendparm);
00336  
00337   }
00338   if (ivendfit>1) {
00339     LinearFit::Weighted(ivendfit,vendzfit,vendfit,vendphfit,vendparm);
00340  
00341   }
00342 
00343   // calculate direction cosines from linear fits above
00344   Double_t cosz=1./sqrt(1.+dudz*dudz+dvdz*dvdz);
00345   Double_t cosu=dudz*cosz;
00346   Double_t cosv=dvdz*cosz;
00347 
00348   csh.SetDirCosU(cosu);
00349   csh.SetDirCosV(cosv);
00350   csh.SetDirCosZ(cosz);
00351 
00352   csh.SetEndDirCosU(cosu);
00353   csh.SetEndDirCosV(cosv);
00354   csh.SetEndDirCosZ(cosz);
00355 
00356  // set the vertex planes and z position.  The vertex Z position
00357   // is the position of the most upstream strip.  This is not
00358   // correct for cosmics!
00359 
00360   PlexStripEndId stripid = begstrip->GetStripEndId();
00361   UgliScintPlnHandle planehandle = ugh.GetScintPlnHandle(stripid);
00362   csh.SetVtxPlane(begstrip->GetPlane());
00363   csh.SetVtxZ(planehandle.GetZ0());
00364 
00365   stripid = endstrip->GetStripEndId();
00366   planehandle = ugh.GetScintPlnHandle(stripid);
00367   csh.SetEndPlane(endstrip->GetPlane());
00368   csh.SetEndZ(planehandle.GetZ0());
00369 
00370   // the vertex time is defined to be the earliest strip begin time in
00371   // the vertex region defined by VtxPlaneSpan
00372   csh.SetVtxT(tvtx);
00373   csh.SetEndT(tend);  
00374   SetUV(&csh);
00375  
00376   // if this is a cosmic event, base direction on timing
00377   if(isCosmic){
00378     Double_t fitparm[2] = {0,0};
00379     Double_t timefitchi2;
00380     FindTimingDirection(csh,
00381                         ac,
00382                         fitparm,
00383                         timefitchi2);
00384     if(fitparm[1]<0){
00385       
00386       MSG("AlgShowerSR",Msg::kDebug) << " swapping vertex " << endl;
00387 
00388       Int_t newvtx = csh.GetEndPlane();
00389       Int_t newend = csh.GetVtxPlane();
00390       Double_t newvtxz = csh.GetEndZ();
00391       Double_t newendz = csh.GetVtxZ();
00392       Double_t newvtxt = csh.GetEndT();
00393       Double_t newendt = csh.GetVtxT();
00394 
00395       csh.SetVtxZ(newvtxz);
00396       csh.SetEndZ(newendz);
00397       csh.SetVtxPlane(newvtx);
00398       csh.SetEndPlane(newend);
00399       csh.SetVtxT(newvtxt);
00400       csh.SetEndT(newendt);  
00401     }
00402   }
00403 
00404   
00405   // the vertex and end U/V positions are based on fits if the fit was done
00406 
00407   csh.SetVtxU(csh.GetU(csh.GetVtxPlane()));
00408   csh.SetVtxV(csh.GetV(csh.GetVtxPlane()));
00409   
00410   csh.SetEndU(csh.GetU(csh.GetEndPlane()));
00411   csh.SetEndV(csh.GetV(csh.GetEndPlane()));
00412 
00413   Calibrate(&csh);
00414   
00415   CandTrackHandle * associatedtrk=0;
00416   csh.CalibrateEnergy(associatedtrk,ac);
00417 
00418   delete [] uvtxfit;
00419   delete [] uvtxphfit;
00420   delete [] uvtxzfit;
00421   delete [] vvtxfit;
00422   delete [] vvtxphfit;
00423   delete [] vvtxzfit;
00424 
00425   delete [] uendfit;
00426   delete [] uendphfit;
00427   delete [] uendzfit;
00428   delete [] vendfit;
00429   delete [] vendphfit;
00430   delete [] vendzfit;
00431 }

void AlgShowerSR::SetT ( CandShowerHandle candshower  )  [protected]

Definition at line 653 of file AlgShowerSR.cxx.

References UgliStripHandle::ClearFiber(), CandHandle::GetDaughterIterator(), UgliStripHandle::GetHalfLength(), UgliGeomHandle::GetStripHandle(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), Msg::kDebug, CandStripHandle::KeyFromPlane(), StripEnd::kNegative, CalTimeType::kNone, CalDigitType::kPE, StripEnd::kPositive, PlaneView::kU, SimFlag::kUnknown, PlaneView::kUnknown, PlaneView::kV, max, min, MSG, Munits::ns, CandShowerHandle::SetT(), PropagationVelocity::Velocity(), ArrivalTime::Weight(), and UgliStripHandle::WlsPigtail().

Referenced by SetUVT().

00654 {
00655 
00656   // we take a weighted average of hits in the same
00657   // plane.  the proper way to do this would be to consider hits coming
00658   // from the same PMT and use the earliest time.
00659 
00660   MSG("Shower",Msg::kDebug) << "starting SetT" << endl;
00661 
00662   TIter stripItr(candshower->GetDaughterIterator());
00663   CandStripHandle *firststrip = 0;
00664   firststrip = dynamic_cast<CandStripHandle*>(stripItr());
00665   UgliGeomHandle *ugh = 0;
00666   SimFlag::SimFlag_t simFlag = SimFlag::kUnknown;
00667   if (firststrip) {
00668     ugh = new UgliGeomHandle(*firststrip->GetVldContext());
00669     simFlag = firststrip->GetVldContext()->GetSimFlag();
00670   }
00671 
00672   Double_t timewgt[2] = {0.,0.};
00673   Double_t time[2] = {0.,0.};
00674   StripEnd::StripEnd_t stripend[2] = {StripEnd::kNegative,StripEnd::kPositive};
00675 
00676   Double_t apos;
00677   ArrivalTime arrtime;
00678   CalTimeType::CalTimeType_t caltimetype = CalTimeType::kNone;
00679   Double_t propagation_velocity = PropagationVelocity::Velocity(simFlag)*Munits::ns;
00680 
00681   CandStripHandleItr orderedstripItr(candshower->GetDaughterIterator());
00682   CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc();
00683   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00684   orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00685   stripKf = 0;
00686 
00687   Int_t oldplane = -1;
00688   PlaneView::PlaneView_t oldview = PlaneView::kUnknown;
00689 
00690   Bool_t first(1);
00691   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) {
00692     if (strip->GetPlaneView()==PlaneView::kU ||
00693         strip->GetPlaneView()==PlaneView::kV) {
00694       if (strip->GetPlaneView()==PlaneView::kU) {
00695         apos = candshower->GetV(strip->GetPlane());
00696       } else {
00697         apos = candshower->GetU(strip->GetPlane());
00698       }
00699       apos=min(apos,4.0); 
00700       apos=max(apos,-4.0);  // sanity check on apos bounds
00701       UgliStripHandle striphandle = ugh->GetStripHandle(strip->GetStripEndId());
00702       caltimetype = strip->GetCalTimeType();
00703       if (first || strip->GetPlane()==oldplane) {
00704         if (first) {
00705           first = 0;
00706           oldplane = strip->GetPlane();
00707           oldview = strip->GetPlaneView();
00708         }
00709         Double_t thisph[2] = {
00710           strip->GetCharge(CalDigitType::kPE,stripend[0]),
00711           strip->GetCharge(CalDigitType::kPE,stripend[1])};
00712         Double_t thiswgt[2] = {
00713           min(0.4,arrtime.Weight(max(1.,thisph[0]))),
00714           min(0.4,arrtime.Weight(max(1.,thisph[1])))};
00715         for (int i=0; i<2; i++) {
00716           if (thisph[i]>0.) {
00717             Double_t plusminus = 1.;
00718             if ((strip->GetPlaneView()==PlaneView::kV && i==1) ||
00719                 (strip->GetPlaneView()==PlaneView::kU && i==0)) {
00720               plusminus = -1.;
00721             }
00722             Double_t thistime = strip->GetBegTime(stripend[i]) - 
00723               1.e-9*(striphandle.ClearFiber(stripend[i]) + 
00724                      striphandle.WlsPigtail(stripend[i]) + 
00725                      (striphandle.GetHalfLength() + plusminus*apos))/propagation_velocity;
00726 
00727             MSG("Shower",Msg::kDebug) 
00728               << "plane " << strip->GetPlane() << " end " << i 
00729               << " raw time " << strip->GetBegTime(stripend[i])*1e9 
00730               << " pulse height " << thisph[i] << " weight " << thiswgt[i] 
00731               << " clearfiber " << striphandle.ClearFiber(stripend[i]) 
00732               << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) 
00733               << " halflength " << striphandle.GetHalfLength() 
00734               << " apos " << apos << " 3D time " << thistime*1e9 << endl;
00735             time[i] += thiswgt[i]*thistime;
00736             timewgt[i] += thiswgt[i];
00737             MSG("Shower",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl;
00738           }
00739         }
00740       }
00741       else {
00742         for (int i=0; i<2; i++) {
00743           if (timewgt[i]>0.) {
00744             time[i] /= timewgt[i];
00745             candshower->SetT(oldplane,stripend[i],time[i]);
00746             MSG("Shower",Msg::kDebug) << "setting time " << oldplane << " " 
00747                                       << stripend[i] << " " << time[i]*1e9 << endl;
00748           }
00749         }
00750         oldplane = strip->GetPlane();
00751         oldview = strip->GetPlaneView();
00752         Double_t thisph[2] = {
00753           strip->GetCharge(CalDigitType::kPE,stripend[0]),
00754           strip->GetCharge(CalDigitType::kPE,stripend[1])};
00755         Double_t thiswgt[2] = {
00756           min(0.4,arrtime.Weight(max(1.,thisph[0]))),
00757           min(0.4,arrtime.Weight(max(1.,thisph[1])))};
00758         for (int i=0; i<2; i++) {
00759           time[i] = 0.;
00760           timewgt[i] = 0.;
00761           if (thisph[i]>0.) {
00762             Double_t plusminus = 1.;
00763             if ((strip->GetPlaneView()==PlaneView::kV && i==1) ||
00764                 (strip->GetPlaneView()==PlaneView::kU && i==0)) {
00765               plusminus = -1.;
00766             }
00767             Double_t thistime = strip->GetBegTime(stripend[i]) - 
00768               1.e-9*(striphandle.ClearFiber(stripend[i]) + 
00769                      striphandle.WlsPigtail(stripend[i]) + 
00770                      (striphandle.GetHalfLength() + plusminus*apos))/propagation_velocity;
00771             MSG("Shower",Msg::kDebug)
00772               << "plane " << strip->GetPlane() << " end " << i 
00773               << " raw time " << strip->GetBegTime(stripend[i])*1e9 
00774               << " pulse height " << thisph[i] << " weight " << thiswgt[i] 
00775               << " clearfiber " << striphandle.ClearFiber(stripend[i]) 
00776               << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) 
00777               << " halflength " << striphandle.GetHalfLength() 
00778               << " apos " << apos << " 3D time " << thistime*1e9 << endl;
00779             time[i] = thiswgt[i]*thistime;
00780             timewgt[i] = thiswgt[i];
00781             MSG("Shower",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl;
00782           }
00783         }
00784       }
00785     }
00786   }
00787   for (int i=0; i<2; i++) {
00788     if (timewgt[i]>0.) {
00789       time[i] /= timewgt[i];
00790       candshower->SetT(oldplane,stripend[i],time[i]);
00791       MSG("Shower",Msg::kDebug) << "setting time " << oldplane << " " 
00792                                 << stripend[i] << " " << time[i]*1e9 << endl;
00793     }
00794   }
00795   if (ugh) delete ugh;
00796 }

void AlgShowerSR::SetUV ( CandShowerHandle candshower  )  [protected]

Definition at line 443 of file AlgShowerSR.cxx.

References det, PlexPlaneId::GetAdjoinScint(), CandRecoHandle::GetBegPlane(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandRecoHandle::GetEndPlane(), PlexPlaneId::GetPlane(), VHS::GetPlane(), UgliGeomHandle::GetScintPlnHandle(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), CandShowerHandle::GetZ(), UgliPlnHandle::GetZ0(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), CandStripHandle::KeyFromPlane(), StripEnd::kNegative, StripEnd::kPositive, PlaneView::kU, PlaneView::kUnknown, PlaneView::kV, Msg::kWarning, max, min, MSG, CandShowerHandle::SetU(), and CandShowerHandle::SetV().

Referenced by RunAlg(), and SetUVT().

00444 {
00445   TIter stripItr(candshower->GetDaughterIterator());
00446   CandStripHandle *firststrip = 0;
00447   firststrip = dynamic_cast<CandStripHandle*>(stripItr());
00448   if (!firststrip) return;  // if no strips, do nothing
00449   
00450   const VldContext *vldc = firststrip->GetVldContext();
00451   Detector::Detector_t det = vldc->GetDetector();
00452 
00453   UgliGeomHandle ugh = UgliGeomHandle(*vldc);
00454 
00455   CandStripHandleItr orderedstripItr(candshower->GetDaughterIterator());
00456   CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc();
00457   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00458   orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00459   stripKf = 0;
00460   Bool_t first(1);
00461   Int_t oldplane = -1;
00462   PlaneView::PlaneView_t oldview = PlaneView::kUnknown;
00463   Float_t tpos = 0.;
00464   Float_t ph = 0.;
00465   Int_t idir=1;
00466   Int_t vtxplane = candshower->GetBegPlane();
00467   Int_t endplane = candshower->GetEndPlane();
00468   if (vtxplane>endplane) idir=-1;
00469   Int_t begplane = min(vtxplane,endplane);
00470   endplane = max(vtxplane,endplane);
00471 
00472   begplane = -1;
00473   endplane = -1;
00474 
00475 // find begin and end planes
00476   while (CandStripHandle *strip =
00477                     dynamic_cast<CandStripHandle*>(orderedstripItr())) {
00478     if (strip->GetPlaneView()==PlaneView::kU ||
00479         strip->GetPlaneView()==PlaneView::kV) {
00480       if (begplane<0 || strip->GetPlane()<begplane) {
00481         begplane = strip->GetPlane();
00482       }
00483       if (endplane<0 || strip->GetPlane()>endplane) {
00484         endplane = strip->GetPlane();
00485       }
00486     }
00487   }
00488   Int_t *planelist = new Int_t[endplane+1];
00489   Int_t *stripendlist = new Int_t[endplane+1];
00490   if(begplane>-1 && endplane>-1){
00491 
00492     for (int i=0; i<=endplane; i++) {
00493       planelist[i] = 0;
00494       stripendlist[i] = 0;
00495     }
00496     
00497     // transverse positions
00498     orderedstripItr.Reset();
00499     
00500     while (CandStripHandle *strip =
00501            dynamic_cast<CandStripHandle*>(orderedstripItr())) {
00502       if (strip->GetPlaneView()==PlaneView::kU ||
00503           strip->GetPlaneView()==PlaneView::kV) {
00504         planelist[strip->GetPlane()] = strip->GetPlaneView();
00505         if (strip->GetNDigit(StripEnd::kNegative)) {
00506           stripendlist[strip->GetPlane()] += (1<<StripEnd::kNegative);
00507         }
00508         if (strip->GetNDigit(StripEnd::kPositive)) {
00509           stripendlist[strip->GetPlane()] += (1<<StripEnd::kPositive);
00510         }
00511         if (first || strip->GetPlane()==oldplane) {
00512           if (first) {
00513             first = 0;
00514             oldplane = strip->GetPlane();
00515             oldview = strip->GetPlaneView();
00516           }
00517           Float_t stripph = strip->GetCharge();
00518           tpos += strip->GetTPos()*stripph;
00519           ph += stripph;
00520         }
00521         else {
00522           tpos /= ph;
00523           if (oldview==PlaneView::kU) {
00524             candshower->SetU(oldplane,tpos);
00525           } else if (oldview==PlaneView::kV) {
00526             candshower->SetV(oldplane,tpos);
00527           }
00528           oldplane = strip->GetPlane();
00529           oldview = strip->GetPlaneView();
00530           Float_t stripph = strip->GetCharge();
00531           tpos = strip->GetTPos()*stripph;
00532           ph = stripph;
00533         }
00534       }
00535     }
00536     if (ph>0.) {
00537       tpos /= ph;
00538       if (oldview==PlaneView::kU) {
00539         candshower->SetU(oldplane,tpos);
00540       } else if (oldview==PlaneView::kV) {
00541         candshower->SetV(oldplane,tpos);
00542       }
00543     }
00544     
00545     // create a starting and stopping PlexPlaneId
00546     Int_t idirAdjoin = (begplane>endplane) ? -1 : 1; // direction
00547     PlexPlaneId planeid(det,begplane,false);
00548     PlexPlaneId planeidend(det,endplane,false);
00549     PlexPlaneId planeid_toofar = planeidend.GetAdjoinScint(idirAdjoin);
00550     
00551     // after each step, if not reached end, move to adjoining *scint* plane
00552     for ( ; planeid != planeid_toofar ; 
00553           planeid  = planeid.GetAdjoinScint(idirAdjoin) ) {
00554       
00555       // quit if hit something meaningless
00556       if ( ! planeid.IsValid() ) break;
00557       
00558       UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(planeid);
00559       
00560       if ( ! scintpln.IsValid() ) continue;
00561       
00562       Int_t iplane = planeid.GetPlane();
00563       
00564       Float_t zpos = scintpln.GetZ0();
00565       for (Int_t planeview=PlaneView::kU; planeview<=PlaneView::kV;
00566            planeview++) {
00567         if (planelist[iplane]!=planeview) {
00568           Int_t found = 0;
00569           Int_t nfound = 0;
00570           Int_t jplane[2];
00571           
00572           // search upstream
00573           for (Int_t iplane1=iplane-1; iplane1>=begplane && !found;
00574                iplane1--) {
00575             if (planelist[iplane1]==planeview) {
00576               found++;
00577               jplane[nfound] = iplane1;
00578               nfound++;
00579             }
00580           }
00581           Int_t ndown = 1;
00582           if (!found) {
00583             ndown++;
00584           }
00585           found = 0;
00586           
00587           // search downstream
00588           for (Int_t iplane1=iplane+1; iplane1<=endplane && found<ndown;
00589                iplane1++) {
00590             if (planelist[iplane1]==planeview) {
00591               found++;
00592               jplane[nfound] = iplane1;
00593               nfound++;
00594             }
00595           }
00596           if (!found) {
00597             for (Int_t iplane1=jplane[0]-1; iplane1>=begplane && !found;
00598                  iplane1--) {
00599               if (planelist[iplane1]==planeview) {
00600                 found++;
00601                 jplane[nfound] = iplane1;
00602                 nfound++;
00603               }
00604             }
00605           }
00606           if (nfound==1) {
00607             if (planeview==PlaneView::kU) {
00608               candshower->SetU(iplane,candshower->GetU(jplane[0]));
00609             } else {
00610               candshower->SetV(iplane,candshower->GetV(jplane[0]));
00611             }
00612           }
00613           else {
00614             Float_t z[2];
00615             z[0] = candshower->GetZ(jplane[0]);
00616             z[1] = candshower->GetZ(jplane[1]);
00617             Float_t x[2];
00618             if (planeview==PlaneView::kU) {
00619               x[0] = candshower->GetU(jplane[0]);
00620               x[1] = candshower->GetU(jplane[1]);
00621               Float_t denom = (z[1]-z[0]);
00622               if (denom) candshower->SetU(iplane,
00623                                           x[0] + (zpos-z[0]) * ((x[1]-x[0])/denom));
00624               else {
00625                 MSG("Shower",Msg::kWarning) << endl
00626                                             << "Denominator (z[1]-z[0]) is zero."
00627                                             << "  SetU not called." << endl;
00628               }
00629             } else {
00630               x[0] = candshower->GetV(jplane[0]);
00631               x[1] = candshower->GetV(jplane[1]);
00632               Float_t denom = (z[1]-z[0]);
00633               if (denom) candshower->SetV(iplane,
00634                                           x[0] + (zpos-z[0]) * ((x[1]-x[0])/denom));
00635               else {
00636                 MSG("Shower",Msg::kWarning) << endl
00637                                             << "Denominator is (z[1]-z[0]) zero."
00638                                             << "  SetV not called." << endl;
00639               }
00640             }
00641           }
00642         }
00643       }
00644     }
00645   }
00646   delete [] planelist;
00647   delete [] stripendlist;
00648 
00649 }

void AlgShowerSR::SetUVT ( CandShowerHandle candshower  )  [protected]

Definition at line 435 of file AlgShowerSR.cxx.

References SetT(), and SetUV().

Referenced by AlgShowerSS::RunAlg().

00436 {
00437   this->SetUV(candshower);
00438   this->SetT(candshower);
00439 }

void AlgShowerSR::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Reimplemented in AlgShowerSS.

Definition at line 1021 of file AlgShowerSR.cxx.

01022 {
01023 }


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1