AtmosCalculator.cxx

Go to the documentation of this file.

#include "AtmosCalculator.h"
#include "ShowerMOI.h"
#include "ShowerTrace.h"
#include "AtNuEvent/AtmosStrip.h"

#include "MessageService/MsgService.h"
#include "TClonesArray.h"
#include "TObject.h"

#include <cmath>
#include <cstdlib>

ClassImp(AtmosCalculator)

CVSID("$Id: AtmosCalculator.cxx,v 1.18 2009/02/28 21:46:10 gmieg Exp $");

AtmosCalculator::AtmosCalculator()
{
  MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::AtmosCalculator()  " << endl;
}

AtmosCalculator::~AtmosCalculator()
{
  MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::~AtmosCalculator()  " << endl;
}

void AtmosCalculator::EventProperties(AtmosEvent* MyEvent, TClonesArray* StripList)
{
  MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::EventProperties(...)  " << endl;
  
  // pulse height stuff

  double dQ(0.0),dQe(0.0),dQw(0.0);
  double totalQ(0.0),eastQ(0.0),westQ(0.0);
  double singleQ(0.0),singleQU(0.0),singleQV(0.0);
  double doubleQ(0.0),doubleQU(0.0),doubleQV(0.0);

  int ListEnd = StripList->GetLast()+1;
  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      if(strip.Plane>-1 && strip.Plane<486)
        {
          dQe=strip.Sigcorr[0]/65.0;
          dQw=strip.Sigcorr[1]/65.0;
          dQ=(strip.Sigcorr[0]+strip.Sigcorr[1])/65.0; 
            // hack a rough conversion from SigCorr -> PEcorr

          totalQ += dQ;
          eastQ += dQe; westQ += dQw;

          if( strip.Ndigits==1 )
            {
              singleQ += dQ;
              if( strip.View==0 ) singleQU += dQ;
              if( strip.View==1 ) singleQV += dQ;
            }
          if( strip.Ndigits==2 )
            {
              doubleQ += dQ;
              if( strip.View==0 ) doubleQU += dQ;
              if( strip.View==1 ) doubleQV += dQ;
            }

        }
    }

  MyEvent->RecoInfo.TotalCharge = totalQ;
  MyEvent->RecoInfo.EastSideCharge = eastQ;
  MyEvent->RecoInfo.WestSideCharge = westQ;
  MyEvent->RecoInfo.SingleEndedCharge = singleQ;
  MyEvent->RecoInfo.SingleEndedChargeUview = singleQU;
  MyEvent->RecoInfo.SingleEndedChargeVview = singleQV;
  MyEvent->RecoInfo.DoubleEndedCharge = doubleQ;
  MyEvent->RecoInfo.DoubleEndedChargeUview = doubleQU;
  MyEvent->RecoInfo.DoubleEndedChargeVview = doubleQV;

  // multiple muon reconstruction goes here
  MyEvent->RecoInfo.MultiMuReco = 0;
  MyEvent->RecoInfo.MultiMuTracks = 0;
  MyEvent->RecoInfo.MultiMuDirCosU = 0.0;
  MyEvent->RecoInfo.MultiMuDirCosV = 0.0;
  MyEvent->RecoInfo.MultiMuDirCosX = 0.0;
  MyEvent->RecoInfo.MultiMuDirCosY = 0.0;
  MyEvent->RecoInfo.MultiMuDirCosZ = 0.0;

  return;
}

