////////////////////////////////////////////////////////////////////////
//
//  JSFCDCTracks
//
//  Utilities for CDC Tracks 
//
// (Update)
// 2-May-1999 A.Miyamoto  Bug in MovePivot is fixed.
// 5-May-1999 A.Miyamoto  Add MovePivotToIP(...) and AddMSError(...)
// 14-May-1999 A.Miyamoto Avoid negative costh in MovePivotToIP.
// 17-May-1999 A.Miyamoto Put "delete helix" in MovePivotToIP to avoid memory leak
//                      problem.  Thanks Hyunwoo Kim for pointed out the problem.
//  1-Feb-2000 A.Miyamoto Increase limitter of #VTX hits in AddVTXHits from 10 to kMaxVTXAssoc. 
//                        and avoid invalid memory access when hits exceds kMaxVTXAssoc.
// 
//$Id: JSFCDCTrack.cxx,v 1.20 2003/02/14 07:53:52 miyamoto Exp $
//
////////////////////////////////////////////////////////////////////////

#include "JSFSteer.h"
#include "JSFQuickSimParam.h"
#include "JSFCDCTrack.h"
#include "JSFHelicalTrack.h"

ClassImp(JSFCDCTrack)

extern "C" {
extern void utrkmv_(Int_t *lnxhlx, Float_t helxin[], 
		    Float_t xp[], Float_t helxot[]);
};

//_____________________________________________________________________________
 JSFCDCTrack::JSFCDCTrack() : fNVTX(0)
{
}

//_____________________________________________________________________________
 JSFCDCTrack::JSFCDCTrack(Int_t itrkp[])
{
  // Make a JSFCDCTrack class from a data from Production:CDC;Track_Parameter


  fGenID=itrkp[56];
  Float_t *trk=(Float_t*)&itrkp[0];
  fP[0]=trk[0];  fP[1]=trk[1];  fP[2]=trk[2];
  fE=trk[3];
  fX[0]=trk[4]; fX[1]=trk[5]; fX[2]=trk[6];

  //  if( TMath::Abs(trk[6]) > 500.0 ) {
  //    printf(" Strange track in CDCTrack bank... closest approach is x,y,z=%g,%g,%gn",
  //	   trk[4],trk[5],trk[6]);
  //  }

  fCharge = itrkp[8];
  Int_t i;
  for(i=0;i<5;i++){ fHelix[i]=trk[10+i] ;}
  fPivot[0]=trk[15] ; fPivot[1]=trk[16]; fPivot[2]=trk[17] ;
  Short_t *ndf=(Short_t*)&itrkp[50];
  fNDF=*ndf;

  fPosAtEMC[0]=-1.0; fPosAtEMC[1]=0.0 ; fPosAtEMC[2]=0.0 ;
  fEPosAtEMC[0]=0.0 ;    fEPosAtEMC[1]=0.0 ;  

#ifdef R__ACC
  memcpy(&fError[0],&itrkp[18],120);
#else
  Double_t *err=(Double_t*)&itrkp[18];
  for(i=0;i<15;i++){ fError[i]=err[i]; }
#endif

  fNVTX = 0;
  for(i=0;i<kMaxVTXAssoc;i++){ fVTXHits[i]=NULL; }

}


//_____________________________________________________________________________
 JSFCDCTrack::JSFCDCTrack(Float_t trkf[], Double_t trkd[])
{
  // Make a JSFCDCTrack class from a data from Production:CDC;Track_Parameter
  fP[0]=trkf[0]; fP[1]=trkf[1]; fP[2]=trkf[2];
  fE=trkf[3];
  fX[0]=trkf[4]; fX[1]=trkf[5] ; fX[2]=trkf[6];
  fCharge=(Int_t)trkf[7];  fGenID=(Int_t)trkf[8];
  Int_t i;
  for(i=0;i<5;i++){ fHelix[i]=trkf[9+i] ;}
  fPivot[0]=trkf[14] ; fPivot[1]=trkf[15]; fPivot[2]=trkf[16] ;
  fNDF=(Int_t)trkf[17];
  fPosAtEMC[0]=trkf[18]; fPosAtEMC[1]=trkf[19] ; fPosAtEMC[2]=trkf[20] ;
  fEPosAtEMC[0]=trkf[21] ;    fEPosAtEMC[1]=trkf[22];  
  for(i=0;i<15;i++){ fError[i]=trkd[i]; }

  fNVTX = 0;
  for(i=0;i<kMaxVTXAssoc;i++){ fVTXHits[i]=NULL; }

}


//_____________________________________________________________________________
 JSFCDCTrack::JSFCDCTrack(JSFCDCTrack& t)
{
  // Make a JSFCDCTrack class from a data from Production:CDC;Track_Parameter

  fE=t.fE;  fCharge=t.fCharge ; fGenID=t.fGenID; fNDF=t.fNDF;
  Int_t i;
  for(i=0;i<3;i++){
    fP[i]=t.fP[i] ; fX[i]=t.fX[i]; 
    fPivot[i]=t.fPivot[i] ;  fPosAtEMC[i]=t.fPosAtEMC[i];
  }
  for(i=0;i<5;i++){ fHelix[i]=t.fHelix[i]; }
  for(i=0;i<15;i++){ fError[i]=t.fError[i]; }
  fNVTX=t.fNVTX;
  for(i=0;i<fNVTX;i++){ fVTXHits[i]=t.fVTXHits[i]; }

}


//_____________________________________________________________________________
 JSFCDCTrack::~JSFCDCTrack()
{
}

//_____________________________________________________________________________
 void JSFCDCTrack::AddVTXHit(JSFVTXHit *v)
{ 
  static Int_t ntrue=0;
  if( fNVTX >= kMaxVTXAssoc ) { 
    ntrue++;
    printf("Too many VTX hits linked to CDC tracks in JSFCDCTrack::AddVTXHit.n");
    printf("Track Momentum (E,Px,Py,Pz)=(%g,%g,%g,%g)n",fE,fP[0],fP[1],fP[2]);
    printf("This is %d-th hit, but not saved in the buffer.n",ntrue);
    return;
  }
  fVTXHits[fNVTX]=v; 
  fNVTX++ ; 
  ntrue=fNVTX;
}


//_____________________________________________________________________________
 void JSFCDCTrack::Print()
{
  printf(" CDCTrack parameter n");
  printf(" fHelix=%g %g %g %g %g n",
	 fHelix[0],fHelix[1],fHelix[2],fHelix[3],fHelix[4]);
  printf(" fPivot[0:2]=%g %g %gn",fPivot[0],fPivot[1],fPivot[2]);
  printf(" fError[0,2,5,9,14]=%g %g %g %g %gn",
	 fError[0],fError[2],fError[5],fError[9],fError[14]);

}

//_____________________________________________________________________________
 void JSFCDCTrack::GetErrorMatrix(Double_t trkd[])
{
  Int_t i;
  for(i=0;i<15;i++){ trkd[i]=fError[i]; }
}

//_____________________________________________________________________________
 void JSFCDCTrack::GetTrackParameter(Float_t trkf[])
{
  // Get the Track parameter with the SIMDST format.
  
  trkf[0]=fP[0]; trkf[1]=fP[1];  trkf[2]=fP[2];
  trkf[3]=fE;
  trkf[4]=fX[0]; trkf[5]=fX[1]; trkf[6]=fX[2];
  trkf[7]=(Float_t)fCharge;   trkf[8]=(Float_t)fGenID;
  Int_t i;
  for(i=0;i<5;i++){ trkf[9+i]=fHelix[i] ;}
  trkf[14]=fPivot[0] ; trkf[15]=fPivot[1]; trkf[16]=fPivot[2] ;
  trkf[17]=fNDF;
  trkf[18]=fPosAtEMC[0];    trkf[19]=fPosAtEMC[1];  trkf[20]=fPosAtEMC[2];
  trkf[21]=fEPosAtEMC[0];   trkf[22]=fEPosAtEMC[1];
}

//_____________________________________________________________________________
 void JSFCDCTrack::GetPositionAtEMC(Float_t pos[3], Float_t poserr[2])
{
  // Return position and its error at EMC surface.

  for(Int_t i=0;i<3;i++){ pos[i]=fPosAtEMC[i]; }
  for(Int_t i=0;i<2;i++){ poserr[i]=fEPosAtEMC[i]; }

}