void AtmosCalculator::TrackProperties(AtmosTrack* MyTrack, TClonesArray* StripList)
{
  MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::TrackProperties(...)  " << endl;

  // Get Track Index
  int index=MyTrack->Index;

  // minimum/maximum plane number
  int vtxpln=MyTrack->VtxPlane; 
  int endpln=MyTrack->EndPlane;
  int minpln=(MyTrack->VtxPlane<MyTrack->EndPlane)?MyTrack->VtxPlane:MyTrack->EndPlane;
  int maxpln=(MyTrack->VtxPlane<MyTrack->EndPlane)?MyTrack->EndPlane:MyTrack->VtxPlane;
  double dir=(MyTrack->VtxPlane>MyTrack->EndPlane)?-1.0:1.0;
  MyTrack->MinPlaneNumber = minpln;
  MyTrack->MaxPlaneNumber = maxpln;

  // Set vertex/end position/direction
  double vtxdiru=MyTrack->VtxDirCosU; double enddiru=MyTrack->EndDirCosU;
  double vtxdirv=MyTrack->VtxDirCosV; double enddirv=MyTrack->EndDirCosV;
  double vtxdirz=MyTrack->VtxDirCosZ; double enddirz=MyTrack->EndDirCosZ;
  double vtxu=MyTrack->VtxU; double endu=MyTrack->EndU;
  double vtxv=MyTrack->VtxV; double endv=MyTrack->EndV;
  double vtxx=MyTrack->VtxX; double endx=MyTrack->EndX;
  double vtxy=MyTrack->VtxY; double endy=MyTrack->EndY;
  double vtxz=MyTrack->VtxZ; double endz=MyTrack->EndZ;
  
  // Closest distance to edge of detector
  MSG("AtmosCalculator",Msg::kVerbose) << "    Distance from Edge of Detector" << endl;
  double rmin,temp;

  rmin=4.0; 
  temp=4.0-vtxu; if(temp<rmin) rmin=temp;
  temp=4.0+vtxu; if(temp<rmin) rmin=temp;
  temp=4.0-vtxv; if(temp<rmin) rmin=temp;
  temp=4.0+vtxv; if(temp<rmin) rmin=temp;
  temp=4.0-vtxx; if(temp<rmin) rmin=temp;
  temp=4.0+vtxx; if(temp<rmin) rmin=temp;
  temp=4.0-vtxy; if(temp<rmin) rmin=temp;
  temp=4.0+vtxy; if(temp<rmin) rmin=temp;
  MyTrack->VtxDistToEdge = rmin; // VtxDistToEdge

  rmin=4.0; 
  temp=4.0-endu; if(temp<rmin) rmin=temp;
  temp=4.0+endu; if(temp<rmin) rmin=temp;
  temp=4.0-endv; if(temp<rmin) rmin=temp;
  temp=4.0+endv; if(temp<rmin) rmin=temp;
  temp=4.0-endx; if(temp<rmin) rmin=temp;
  temp=4.0+endx; if(temp<rmin) rmin=temp;
  temp=4.0-endy; if(temp<rmin) rmin=temp;
  temp=4.0+endy; if(temp<rmin) rmin=temp; 
  MyTrack->EndDistToEdge = rmin; // EndDistToEdge
 
  // Calculate Traces
  // ================

  double invSqrt2 = pow(1./2.,0.5);
  double m[2]; double c[2]; double Sign;
  double Coord[3]; double Trace[3];
  bool TraceSuccess;

  // For vertex
  if(vtxdirz!=0.) {
    m[0]=vtxdiru/vtxdirz;
    m[1]=vtxdirv/vtxdirz;
    c[0]=vtxu-(m[0]*vtxz);
    c[1]=vtxv-(m[1]*vtxz);
    Coord[0] = vtxu; Coord[1] = vtxv; Coord[2] = vtxz;
    Trace[0] = 1.e99; Trace[1] = 1.e99; Trace[2] = 1.e99;

    TraceSuccess = CalculateTrace(m,c,Coord,Trace);
    if(TraceSuccess==true) {
      Sign=1.;
      if(fabs(Coord[0])>4. || fabs(Coord[1])>4. || fabs(invSqrt2*(Coord[0]+Coord[1]))>4. 
         || fabs(invSqrt2*(Coord[0]-Coord[1]))>4.) {Sign=-1.;}
        
      MyTrack->VtxTrace=(Sign*pow(pow(Trace[0],2) + pow(Trace[1],2) + pow(Trace[2],2), 0.5));
      MyTrack->VtxTraceZ=(Sign*fabs(Trace[2]));
    }
  }
    
  // For end
  if(enddirz!=0.) {
    m[0] = enddiru/enddirz;
    m[1] = enddirv/enddirz;
    c[0] = endu-(m[0]*endz);
    c[1] = endv-(m[1]*endz);
    Coord[0] = endu; Coord[1] = endv; Coord[2] = endz;
    Trace[0] = 1.e99; Trace[1]=1.e99; Trace[2]=1.e99;
    
    TraceSuccess = CalculateTrace(m,c,Coord,Trace);
    if(TraceSuccess==true) {
      Sign=1.;
      if(fabs(Coord[0])>4. || fabs(Coord[1])>4. || fabs(invSqrt2*(Coord[0]+Coord[1]))>4. 
         || fabs(invSqrt2*(Coord[0]-Coord[1]))>4.) {Sign=-1.;}

      MyTrack->EndTrace=(Sign*pow(pow(Trace[0],2) + pow(Trace[1],2) + pow(Trace[2],2), 0.5));
      MyTrack->EndTraceZ=(Sign*fabs(Trace[2]));
    }
  }



  // Track containment variables
  // ===========================

  MSG("AtmosCalculator",Msg::kVerbose) << "    Track Containment Variables" << endl;
  Int_t nSingleEnds(0),nDoubleEnds(0),nDigits(0);
  Double_t dT,dU,dV,dQ,vtxdUmax(0.),vtxdVmax(0.),enddUmax(0.),enddVmax(0.);
  Double_t vtxUwt2(0.),vtxUwt(0.),vtxUw(0.);
  Double_t vtxVwt2(0.),vtxVwt(0.),vtxVw(0.);
  Double_t endUwt2(0.),endUwt(0.),endUw(0.);
  Double_t endVwt2(0.),endVwt(0.),endVw(0.);
  Double_t totTrkPH(0.),totTrkAssocPH(0.),totPlnPH(0.);
  Double_t trkAssocPH[500],trkPH[500],plnPH[500];
  Int_t bstrp[500],estrp[500];
  Int_t nstrps[500],ntrkstrps[500];
  Int_t planeview[500];
  memset(trkAssocPH, 0, 500*sizeof(Double_t));
  memset(trkPH, 0, 500*sizeof(Double_t));
  memset(plnPH, 0, 500*sizeof(Double_t));
  memset(bstrp, -1, 500*sizeof(Int_t));
  memset(estrp, -1, 500*sizeof(Int_t));
  memset(nstrps, 0, 500*sizeof(Int_t));
  memset(ntrkstrps, 0, 500*sizeof(Int_t));
  memset(planeview, -1, 500*sizeof(Int_t));
  int thispln;

  int ListEnd = StripList->GetLast()+1;
  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln>-1 && thispln<486)
        {

          ++nstrps[thispln];
 
          if(strip.View==0)
            {
              if(planeview[thispln]<0) planeview[thispln]=0;
            }

          if(strip.View==1)
            {
              if(planeview[thispln]<0) planeview[thispln]=1;
            }

          if((strip.Trk&index)==index) 
            {
              if(bstrp[thispln]<0||strip.Strip<bstrp[thispln]) bstrp[thispln]=strip.Strip;
              if(estrp[thispln]<0||strip.Strip>estrp[thispln]) estrp[thispln]=strip.Strip;
              if(strip.Ndigits==1) ++nSingleEnds; if(strip.Ndigits==2) ++nDoubleEnds;
              nDigits += strip.Ndigits;
              ++ntrkstrps[thispln];
            }
        }
    }

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      dQ=(strip.Sigcorr[0]+strip.Sigcorr[1])/65.0; // hack a rough conversion 
                                                   // from SigCorr -> PEcorr

      // Record pulse height information
      if(thispln>-1 && thispln<486)
        {
          // strips associated with track (+/-1 strip window...)
          if( bstrp[thispln]>-1 && estrp[thispln]>-1 && strip.Strip>bstrp[thispln]-2 && strip.Strip<estrp[thispln]+2)
            {
              trkAssocPH[thispln]+=dQ;
              totTrkAssocPH+=dQ;
            }

          // strips on track (+/-0 strip window...)
          if( bstrp[thispln]>-1 && estrp[thispln]>-1 && strip.Strip>bstrp[thispln]-1 && strip.Strip<estrp[thispln]+1)
            {
              trkPH[thispln]+=dQ;
              totTrkPH+=dQ;
            }

          // total pulse height
          plnPH[thispln]+=dQ;
          totPlnPH+=dQ;
        }
      
      // Handle the Track containment information
      // (could require the strips to NOT be on the track here...)
      if( abs(vtxpln-thispln)<5 && !(vtxpln<249 && thispln>249) && !(vtxpln>249 && thispln<249) )
        {//strip must be within four of the vtx and not on the other side of the SM gap - could do this more neatly with Z
          if(strip.View==0)
            { 
              dU=strip.T-vtxu; vtxdUmax=(dU>vtxdUmax)?dU:vtxdUmax;
              dT=(strip.T)-(vtxu+dir*vtxdiru*(strip.Z-vtxz));
              vtxUwt2+=dQ*dT*dT; vtxUwt+=dQ*dT; vtxUw+=dQ;
            }
          if(strip.View==1)
            { 
              dV=strip.T-vtxv; vtxdVmax=(dV>vtxdVmax)?dV:vtxdVmax;
              dT=(strip.T)-(vtxv+dir*vtxdirv*(strip.Z-vtxz));
              vtxVwt2+=dQ*dT*dT; vtxVwt+=dQ*dT; vtxVw+=dQ;
            }
        }

      if( abs(endpln-thispln)<5 && !(endpln<249 && thispln>249) && !(endpln>249 && thispln<249) )
        {//strip must be within four of the end and not on the other side of the SM gap - could do this more neatly with Z
          if(strip.View==0)
            { 
              dU=strip.T-endu; enddUmax=(dU>enddUmax)?dU:enddUmax;
              dT=(strip.T)-(endu+dir*enddiru*(strip.Z-endz));
              endUwt2+=dQ*dT*dT; endUwt+=dQ*dT; endUw+=dQ;
            }
          if(strip.View==1)
            { 
              dV=strip.T-endv; enddVmax=(dV>enddVmax)?dV:enddVmax;
              dT=(strip.T)-(endv+dir*enddirv*(strip.Z-endz));
              endVwt2+=dQ*dT*dT; endVwt+=dQ*dT; endVw+=dQ;
            }
        }
    }

  // VtxRmax, EndRmax
  MyTrack->VtxRmax = pow( vtxdUmax*vtxdUmax + vtxdVmax*vtxdVmax, 0.5 );
  MyTrack->EndRmax = pow( enddUmax*enddUmax + enddVmax*enddVmax, 0.5 );

  // VtxUmean, VtxUwidth
  if(vtxUw>0.0)
    {
      MyTrack->VtxUmean=vtxUwt/vtxUw;
      if((vtxUwt2/vtxUw)>(vtxUwt/vtxUw)*(vtxUwt/vtxUw))
        {
          MyTrack->VtxUwidth=sqrt((vtxUwt2/vtxUw)-(vtxUwt/vtxUw)*(vtxUwt/vtxUw));
        }
    }

  // VtxVmean, VtxVwidth
  if(vtxVw>0.0)
    {
      MyTrack->VtxVmean=vtxVwt/vtxVw;
      if((vtxVwt2/vtxVw)>(vtxVwt/vtxVw)*(vtxVwt/vtxVw))
        {
          MyTrack->VtxVwidth = sqrt((vtxVwt2/vtxVw)-(vtxVwt/vtxVw)*(vtxVwt/vtxVw));
        }
    }


  // EndUmean, EndUwidth
  if(endUw>0.0)
    {
      MyTrack->EndUmean=endUwt/endUw;
      if((endUwt2/endUw)>(endUwt/endUw)*(endUwt/endUw))
        {
          MyTrack->EndUwidth=sqrt((endUwt2/endUw)-(endUwt/endUw)*(endUwt/endUw));
        }
    }

  // EndVmean, EndVwidth
  if(endVw>0.0)
    {
      MyTrack->EndVmean=endVwt/endVw;  
      if((endVwt2/endVw)>(endVwt/endVw)*(endVwt/endVw))
        {
          MyTrack->EndVwidth=sqrt((endVwt2/endVw)-(endVwt/endVw)*(endVwt/endVw));
        }
    }

  // Pulse height variables
  MSG("AtmosCalculator",Msg::kVerbose) << "    Pulse Height Variables" << endl;
  double vtxqmax=0.0,endqmax=0.0;
  for(int ipln=0; ipln<500; ++ipln)
    {
      if( abs(vtxpln-ipln)<5 && !(vtxpln<249 && ipln>249) && !(vtxpln>249 && ipln<249) )
        {
          if( plnPH[ipln]>vtxqmax ) vtxqmax=plnPH[ipln];
        }
      if( abs(endpln-ipln)<5 && !(endpln<249 && ipln>249) && !(endpln>249 && ipln<249) )
        {
          if( plnPH[ipln]>endqmax ) endqmax=plnPH[ipln];
        }
    }

  MyTrack->VtxQmax = vtxqmax;          // VtxQmax
  MyTrack->EndQmax = endqmax;          // EndQmax
  MyTrack->TrkPH = totTrkPH;           // TrkPH
  MyTrack->AssocTrkPH = totTrkAssocPH; // AssocTrkPH
  if(totPlnPH>0.0) MyTrack->AssocTrkPHfrac = totTrkAssocPH/totPlnPH;

  // Additional plane/strip variables;
  MSG("AtmosCalculator",Msg::kVerbose) << "    Additional Plane/Strip Variables" << endl;
  int uPlanes(0),vPlanes(0);
  int nPlanesTrkOnly(0),nPlanesTrkGaps(0);
  for(int ipln=0; ipln<500; ++ipln)
    {
      if(ipln>=minpln && ipln<=maxpln)
        { 
          if(ntrkstrps[ipln]>0)
            {
              if(planeview[ipln]==0) ++uPlanes;
              if(planeview[ipln]==1) ++vPlanes;
              if(ntrkstrps[ipln]==nstrps[ipln]) ++nPlanesTrkOnly;
            }
          else
            {
              ++nPlanesTrkGaps;
            }
        }
    }

  MyTrack->Ndigits = nDigits;                 // Ndigits
  MyTrack->NstripsSingleEnded = nSingleEnds;  // NstripsSingleEnded
  MyTrack->NstripsDoubleEnded = nDoubleEnds;  // NstripsDoubleEnded
  MyTrack->NplanesTrackOnly = nPlanesTrkOnly; // NplanesTrackOnly
  MyTrack->NplanesTrackGaps = nPlanesTrkGaps; // NplanesTrackGaps
  MyTrack->NplanesUview = uPlanes;            // NplanesUview
  MyTrack->NplanesVview = vPlanes;            // NplanesVview
  if(uPlanes+vPlanes>0) MyTrack->UVassymetry = abs(uPlanes-vPlanes)/(0.5*(uPlanes+vPlanes));                                 // UVassymetry
  
  // Number of track-like planes
  Int_t ntrkplns(0);
  for(int ipln=0; ipln<500; ++ipln)
    {
      if( plnPH[ipln]>0.0 && plnPH[ipln]<80.0 ) 
        {
          if(trkAssocPH[ipln]/plnPH[ipln]>0.8) ++ntrkplns;
        }
    }

  MyTrack->TrkLikePlanes = ntrkplns; // TrkLikePlanes

  
  // Linear Fits
  // ===========

  //double dstrp = 0.041/sqrt(12.0);
  double dstrpSQ = 0.041*0.041/12.0;

  double u,v,x,y,z,w,du,dv;
  double mu,mv,cu,cv,denom;
  double Uchi2,Vchi2;
  int Undf,Vndf;
  int Upts,Vpts;

  // Linear Fit 
  // (Separate Linear Fit In Each View)
  MSG("AtmosCalculator",Msg::kVerbose) << "    Linear Fit " << endl;
  double Uw(0.0), Uwz(0.0), Uwzz(0.0), Uwu(0.0), Uwuu(0.0), Uwzu(0.0);
  double Vw(0.0), Vwz(0.0), Vwzz(0.0), Vwv(0.0), Vwvv(0.0), Vwzv(0.0);
  
  Upts=0; Vpts=0; 

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln>-1 && thispln<486)
        {
          if((strip.Trk&index)==index) 
            {
              z=strip.Z; 
              w=1.0;

              if(strip.View==0)
                {
                  u=strip.T;
                  Uw+=w;
                  Uwu+=w*u; Uwz+=w*z; Uwuu+=w*u*u;
                  Uwzu+=w*u*z; Uwzz+=w*z*z;
                  ++Upts;
                }
              
              if(strip.View==1)
                {
                  v=strip.T;
                  Vw+=w;
                  Vwv+=w*v; Vwz+=w*z; Vwvv+=w*v*v;
                  Vwzv+=w*v*z; Vwzz+=w*z*z;
                  ++Vpts;
                }
              
            }
        }
    }

  mu=0.0; mv=0.0; cu=0; cv=0;

  if(Upts>1 && Vpts>1)
    {
      denom = Uw*Uwzz-Uwz*Uwz;
      if(fabs(denom)>0.0) { mu=(Uw*Uwzu-Uwz*Uwu)/denom; cu=(Uwu*Uwzz-Uwz*Uwzu)/denom; }
      denom = Vw*Vwzz-Vwz*Vwz;
      if(fabs(denom)>0.0) { mv=(Vw*Vwzv-Vwz*Vwv)/denom; cv=(Vwv*Vwzz-Vwz*Vwzv)/denom; }
      // (calculate error properly now)
      // err=(Uwuu+mu*mu*Uwzz+Uw*cu*cu+2.0*(mu*cu*Uwz-mu*Uwzu-cu*Uwu))+(Vwvv+mv*mv*Vwzz+Vw*cv*cv+2.0*(mv*cv*Vwz-mv*Vwzv-cv*Vwv));
      // chi2=err;
      // ndf=Upts+Vpts-4;
    }
  
  double precou=mu; double precov=mv; double precoz=1.0;
  double precos=pow(precou*precou+precov*precov+1.0, -0.5);
  precou*=precos; precov*=precos; precoz*=precos;
  MyTrack->LinearDirCosU=precou*dir; // LinearDirCosU
  MyTrack->LinearDirCosV=precov*dir; // LinearDirCosV
  MyTrack->LinearDirCosZ=precoz*dir; // LinearDirCosZ

  // Unconstrained Linear Fit 
  //  (chi2 about linear fit to track)
  MSG("AtmosCalculator",Msg::kVerbose) << "    Linear Fit (unconstrained) " << endl;

  mu = 0.0;
  if( MyTrack->LinearDirCosZ!=0 )
    mu = MyTrack->LinearDirCosU/MyTrack->LinearDirCosZ; 
  cu = vtxu-mu*vtxz;

  if( MyTrack->LinearDirCosZ!=0 )
    mv = MyTrack->LinearDirCosV/MyTrack->LinearDirCosZ; 
  cv = vtxv-mv*vtxz;

  Uchi2=0.0; Vchi2=0.0; 
  Undf=0; Vndf=0;

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln>-1 && thispln<486)
        {
          if((strip.Trk&index)==index) 
            {
              if( strip.View==0 )
                {
                  du=strip.T-(mu*strip.Z+cu);
                  Uchi2+=(du*du)/dstrpSQ; //(dstrp*dstrp); 
                  ++Undf;
                }

              if( strip.View==1 )
                {
                  dv=strip.T-(mv*strip.Z+cv);
                  Vchi2+=dv*dv/dstrpSQ; //(dstrp*dstrp); 
                  ++Vndf;
                }

            }
        }
    }
  
  MyTrack->LinearDirFitChisqU = Uchi2;       // LinearDirFitChisqU
  MyTrack->LinearDirFitNdfU   = Undf;        // LinearDirFitNdfU
  MyTrack->LinearDirFitChisqV = Vchi2;       // LinearDirFitChisqV
  MyTrack->LinearDirFitNdfV   = Vndf;        // LinearDirFitNdfV
  MyTrack->LinearDirFitChisq  = Uchi2+Vchi2; // LinearDirFitChisq
  MyTrack->LinearDirFitNdf    = Undf+Vndf;   // LinearDirFitNdf

  // Constrained Linear Fits 
  //  (chi2 about linear fit using reconstructed track direction)
  MSG("AtmosCalculator",Msg::kVerbose) << "    Linear Fit (vertex) " << endl;

  mu = 0.0;
  if( MyTrack->VtxDirCosZ!=0 )
    mu = MyTrack->VtxDirCosU/MyTrack->VtxDirCosZ;
  cu = vtxu-mu*vtxz;

  mv = 0.0;
  if( MyTrack->VtxDirCosZ!=0 ) 
    mv = MyTrack->VtxDirCosV/MyTrack->VtxDirCosZ; 
  cv = vtxv-mv*vtxz;

  Uchi2=0.0; Vchi2=0.0; 
  Undf=0; Vndf=0;

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln>-1 && thispln<486)
        {
          if((strip.Trk&index)==index) 
            {
              if( strip.View==0 )
                {
                  du=strip.T-(mu*strip.Z+cu);
                  Uchi2+=(du*du)/dstrpSQ; //(dstrp*dstrp); 
                  ++Undf;
                }

              if( strip.View==1 )
                {
                  dv=strip.T-(mv*strip.Z+cv);
                  Vchi2+=dv*dv/dstrpSQ; //(dstrp*dstrp); 
                  ++Vndf;
                }

            }
        }
    }

  MyTrack->VtxLinearDirFitChisqU = Uchi2; // VtxLinearDirFitChisqU
  MyTrack->VtxLinearDirFitNdfU   = Undf;  // VtxLinearDirFitNdfU
  MyTrack->VtxLinearDirFitChisqV = Vchi2; // VtxLinearDirFitChisqV
  MyTrack->VtxLinearDirFitNdfV   = Vndf;  // VtxLinearDirFitNdfV

  mu = 0.0;
  if( MyTrack->EndDirCosZ!=0 )
    mu = MyTrack->EndDirCosU/MyTrack->EndDirCosZ;
  cu = endu-mu*endz;

  mv = 0.0;
  if( MyTrack->EndDirCosZ!=0 ) 
    mv = MyTrack->EndDirCosV/MyTrack->EndDirCosZ; 
  cv = endv-mv*endz;

  Uchi2=0.0; Vchi2=0.0; 
  Undf=0; Vndf=0;

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln>-1 && thispln<486)
        {
          if((strip.Trk&index)==index) 
            {
              if( strip.View==0 )
                {
                  du=strip.T-(mu*strip.Z+cu);
                  Uchi2+=(du*du)/dstrpSQ; // (dstrp*dstrp); 
                  ++Undf;
                }

              if( strip.View==1 )
                {
                  dv=strip.T-(mv*strip.Z+cv);
                  Vchi2+=dv*dv/dstrpSQ;   // (dstrp*dstrp); 
                  ++Vndf;
                }

            }
        }
    }

  MyTrack->EndLinearDirFitChisqU = Uchi2; // EndLinearDirFitChisqU
  MyTrack->EndLinearDirFitNdfU   = Undf;  // EndLinearDirFitNdfU
  MyTrack->EndLinearDirFitChisqV = Vchi2; // EndLinearDirFitChisqV
  MyTrack->EndLinearDirFitNdfV   = Vndf;  // EndLinearDirFitNdfV
  
  // Zero Curvature Fit
  //  (join up the track vertex and end positions - 
  //   calculate the RMS of the strips about this line)
  MSG("AtmosCalculator",Msg::kVerbose) << "    Zero Curvature Fit" << endl;

  mu = (endu-vtxu)/(endz-vtxz); 
  cu = vtxu-mu*vtxz;
  mv = (endv-vtxv)/(endz-vtxz);
  cv = vtxv-mv*vtxz;
 
  Uchi2=0.0; Vchi2=0.0; 
  Undf=0; Vndf=0;

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln>-1 && thispln<486)
        {
          if((strip.Trk&index)==index) 
            {
              if( strip.View==0 )
                {
                  du=strip.T-(mu*strip.Z+cu);
                  Uchi2+=(du*du)/dstrpSQ; //(dstrp*dstrp); 
                  ++Undf;
                }

              if( strip.View==1 )
                {
                  dv=strip.T-(mv*strip.Z+cv);
                  Vchi2+=dv*dv/dstrpSQ; //(dstrp*dstrp); 
                  ++Vndf;
                }

            }
        }
    }

  MyTrack->ZeroCurveChi2 = Uchi2 + Vchi2;
  MyTrack->ZeroCurveNdf = Undf + Vndf;


  // Digit Containment Variables
  // ===========================

  MSG("AtmosCalculator",Msg::kVerbose) << "    Digit Containment Variables" << endl;
  double wPlaneChargeSigCorr[500],wPlaneMeanTPos[500],wPlaneZPos[500];
  memset(wPlaneChargeSigCorr, 0, 500*sizeof(double));
  memset(wPlaneMeanTPos, 0, 500*sizeof(double));
  memset(wPlaneZPos, 0, 500*sizeof(double));
  double stripU(0.0),stripV(0.0);
  double nonfidTrkQ(0.0),totTrkQ(0.0);
  bool UCsides(false),UCends(false);
  double Xdis(0.3);
  int Xpln(5);
  double VtxZdigit[3]={vtxx,vtxy,MyTrack->VtxPlane};
  double EndZdigit[3]={endx,endy,MyTrack->EndPlane};
  double VtxRdigit[3]={vtxx,vtxy,MyTrack->VtxPlane};
  double EndRdigit[3]={endx,endy,MyTrack->EndPlane};
  int fPlaneMin = 500,fPlaneMax = 0;
 
  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      wPlaneZPos[strip.Plane] = strip.Z;
      wPlaneChargeSigCorr[strip.Plane] += (strip.Sigcorr[0]+strip.Sigcorr[1]);
      wPlaneMeanTPos[strip.Plane] += (strip.T*(strip.Sigcorr[0]+strip.Sigcorr[1]));
      if(strip.Plane>fPlaneMax) fPlaneMax = strip.Plane;
      if(strip.Plane<fPlaneMin) fPlaneMin = strip.Plane;
    }

  for(int ipln=fPlaneMin; ipln<=fPlaneMax; ++ipln)
    {
      if(wPlaneChargeSigCorr[ipln]>0.0)
        {
          wPlaneMeanTPos[ipln]/=wPlaneChargeSigCorr[ipln];
        }
    }
 
  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));

      const int plane = strip.Plane;
      const double tpos = strip.T;
      
      double nz[2]={-1.0, -1.0};
      double ntpos[2]={-1.0, -1.0};
      double stp_guess_z=wPlaneZPos[plane];
      double stp_guess_tpos=-9999.0;

      if(plane>0 && plane<500 && wPlaneChargeSigCorr[plane-1]>0.0 && wPlaneChargeSigCorr[plane+1]>0.0)
        {
          nz[0] = wPlaneZPos[plane-1];
          nz[1] = wPlaneZPos[plane+1];
          ntpos[0] = wPlaneMeanTPos[plane-1];
          ntpos[1] = wPlaneMeanTPos[plane+1];
        }
      else if(plane>2  && plane<500 && wPlaneChargeSigCorr[plane-1]>0.0 && wPlaneChargeSigCorr[plane-3]>0.0)
        {
          nz[0] = wPlaneZPos[plane-1];
          nz[1] = wPlaneZPos[plane-3];
          ntpos[0] = wPlaneMeanTPos[plane-1];
          ntpos[1] = wPlaneMeanTPos[plane-3];
        }
      else if(plane<497  && plane>0 && wPlaneChargeSigCorr[plane+1]>0.0 && wPlaneChargeSigCorr[plane+3]>0.0)
        {
          nz[0] = wPlaneZPos[plane+1];
          nz[1] = wPlaneZPos[plane+3];
          ntpos[0] = wPlaneMeanTPos[plane+1];
          ntpos[1] = wPlaneMeanTPos[plane+3];      
        }
      else if(plane>1 && wPlaneChargeSigCorr[plane-1]>0.0)
        {stp_guess_tpos = wPlaneMeanTPos[plane-1];}
      else if(plane<497 && wPlaneChargeSigCorr[plane+1]>0.0)
        {stp_guess_tpos = wPlaneMeanTPos[plane+1];}
      else{/*lone hit*/ } 
       
      //
      //Extrapolate position
      //
      if(stp_guess_tpos<-1000.0 &&  nz[0]>0.0)
        {
          stp_guess_tpos =((ntpos[0]-ntpos[1])/(nz[0]-nz[1]))*stp_guess_z + 
            ((ntpos[1]*nz[0] - ntpos[0]*nz[1]) / (nz[0]-nz[1]));
        }
    
      //calc XY coordinates
      if(stp_guess_tpos>-1000.0)
        {
          u = (strip.View==0)?tpos:stp_guess_tpos;
          v = (strip.View==1)?tpos:stp_guess_tpos;
          y = (u+v)*0.707106781;
          x = (u-v)*0.707106781;
          
          /*
            Want to record pretrackPH, posttrackPH, first hit before vtx, last hit after end
          */
          if((strip.Trk&index)==index)
            {
              stripU = (strip.View==0)?strip.T:strip.L;
              stripV = (strip.View==1)?strip.T:strip.L;
              UCsides = ( (fabs(stripU)>(4-Xdis)) || (fabs(stripV)>(4-Xdis)) || (((fabs(stripU)+fabs(stripV))/sqrt(2.0))>(4-Xdis)) );
              UCends  = ( (strip.Plane<=Xpln) || ((strip.Plane>=(249-Xpln)) && (strip.Plane<=(249+Xpln))) || (strip.Plane>=(486-Xpln)) );
              if( UCsides || UCends ) nonfidTrkQ+=strip.Sigcorr[0]+strip.Sigcorr[1];
              totTrkQ+=strip.Sigcorr[0]+strip.Sigcorr[1];
            }

          // don't bother using very low PH hits
          if( (strip.Sigcorr[0]+strip.Sigcorr[1])<100.0 ) continue; 

          if(MyTrack->VtxPlane<MyTrack->EndPlane)
            {
              if(strip.Plane<=MyTrack->VtxPlane)
                {
                  if(VtxZdigit[2]>double(strip.Plane))
                    {
                      VtxZdigit[0]=x;
                      VtxZdigit[1]=y;
                      VtxZdigit[2]=double(strip.Plane);
                    }
                  if( ( pow(VtxRdigit[0]*VtxRdigit[0] + VtxRdigit[1]*VtxRdigit[1] ,0.5)<pow(u*u + v*v ,0.5)) )
                    {
                      VtxRdigit[0]=x;
                      VtxRdigit[1]=y;
                      VtxRdigit[2]=double(strip.Plane);
                    }
                }
              
              if(strip.Plane>=MyTrack->EndPlane)
                {
                  if(EndZdigit[2]<double(strip.Plane))
                    {
                      EndZdigit[0]=x;
                      EndZdigit[1]=y;
                      EndZdigit[2]=double(strip.Plane);
                    }
                  if(pow(EndRdigit[0]*EndRdigit[0] + EndRdigit[1]*EndRdigit[1] ,0.5)<pow(u*u + v*v ,0.5))
                    {
                      EndRdigit[0]=x;
                      EndRdigit[1]=y;
                      EndRdigit[2]=double(strip.Plane);
                    }
                }
              
            }
          else//neglect the equality case
            {
              if(strip.Plane>=MyTrack->VtxPlane)
                {
                  if(VtxZdigit[2]<double(strip.Plane))
                    {
                      VtxZdigit[0]=x;
                      VtxZdigit[1]=y;
                      VtxZdigit[2]=double(strip.Plane);
                    }
                  if(pow(VtxRdigit[0]*VtxRdigit[0] + VtxRdigit[1]*VtxRdigit[1] ,0.5)<pow(u*u + v*v ,0.5))
                    {
                      VtxRdigit[0]=x;
                      VtxRdigit[1]=y;
                      VtxRdigit[2]=double(strip.Plane);
                    }
                }
              if(strip.Plane<=MyTrack->EndPlane)
                {
                  //Trace stuff in here
                  if(EndZdigit[2]>double(strip.Plane))
                    {
                      EndZdigit[0]=x;
                      EndZdigit[1]=y;
                      EndZdigit[2]=double(strip.Plane);
                    }
                  if(pow(EndRdigit[0]*EndRdigit[0] + EndRdigit[1]*EndRdigit[1] ,0.5)<pow(u*u + v*v ,0.5))
                    {
                      EndRdigit[0]=x;
                      EndRdigit[1]=y;
                      EndRdigit[2]=double(strip.Plane);
                    }
                }
            }
        }           
    }
  //here decide on the values to use...
  MyTrack->VtxDistToEdgeDigits = 4.0 - pow(VtxRdigit[0]*VtxRdigit[0] + VtxRdigit[1]*VtxRdigit[1] ,0.5); // VtxDistToEdgeDigits
  MyTrack->VtxPlaneDigits = int(VtxZdigit[2]); // VtxPlaneDigits

  MyTrack->EndDistToEdgeDigits = 4.0 - pow(EndRdigit[0]*EndRdigit[0] + EndRdigit[1]*EndRdigit[1] ,0.5); // EndDistToEdgeDigits
  MyTrack->EndPlaneDigits = int(EndZdigit[2]); // EndPlaneDigits

  //calc the amount of track outside the fiducial volume
  if( totTrkQ>0.0 ) MyTrack->NonFidFrac=nonfidTrkQ/totTrkQ; // NonFidFrac
  return;
}