//_____________________________________________________________________________
 void JSFCDCTrack::SetPositionAtEMC(Float_t xp[])
{
  // Calculate position and its error of CDC Track at EMC entrance.
  // and set it to the data member of the class.

  Int_t j;
  const Int_t lenhelix=40;
  Float_t helixini[lenhelix];
  for(j=0;j<5;j++){ helixini[j]=fHelix[j]; }  
  helixini[5]=fPivot[0];  helixini[6]=fPivot[1];  helixini[7]=fPivot[2]; 
  Double_t *err=(Double_t*)&helixini[8];
  Int_t map[15]={0, 1,5,  2,6,9,  3,7,10,12,  4,8,11,13,14 };
   //  LC-packing   ------> Packing used by utrkmv.
   //  0                       0
   //  1   2                   1  5
   //  3   4   5               2  6   9
   //  6   7   8   9           3  7  10  12
   // 10  11  12  13  14       4  8  11  13  14
  for(j=0;j<15;j++){err[map[j]]=fError[j];}

  Float_t helixnew[lenhelix];
  Int_t il=lenhelix;
  utrkmv_(&il, helixini, xp, helixnew);
  ExtrapolateErrorAtEMC(helixnew, fPosAtEMC, fEPosAtEMC);
 
}


//_____________________________________________________________________________
 void JSFCDCTrack::ExtrapolateErrorAtEMC(Float_t helix[], Float_t x[], Float_t dx[])
{
// 
//   Calculate Pivot's polar coordinate (theta, phi) and its error.
//   (Input)
//       Float_t helix[] ; helix parameter and error matrix
//   (Output)
//       Float_t x[0] ; radius of the pivot ( distance from the beam line)
//               x[1] ; theta of pivot
//               x[2] ; phi 
//       Float_t dx[0] ; delta-costh
//               dx[1] ; delta-phi

   Double_t dr=helix[0];  Double_t f0=helix[1];  Double_t dz=helix[3];
   struct XYZ { Double_t x; Double_t y; Double_t z; } xyz;

   Double_t csf0 = TMath::Cos(f0);
   Double_t snf0 = TMath::Sin(f0);
   xyz.x=helix[5] + dr*csf0;  
   xyz.y=helix[6] + dr*snf0;  
   xyz.z=helix[7] + dz ;

   Double_t r2   = xyz.x*xyz.x + xyz.y*xyz.y;
   Double_t r    = TMath::Sqrt(r2);
   Double_t th   = TMath::ATan2(r,xyz.z);
   Double_t phi  = TMath::ATan2(xyz.y,xyz.x);
   
   Double_t csphi = xyz.x/r;
   Double_t snphi = xyz.y/r;
   
   Double_t dphidr = (-xyz.y*csf0 + xyz.x*snf0)/r2;
   Double_t dphif0 = ( xyz.y*snf0 + xyz.x*csf0)/r2*dr;

   Double_t R2 = xyz.x*xyz.x + xyz.y*xyz.y + xyz.z*xyz.z; 
   Double_t R  = TMath::Sqrt(R2);
   Double_t dthdr = xyz.z/R2*(csphi*csf0+snphi*snf0);
   Double_t dthf0 = xyz.z/R2*(-csphi*snf0 + snphi*csf0 )*dr;
   Double_t dthdz = -TMath::Sin(th)/R;

   const Int_t iPdrdr=0;  
   const Int_t iPdrf0=1;  const Int_t iPf0f0=5;  
   const Int_t iPdrdz=3;  const Int_t iPf0dz=7;  const Int_t iPdzdz=12;
   //  LC-packing   ------> Packing used by utrkmv.
   //  0                       0
   //  1   2                   1  5
   //  3   4   5               2  6   9
   //  6   7   8   9           3  7  10  12
   // 10  11  12  13  14       4  8  11  13  14
  


#ifdef R__ACC
   Double_t erm[15];
   memcpy(erm,&helix[8],120);
#else
   Double_t *erm = (Double_t*)&helix[8];
#endif

   Double_t dphi2 = dphidr*dphidr/2.*erm[iPdrdr] 
                  + dphidr*dphif0*erm[iPdrf0] + dphif0*dphif0/2.0*erm[iPf0f0];
   Double_t dth2  = dthdr*dthdr/2.*erm[iPdrdr] 
                  + dthdr*dthf0*erm[iPdrf0] + dthdr*dthdz*erm[iPdrdz]
                  + dthf0*dthf0/2.0*erm[iPf0f0] +dthf0*dthdz*erm[iPf0dz]
                  + dthdz*dthdz/2.0*erm[iPdzdz];

   if( dth2 < 0.0 ) {
     printf(" Invalid dth2 .. \n");
     
     printf(" dth2=%g",dth2);
     printf(" dphi2 =%g dr,f0,dz=%g %g %g\n",dphi2, dr,f0,dz);
     printf(" helix[0:4]=%g %g %g %g %g\n",helix[0],helix[1],helix[2],
	    helix[3],helix[4]);
     printf(" helix[5:7]=%g %g %g\n",helix[5],helix[6],helix[7]);
     printf(" erm=%g %g %g %g %g %g\n",erm[iPdrdr],erm[iPdrf0],
	    erm[iPdrdz],erm[iPf0f0],erm[iPf0dz],erm[iPdzdz]);
     printf(" dthdr=%g",dthdr);
     printf(" dthf0=%g",dthf0);
     printf(" dthdz=%g",dthdz);
     printf("\n");

    Double_t t1  = dthdr*dthdr/2.*erm[iPdrdr] ;
    Double_t t2  = dthdr*dthf0*erm[iPdrf0] ;
    Double_t t3  = dthdr*dthdz*erm[iPdrdz] ;
    Double_t t4  = dthf0*dthf0/2.0*erm[iPf0f0] ;
    Double_t t5  = dthf0*dthdz*erm[iPf0dz];
    Double_t t6  = dthdz*dthdz/2.0*erm[iPdzdz];

    printf(" t1:t6=%g %g %g %g %g %g \n",t1,t2,t3,t4,t5,t6);

     printf("\n");
   }

   x[0] = r;
   x[1] = th;
   x[2] = phi;
   dx[0] = TMath::Sqrt(dth2);
   dx[1] = TMath::Sqrt(dphi2);

}