void AtmosCalculator::ShowerProperties(AtmosShower* MyShower, TClonesArray* StripList)
{
  MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::ShowerProperties(...)  " << endl;

  ShowerTrace ShwTrace(MyShower);
  MyShower->VtxForTrace = ShwTrace.Trace;
  MyShower->VtxForTraceZ = ShwTrace.TraceZ;
  MyShower->VtxRevTrace = ShwTrace.ReTrace;
  MyShower->VtxRevTraceZ = ShwTrace.ReTraceZ;
  MyShower->VtxUpTrace = ShwTrace.UpTrace;
  MyShower->VtxUpTraceZ = ShwTrace.UpTraceZ;

  // Get Shower Index
  Int_t index = MyShower->Index;

  // Closest distance to edge of detector
  MSG("AtmosCalculator",Msg::kVerbose) << "    Distance from Edge of Detector" << endl;
  double vtxdiru=MyShower->VtxDirCosU;
  double vtxdirv=MyShower->VtxDirCosV;
  double vtxdirz=MyShower->VtxDirCosZ;
  double vtxu=MyShower->VtxU;  
  double vtxv=MyShower->VtxV;  
  double vtxx=MyShower->VtxX;  
  double vtxy=MyShower->VtxY;
  double vtxz=MyShower->VtxZ;

  Double_t rmin,temp;

  rmin=4.0;
  temp=4.0-vtxu; if(temp<rmin) rmin=temp;
  temp=4.0+vtxu; if(temp<rmin) rmin=temp;
  temp=4.0-vtxv; if(temp<rmin) rmin=temp;
  temp=4.0+vtxv; if(temp<rmin) rmin=temp;
  temp=4.0-vtxx; if(temp<rmin) rmin=temp;
  temp=4.0+vtxx; if(temp<rmin) rmin=temp;
  temp=4.0-vtxy; if(temp<rmin) rmin=temp;
  temp=4.0+vtxy; if(temp<rmin) rmin=temp;
  MyShower->VtxDistToEdge = rmin;//*************

  // Highest and lowest plane number in shower
  MSG("AtmosCalculator",Msg::kVerbose) << "    Highest and Lowest Plane Number" << endl;
  int thispln;
  int minplane=-1,maxplane=-1;
  Int_t nSingleEnds(0), nDoubleEnds(0), nDigits(0);
  Double_t dQ;
  Double_t nonfidShwPH(0.), totShwPH(0.), totPlnPH(0.);
  Int_t shwplnstp[500];
  Double_t shwplnchg[500];
  Double_t plnchg[500];

  memset(shwplnstp, 0, 500*sizeof(Int_t));
  memset(shwplnchg, 0, 500*sizeof(Double_t));
  memset(plnchg, 0, 500*sizeof(Double_t));

  vector<double> UStripsU, UStripsZ, UStripsE;
  vector<double> VStripsV, VStripsZ, VStripsE;
  vector<double> AllStripsU, AllStripsV, AllStripsZ, AllStripsE;

  double UStpVtxShift, VStpVtxShift, ZStpVtxShift, RStpVtxShift;
  double DotProd, DCosHit, DCoreHit;
  double SumDCosHit(0.0), SumSqDCosHit(0.0), SumDCoreHit(0.), SumSqDCoreHit(0.0);
  
  int ListEnd = StripList->GetLast()+1;

  double ucShwStrips(0.0), ShwStrips(0.0);
  bool UCsides(false);
  bool UCends(false);
  double Xdis(0.3);
  int Xpln(5);
  double stripX, stripY, stripMaxUVXY(-1.0);

  for(int strp=0; strp<ListEnd; ++strp)
    {
      AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
      thispln=strip.Plane;
      if(thispln<0 || thispln>485) continue;

      dQ=(strip.Sigcorr[0]+strip.Sigcorr[1])/65.0; // hack a rough conversion 
      plnchg[thispln] += dQ;

      totPlnPH += dQ;

      //Everything past this point is only for strips associated with
      //the shower
      if((strip.Shw&index)!=index) continue;

      if (strip.View == 0) {
        AllStripsU.push_back(strip.T);
        AllStripsV.push_back(strip.L);

        UStripsU.push_back(strip.T);
        UStripsZ.push_back(strip.Z);
        UStripsE.push_back(dQ);
      }
      else {
        AllStripsU.push_back(strip.L);
        AllStripsV.push_back(strip.T);

        VStripsV.push_back(strip.T);
        VStripsZ.push_back(strip.Z);
        VStripsE.push_back(dQ);
      }

      AllStripsE.push_back(dQ);
      AllStripsZ.push_back(strip.Z);

      shwplnstp[thispln]++;
      shwplnchg[thispln] += dQ;

      totShwPH += dQ;
      if(strip.Ndigits==1) ++nSingleEnds;
      if(strip.Ndigits==2) ++nDoubleEnds;
      nDigits += strip.Ndigits;

      if(minplane<0||thispln<minplane) minplane=thispln;
      if(maxplane<0||thispln>maxplane) maxplane=thispln;

      if(fabs(strip.T)>stripMaxUVXY) stripMaxUVXY = fabs(strip.T);
      if(fabs(strip.L)>stripMaxUVXY) stripMaxUVXY = fabs(strip.L);
      stripX = fabs(strip.T + strip.L) / sqrt(2.0);
      if(fabs(stripX)>stripMaxUVXY) stripMaxUVXY = fabs(stripX);
      stripY = fabs(strip.T - strip.L) / sqrt(2.0);
      if(fabs(stripY)>stripMaxUVXY) stripMaxUVXY = fabs(stripY);
      UCsides = ( stripMaxUVXY>(4-Xdis) );
      UCends  = ( (strip.Plane<=Xpln) || ((strip.Plane>=(249-Xpln)) &&
                (strip.Plane<=(249+Xpln))) || (strip.Plane>=(486-Xpln)) );
      if( UCsides || UCends ) nonfidShwPH += dQ;
      if( UCsides || UCends ) ucShwStrips += 1.0;
      ShwStrips += 1.0;

      //Shift to vertex frame
      UStpVtxShift = (strip.View==0?strip.T:strip.L) - vtxu;
      VStpVtxShift = (strip.View==1?strip.T:strip.L) - vtxv;
      ZStpVtxShift = strip.Z - vtxz;
      RStpVtxShift = sqrt(UStpVtxShift*UStpVtxShift +
                          VStpVtxShift*VStpVtxShift +
                          ZStpVtxShift*ZStpVtxShift);

      if (RStpVtxShift != 0.0) {
        //Find the angle between and Shower Direction
        DotProd = fabs(UStpVtxShift*vtxdiru +
                       VStpVtxShift*vtxdirv + ZStpVtxShift*vtxdirz);
        DCosHit = DotProd / RStpVtxShift;
        if(RStpVtxShift <= DotProd) DCoreHit = 0.0;
        else DCoreHit = sqrt(RStpVtxShift*RStpVtxShift - DotProd*DotProd);

        SumDCosHit +=  DCosHit * dQ;
        SumSqDCosHit += DCosHit * DCosHit * dQ;
        SumDCoreHit +=  DCoreHit * dQ;
        SumSqDCoreHit += DCoreHit * DCoreHit * dQ;
      }
    }

  if( totShwPH>0.0 ) MyShower->NonFidFrac = nonfidShwPH/totShwPH;

  double ShwMeanDCosHit = SumDCosHit / totPlnPH;
  double ShwRMSDCosHit = (SumSqDCosHit/totPlnPH) - (ShwMeanDCosHit*ShwMeanDCosHit);
  if(ShwRMSDCosHit<=0) ShwRMSDCosHit = 0.0;
  else ShwRMSDCosHit = sqrt(ShwRMSDCosHit);

  MyShower->MeanDCosHit = ShwMeanDCosHit;
  MyShower->RMSDCosHit = ShwRMSDCosHit;

  double ShwMeanDCoreHit = SumDCoreHit / totPlnPH;
  double ShwRMSDCoreHit = (SumSqDCoreHit/totPlnPH) - (ShwMeanDCoreHit*ShwMeanDCoreHit);
  if(ShwRMSDCoreHit<=0) ShwRMSDCoreHit = 0.0;
  else ShwRMSDCoreHit = sqrt(ShwRMSDCoreHit);

  MyShower->MeanDCoreHit = ShwMeanDCoreHit;
  MyShower->RMSDCoreHit = ShwRMSDCoreHit;

  //Solve the 3D UVZ MOI
  ShowerMOI MOIUVZ(AllStripsU, AllStripsV, AllStripsZ, AllStripsE);
  MyShower->MOIUVZEigenValue0 = MOIUVZ.EigenValues[0];
  MyShower->MOIUVZEigenValue1 = MOIUVZ.EigenValues[1];
  MyShower->MOIUVZEigenValue2 = MOIUVZ.EigenValues[2];
  for (int i=0; i<3; i++) {
    MyShower->MOIUVZEigenVector0[i] = MOIUVZ.EigenVectors[0][i];
    MyShower->MOIUVZEigenVector1[i] = MOIUVZ.EigenVectors[1][i];
    MyShower->MOIUVZEigenVector2[i] = MOIUVZ.EigenVectors[2][i];
  }

  ShowerMOI MOIUZ(UStripsU, UStripsZ, UStripsE);
  MyShower->MOIUZEigenValue0 = MOIUZ.EigenValues[0];
  MyShower->MOIUZEigenValue1 = MOIUZ.EigenValues[1];
  for (int i=0; i<2; i++) {
    MyShower->MOIUZEigenVector0[i] = MOIUZ.EigenVectors[0][i];
    MyShower->MOIUZEigenVector1[i] = MOIUZ.EigenVectors[1][i];
  }

  ShowerMOI MOIVZ(VStripsV, VStripsZ, VStripsE);
  MyShower->MOIVZEigenValue0 = MOIVZ.EigenValues[0];
  MyShower->MOIVZEigenValue1 = MOIVZ.EigenValues[1];
  for (int i=0; i<2; i++) {
    MyShower->MOIVZEigenVector0[i] = MOIVZ.EigenVectors[0][i];
    MyShower->MOIVZEigenVector1[i] = MOIVZ.EigenVectors[1][i];
  }

  MyShower->MinPlaneNumber = minplane;
  MyShower->MaxPlaneNumber = maxplane;
            
  MyShower->Ndigits = nDigits;
  MyShower->NstripsSingleEnded = nSingleEnds;
  MyShower->NstripsDoubleEnded = nDoubleEnds;

  MyShower->AssocShwPH = totShwPH;
  MyShower->AssocShwPHfrac = totShwPH / totPlnPH;
  
  int uPlanes(0),vPlanes(0);
  double nPlanes(0.);
  int maxplnstp(0); double maxplnchg(0.0); 
  int sumstp(0); double sumsqstp(0.0); 
  double sumchg(0); double sumsqchg(0.0); 
  int ntrkplns(0);
  for(int ipln=minplane; ipln<=maxplane; ipln++)
    {
      if(shwplnstp[ipln] == 0) continue;
      if((ipln/249)^(ipln%2)) uPlanes++;
      else vPlanes++;
      nPlanes += 1.0;

      //Count up the tracklike planes
      if(plnchg[ipln]>0.0 && plnchg[ipln]<80.0) {
        if(shwplnstp[ipln] <=2) ntrkplns++;
      }

      if(shwplnstp[ipln] > maxplnstp) maxplnstp = shwplnstp[ipln];
      if(shwplnchg[ipln] > maxplnchg) maxplnchg = shwplnchg[ipln];

      sumstp += shwplnstp[ipln];
      sumsqstp += shwplnstp[ipln]*shwplnstp[ipln];

      sumchg += shwplnchg[ipln];
      sumsqchg += shwplnchg[ipln]*shwplnchg[ipln];
    }

  MyShower->NplanesUview = uPlanes;
  MyShower->NplanesVview = vPlanes;
  if(uPlanes+vPlanes>0) MyShower->UVassymetry = abs(uPlanes-vPlanes)/(0.5*(uPlanes+vPlanes));
  MyShower->TrkLikePlanes = ntrkplns;

  MyShower->MaxStripsInPlane = maxplnstp;
  MyShower->MaxChargeInPlane = maxplnchg;

  double ShwMeanStripsPerPlane = sumstp/nPlanes;
  double ShwRMSStripsPerPlane = (sumsqstp/nPlanes) - ShwMeanStripsPerPlane*ShwMeanStripsPerPlane;
  if(ShwRMSStripsPerPlane <= 0.0) ShwRMSStripsPerPlane = 0.0;
  else ShwRMSStripsPerPlane = sqrt(ShwRMSStripsPerPlane);
  MyShower->MeanStripsPerPlane = ShwMeanStripsPerPlane;
  MyShower->RMSStripsPerPlane = ShwRMSStripsPerPlane;
    
  double ShwMeanChargePerPlane = sumchg/nPlanes;
  double ShwRMSChargePerPlane = (sumsqchg/nPlanes) - ShwMeanChargePerPlane*ShwMeanChargePerPlane;
  if(ShwRMSChargePerPlane <= 0.0) ShwRMSChargePerPlane = 0.0;
  else ShwRMSChargePerPlane = sqrt(ShwRMSChargePerPlane);
  MyShower->MeanChargePerPlane = ShwMeanChargePerPlane;
  MyShower->RMSChargePerPlane = ShwRMSChargePerPlane;
    
  return;
}
bool AtmosCalculator::CalculateTrace(double* m, double* c, double* Coord, 
                                    double* Trace)
{ 
  MSG("AtmosCalculator",Msg::kDebug) << "CalculateTrace" << endl;

  bool TraceSuccess=false;
  double TempTrace[3]={0.,0.,0.}; double TempSum=0.;
  double TraceLength(1.e99);

  for(int i=0; i<8; ++i) { // Consider the eight sides in turn
    TempSum=0.;

    DetectorSides(m,c,TempTrace, i);

    for(int j=0; j<3; ++j) {TempSum+=pow(Coord[j]-TempTrace[j],2);}
    TempSum=pow(TempSum,0.5);

    if(TempSum<TraceLength) {
      for(int j=0; j<3; ++j) {Trace[j]=TempTrace[j]-Coord[j];}
      TraceLength=TempSum;
      TraceSuccess=true;
    }
  }

  return TraceSuccess;
}