//_____________________________________________________________________________
 void JSFCDCTrack::MovePivot(Float_t pivot[], Float_t bfield)
{
  // Current helix parameter is expressed in terms of new pivot.

  Double_t ptor=333.5640952*10.0/bfield;
  Double_t x2pid = 2.0*TMath::Pi();
  Double_t xpid  = TMath::Pi();
  Double_t dr  = fHelix[0];
  Double_t fi0 = fHelix[1];
  Double_t cpa = fHelix[2];
  Double_t dz  = fHelix[3];
  Double_t tnl = fHelix[4];
  Double_t x0  = fPivot[0];
  Double_t y0  = fPivot[1];
  Double_t z0  = fPivot[2];
  Double_t xv  = pivot[0];
  Double_t yv  = pivot[1];
  Double_t zv  = pivot[2];
  //
  // Transform helix parameters
  //
  Double_t r   = ptor/cpa;
  Double_t rdr = r + dr;
  Double_t fip = (fi0+2*x2pid)-x2pid*((Double_t)((Int_t)((fi0+2*x2pid)/x2pid)));
  Double_t csf0 = TMath::Cos(fip);
  Double_t snf0 = TMath::Sqrt( TMath::Max(0.0, (1.0-csf0)*(1.0+csf0)) );
  if( fip > xpid ) snf0 = -snf0 ;

  Double_t xc  = x0 + rdr*csf0 ;
  Double_t yc  = y0 + rdr*snf0 ;
  Double_t csf = ( xc-xv)/r;
  Double_t snf = (yc-yv)/r;
  Double_t anrm = 1.0/TMath::Sqrt(csf*csf+snf*snf);
           csf  *= anrm;
	   snf  *= anrm;
  Double_t csfd  = csf*csf0 + snf*snf0;
  Double_t snfd  = snf*csf0 - csf*snf0;
           fip   = TMath::ATan2(snf, csf);
  Double_t fid   = (fip-fi0 + 4*x2pid)
                 - x2pid*((Double_t)(Int_t)((fip-fi0 + 4*x2pid)/x2pid));
  if( fid > xpid ) fid -= x2pid;
  Double_t drp   = (x0+dr*csf0+r*(csf0-csf)-xv)*csf
                 + (y0+dr*snf0+r*(snf0-snf)-yv)*snf;
  Double_t dzp   = z0 + dz - r*tnl*fid - zv;
  //C--
  //C  Calculate @AP/@A.
  //C     AP = ( DRP, FIP, CPA, DZP, TNL )
  //C     A  = ( DR , FI0, CPA, DZ , TNL )
  //C--
  Double_t rdrpr = 1.0/(r+drp);
  Double_t rcpar = r/cpa;
  Double_t dapda[5][5];
  memset(dapda,0,200);
  
  dapda[0][0] = csfd ; 
  dapda[1][0] = rdr*snfd;  
  dapda[2][0] = rcpar*(1.0-csfd);   
   
  dapda[0][1] = - rdrpr*snfd;
  dapda[1][1] =   rdr*rdrpr*csfd;
  dapda[2][1] =   rcpar*rdrpr*snfd;
   
  dapda[2][2] =   1.0;
   
  dapda[0][3] =   r*rdrpr*tnl*snfd;
  dapda[1][3] =   r*tnl*(1.0-rdr*rdrpr*csfd);
  dapda[2][3] =   rcpar*tnl*(fid-r*rdrpr*snfd);
  dapda[3][3] =   1.0;
  dapda[4][3] = - r*fid;
   
  dapda[4][4] =   1.0;

  //C--
  //C  Copy error matrix to EEP and symmetrize it into EE.
  //C--
  Double_t ee[5][5];
  Int_t i,j,n;
  n=0;
  for(i=0;i<5;i++){ for(j=0;j<=i;j++) {
    ee[i][j]=fError[n];
    ee[j][i]=ee[i][j];
    n++;
  }}
  //C--
  //C  Transform error matrix EEP to that of XP.
  //C--
  n=0;
  for(i=0;i<5;i++){ for(j=0;j<=i;j++) {
    fError[n]=0.0;
    Int_t k,l;
    for(k=0;k<5;k++){ for(l=0;l<5;l++){
      fError[n]+= dapda[k][i]*ee[l][k]*dapda[l][j];
    }}
    n++;
  }}
  //C--
  //C  Fill HELXOT array.
  //C--
  fHelix[0]=drp;
  fHelix[1]=fip;
  fHelix[2]=cpa;
  fHelix[3]=dzp;
  fHelix[4]=tnl;
  fPivot[0]=xv;
  fPivot[1]=yv;
  fPivot[2]=zv;

  return;

}

//______________________________________________________________________________
 void JSFCDCTrack::AddMSError(Float_t xrad, Float_t deltakappa)
{
  // Increase the error matrix to include the effect of 
  // the multiple scattering in the matterinal of radiation length xrad.
  //
  // deltakappa is subtracted from Kappa of helix parameter for 
  // energy loss correction

  Double_t tnlsq=fHelix[4]*fHelix[4];
  Double_t tnlsqone=1.0+tnlsq;
  Double_t pt=1.0/TMath::Abs(fHelix[2]);
  Double_t p =pt*TMath::Sqrt(tnlsqone);
  Double_t radx=xrad;
  Double_t sigms=0.0141*(1.0+TMath::Log10(radx)/9.0)*TMath::Sqrt(radx)/p;
  Double_t sigmsq=sigms*sigms;

  fError[2] = fError[2] +sigmsq*tnlsqone;
  fError[5] = fError[5] +sigmsq*(fHelix[2]*fHelix[2]*tnlsq);
  fError[12]= fError[12]+sigmsq*fHelix[2]*fHelix[4]*tnlsqone;
  fError[14]= fError[14]+sigmsq*tnlsqone*tnlsqone;

  //C .. E(2,2)=EDAT(3), E(3,3)=Edat(6), E(3,5)=Edat(13), E(5,5)=Edat(15)

  fHelix[2]-=deltakappa;

}