void AtmosCalculator::DetectorSides(double* m, double* c, double* position, int side)
{
  position[2]=-1.e99;
  
  switch(side) {
  case 0: if(m[0]!=-m[1]) {position[2] = (pow(32.,0.5)-c[0]-c[1])/(m[0]+m[1]);}  break;
  case 1: if(m[0]!=0.)    {position[2] = (4.-c[0])/m[0];}                        break;
  case 2: if(m[0]!=m[1])  {position[2] = (pow(32.,0.5)-c[0]+c[1])/(m[0]-m[1]);}  break;
  case 3: if(m[1]!=0.)    {position[2] = (-4.-c[1])/m[1];}                       break;
  case 4: if(m[0]!=-m[1]) {position[2] = (-pow(32.,0.5)-c[0]-c[1])/(m[0]+m[1]);} break;
  case 5: if(m[0]!=0.)    {position[2] = (-4.-c[0])/m[0];}                       break;
  case 6: if(m[0]!=m[1])  {position[2] = (-pow(32.,0.5)-c[0]+c[1])/(m[0]-m[1]);} break;
  case 7: if(m[1]!=0.)    {position[2] = (4.-c[1])/m[1];}                        break;
  default: position[2] = -1.e99; break;
  }
  
  if(position[2]!=-1.e99) {
    position[0]=m[0]*position[2]+c[0];
    position[1]=m[1]*position[2]+c[1];  
  }
  else {position[0]=-1.e99; position[1]=-1.e99;}
}

---