//____________________________________________________________________________
 Bool_t JSFCDCTrack::MovePivotToIP(JSFQuickSimParam *spar)
{
  // Move Pivot of Track parameter to IP
  // This program assumes that pivot of the track parameter is at the
  // first layer of VTX.  Procedure is 
  //   (1) Include multiple-scattering effect in first layer of VTX, 
  //   (2) Move pivot to beam pipe,
  //   (3) Include multiple-scatteing efect in the beam pipe.
  //   (4) Move pivot to IP.
  // Geometry information is obtained from JSFQuickSimParam
  
  // First make sure input parameter is correct.

  JSFHelixParameter hp=GetHelix();
  Double_t pnorm=TMath::Sqrt(1.0+hp.tanl*hp.tanl);
  Double_t rnow=TMath::Sqrt(hp.pivot.x*hp.pivot.x+hp.pivot.y*hp.pivot.y);
  Float_t field=spar->GetBField();

  // Make sure Pivot is near 1st layer of VTX and dr is small enough
  if( rnow > spar->GetVTXRadius(2) - 0.01 ) { 
    printf("Warning in JSFCDCTrack::MovePivotToIP() ..");
    printf("Pivot is far away of 1st layer of VTX.  Move to the 1st layer of VTX first.n");
    printf(" Radius of pivot is %g",rnow);
    printf(" radius of 1st VTX layer is %gn",spar->GetVTXRadius(1));
    return kFALSE; 
  }
  else if(  TMath::Abs(hp.dr) > 0.01 || TMath::Abs(rnow - spar->GetVTXRadius(1)) > 0.001 ) {
    //    printf("Warning in JSFCDCTrack::MovePivotToIP() ..");
    //    printf("Pivot is moved to the first layer of VTX since it is not on yet.");
    //    printf(" or dr is too large.n");
    //    printf(" dr=%g",hp.dr);
    //    printf(" rnow=%gn",rnow);
    JSFHelicalTrack ht=GetHelicalTrack();
    ht.SetBfield(field);
    JSF3DV pos=ht.GetCoordinate(0.0);
    JSFRPhiZ ref(pos);
    ref.r=spar->GetVTXRadius(1);
    Double_t defang=0.0;
    Bool_t st=ht.IntersectWithCylinder(ref, defang);
    if( !st ) { return kFALSE; }
    JSF3DV piv=ht.GetCoordinate(defang);
    Float_t pivot[3]={piv.x, piv.y, piv.z};
    MovePivot(pivot, field);

    hp=GetHelix();
    pnorm=TMath::Sqrt(1.0+hp.tanl*hp.tanl);
    rnow=TMath::Sqrt(hp.pivot.x*hp.pivot.x+hp.pivot.y*hp.pivot.y);
    //    printf("After modification rnow=%g",rnow);
    //    printf(" pivot is %g,%gn",hp.pivot.x,hp.pivot.y);
    //    Print();

  }

  Double_t costh=(-hp.pivot.x*TMath::Sin(hp.phi0) + 
		  hp.pivot.y*TMath::Cos(hp.phi0)) / (rnow*pnorm) ;
  if( costh < 0.0 ) {
    /*
	printf("Warning in JSFCDCTrack::MovePivotToIP()...vtx1");
	printf(" costh=%g <=0, though it should be positiven",costh);
	printf("Abs(costh) is used instead..n");	
	printf(" rnow=%g  pnorm=%gn",rnow, pnorm);
        printf(" pivot=(%g,%g) (direction x, direction y)=(%g,%g)n",hp.pivot.x,
               hp.pivot.y,-TMath::Sin(hp.phi0), TMath::Cos(hp.phi0));
	JSFHelicalTrack ht=GetHelicalTrack();
	ht.SetBfield(field);
        JSF3DV pnow=ht.GetMomentum(0.0);
        JSF3DV xnow=ht.GetCoordinate(0.0);
        printf("Momentum=(%g,%g,%g)",pnow.x,pnow.y,pnow.z);
        printf("Position=(%g,%g,%g)",xnow.x,xnow.y,xnow.z);
	printf("Radius=%g",TMath::Sqrt(xnow.x*xnow.x + xnow.y*xnow.y));
        printf("n");
	ht.Print();
    */
	costh= -costh;
  }
  if( costh < 0.001 ) { costh=0.001; }  // Too small costh is unnatural 
  Float_t rad1=spar->GetVTXThickness(1)/costh;
  AddMSError(rad1);  // Include MS at the first layer of VXT

  JSFHelicalTrack helix=GetHelicalTrack();
  helix.SetBfield(field);
  Double_t rcyl=spar->GetVTXRadius(0);
  Double_t zcyl=spar->GetVTXZplus(0);
  Double_t phi0, phi1;
  Int_t maxloop=5;
  Int_t ncros=helix.OriginToCylinder(rcyl, zcyl, phi0, phi1, maxloop);
  if( ncros != 0 ) {
    /*
    printf("Warning in JSFCDCTrack::MovePivotToIP(...)  .. Track does not intersect");
    printf(" with beam pipe.n");
    printf("Track parameter is not changed.n");
        JSF3DV pnow=helix.GetMomentum(0.0);
        JSF3DV xnow=helix.GetCoordinate(0.0);
        printf("Momentum=(%g,%g,%g)",pnow.x,pnow.y,pnow.z);
        printf("Position=(%g,%g,%g)",xnow.x,xnow.y,xnow.z);
	printf("Radius of position =%g",TMath::Sqrt(xnow.x*xnow.x + xnow.y*xnow.y));
        printf("n");
	helix.Print();
    */
    return kFALSE;
  }
  //
  JSF3DV piv=helix.GetCoordinate(phi1);
  Float_t pivot[3]={piv.x, piv.y, piv.z};
  // Move to Beam pipe.
  MovePivot(pivot, field);

  
  //  hp=helix.GetHelixParameter();
  hp=GetHelix();
  pnorm=TMath::Sqrt(1.0+hp.tanl*hp.tanl);
  costh=(-pivot[0]*TMath::Sin(hp.phi0)+pivot[1]*TMath::Cos(hp.phi0))
                / (rcyl*pnorm) ;
  if( costh < 0.0 ) {
    /*
	printf("Warning in JSFCDCTrack::MovePivotToIP()...ip");
	printf(" costh=%g <=0, though it should be positiven",costh);
	printf("Abs(costh) is used instead..n");	
    */
	costh= -costh;
  }
  if( costh < 0.001 ) { costh=0.001; }  // Too small costh is unnatural 
  Float_t rad0=spar->GetVTXThickness(0)/costh;
  AddMSError(rad0);
  Float_t ip[3]={0.0, 0.0, 0.0};
  MovePivot(ip, field);

  return kTRUE;

}

#if __ROOT_FULLVERSION__ >= 30000
//____________________________________________________________________________
 void JSFCDCTrack::Streamer(TBuffer &R__b)
{
   // Stream an object of class JSFCDCTrack.

   if (R__b.IsReading()) {
     UInt_t R__s, R__c;
     Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
     if( R__v > 1 ) {
       JSFCDCTrack::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
     }
     else {
       //=== Process old versionsbefore automatic schema evolution
       TObject::Streamer(R__b);
       R__b.ReadStaticArray(fP);
       R__b >> fE;
       R__b.ReadStaticArray(fX);
       R__b >> fCharge;
       R__b >> fGenID;
       R__b.ReadStaticArray(fHelix);
       R__b.ReadStaticArray(fPivot);
       R__b.ReadStaticArray(fError);
       R__b >> fNDF;
       R__b.ReadStaticArray(fPosAtEMC);
       R__b.ReadStaticArray(fEPosAtEMC);
       R__b.CheckByteCount(R__s, R__c, JSFCDCTrack::IsA());
     }
     //      R__b >> fNVTX;
     fNVTX=0; // fNVTX is set to zero untill VTX address is set in SetPointers()
   } else {
     JSFCDCTrack::Class()->WriteBuffer(R__b, this);
   }
}


#else
//____________________________________________________________________________
 void JSFCDCTrack::Streamer(TBuffer &R__b)
{
   // Stream an object of class JSFCDCTrack.

   if (R__b.IsReading()) {
      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
      TObject::Streamer(R__b);
      R__b.ReadStaticArray(fP);
      R__b >> fE;
      R__b.ReadStaticArray(fX);
      R__b >> fCharge;
      R__b >> fGenID;
      R__b.ReadStaticArray(fHelix);
      R__b.ReadStaticArray(fPivot);
      R__b.ReadStaticArray(fError);
      R__b >> fNDF;
      R__b.ReadStaticArray(fPosAtEMC);
      R__b.ReadStaticArray(fEPosAtEMC);
      //      R__b >> fNVTX;
      fNVTX=0; // fNVTX is set to zero untill VTX address is set in SetPointers()
   } else {
      R__b.WriteVersion(JSFCDCTrack::IsA());
      TObject::Streamer(R__b);
      R__b.WriteArray(fP, 3);
      R__b << fE;
      R__b.WriteArray(fX, 3);
      R__b << fCharge;
      R__b << fGenID;
      R__b.WriteArray(fHelix, 5);
      R__b.WriteArray(fPivot, 3);
      R__b.WriteArray(fError, 15);
      R__b << fNDF;
      R__b.WriteArray(fPosAtEMC, 3);
      R__b.WriteArray(fEPosAtEMC, 2);
      //  R__b << fNVTX;
   }
}


#endif





ROOT page - Home page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.