#include <PndLmdLinFitTask.h>

Inheritance diagram for PndLmdLinFitTask:

Public Member Functions
	PndLmdLinFitTask ()

	PndLmdLinFitTask (TString tTCandBranchName, TString tRecoBranchName, TString tOutputBranchName="LMDPndTrack", TString tOutputFolder="PndLmd")

virtual	~PndLmdLinFitTask ()

virtual InitStatus	Init ()

virtual void	Exec (Option_t *opt)

virtual void	Finish ()

void	SetRadLen (double x)

Protected Member Functions
double	FCN_MS (const double *vars)

double	ScatteredAngle (double radLen)

double	line3DfitMS (TGraph2DErrors gr, const TVector3 &posSeed, const TVector3 &dirSeed, std::vector< double > &fitpar, TMatrixDSym covmatrix)

double	GetSigmaMS (int side)

	ClassDef (PndLmdLinFitTask, 2)

Static Protected Member Functions
static double	distance_MS (double x, double y, double z, double errx, double erry, const double p, double zpr)

Protected Attributes
TClonesArray *	fTCandArray

TClonesArray *	fRecoArray

TString	fTCandBranchName

TString	fRecoBranchName

TString	fOutputBranchName

TString	fOutputFolder

TClonesArray *	fTrackArray

double	fPbeam

double	fsigmaMSa

double	fsigmaMSb

bool	hitMergedfl [4]

double	fPDGCode

int	fCharge

PndGeoHandling *	fGeoH

ROOT::Math::Minimizer *	fmin

TGraph2DErrors *	fGraph2D

double	ftotRadLen

Detailed Description

Definition at line 30 of file PndLmdLinFitTask.h.

Constructor & Destructor Documentation

PndLmdLinFitTask::PndLmdLinFitTask ( )

Definition at line 44 of file PndLmdLinFitTask.cxx.

References fCharge, fmin, fOutputBranchName, fOutputFolder, fPbeam, fPDGCode, fRecoBranchName, fsigmaMSa, fsigmaMSb, fTCandBranchName, ftotRadLen, and hitMergedfl.

     : FairTask("3D-Straight-Line-Fit"),
       fTCandArray(0),
       fRecoArray(0),
       fTrackArray(0),
       fGeoH(0) {
   fTCandBranchName = "LMDTrackCand";
   // fRecoBranchName = "LMDHitsStrip";
   fRecoBranchName = "LMDHitsPixel";
   //  fRecoBranchName = "LmdHits";
   fOutputBranchName = "LMDPndTrack";
   fOutputFolder = "PndLmd";
 
   fPbeam = 0;
   fPDGCode = -2212;  // barp
   fCharge = -1;      // barp
   fsigmaMSa = 0;
   fsigmaMSb = 0;
   ftotRadLen = 0;
 
   fmin = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
 
   fmin->SetMaxFunctionCalls(1000000);
   fmin->SetMaxIterations(100000);
   fmin->SetTolerance(0.001);
 
   for (int ih = 0; ih < 4; ih++) hitMergedfl[ih] = false;
 }

PndLmdLinFitTask::PndLmdLinFitTask	(	TString	tTCandBranchName,
		TString	tRecoBranchName,
		TString	tOutputBranchName = `"LMDPndTrack"`,
		TString	tOutputFolder = `"PndLmd"`
	)

Definition at line 72 of file PndLmdLinFitTask.cxx.

References fCharge, fmin, fOutputBranchName, fOutputFolder, fPbeam, fPDGCode, fRecoBranchName, fsigmaMSa, fsigmaMSb, fTCandBranchName, ftotRadLen, and hitMergedfl.

     : FairTask("3D-Straight-Line-Fit"),
       fTCandArray(0),
       fRecoArray(0),
       fTrackArray(0),
       fGeoH(0) {
   fTCandBranchName = tTCandBranchName;
   fRecoBranchName = tRecoBranchName;
   fOutputBranchName = tOutputBranchName;
   fOutputFolder = tOutputFolder;
 
   fPbeam = 0;
   fPDGCode = -2212;  // barp
   fCharge = -1;      // barp
   //  fsigmaMS = 0;
   fsigmaMSa = 0;
   fsigmaMSb = 0;
   ftotRadLen = 0;
 
   fmin = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
 
   fmin->SetMaxFunctionCalls(1000000);
   fmin->SetMaxIterations(100000);
   fmin->SetTolerance(0.001);
 
   for (int ih = 0; ih < 4; ih++) hitMergedfl[ih] = false;
 }

PndLmdLinFitTask::~PndLmdLinFitTask ( )

virtual

Definition at line 103 of file PndLmdLinFitTask.cxx.

                                     {
   //    delete fmin;
   cout << "PndLmdLinFitTask::~PndLmdLinFitTask()" << endl;
 }

Member Function Documentation

PndLmdLinFitTask::ClassDef	(	PndLmdLinFitTask	,
		2
	)

protected

double PndLmdLinFitTask::distance_MS	(	double	x,
		double	y,
		double	z,
		double	errx,
		double	erry,
		const double *	p,
		double *	zpr
	)

staticprotected

Definition at line 318 of file PndLmdLinFitTask.cxx.

References Double_t, and sqrt().

Referenced by FCN_MS().

                                                   {
   double THfunc[8] = {1, 1, 1, 1, 1, 1, 1, 1};
   //  p[5] = sqrt(1-p[1]*p[1]-p[3]*p[3]);
   // Double_t t_min = p[1]*(x-p[0])+p[3]*(y-p[2])+p[5]*(z-p[4]);
   //  Double_t t_min = p[5]*(z-p[4]);
   // Double_t t_min = p[1]*(x-p[0])+p[3]*(y-p[2]);
   Double_t t_min = p[1] * (x - p[0]) + p[3] * (y - p[2]) +
                    sqrt(1 - p[1] * p[1] - p[3] * p[3]) * (z - p[4]);
   //  Double_t t_min = (z-p[4]);//TEST
   for (int iz = 0; iz < 8; iz++)
     if (((p[4] + t_min) - zpr[iz]) <= 0) THfunc[iz] = 0;
 
   double funcX = p[0] + p[1] * t_min +
                  (-p[6] + p[7]) * ((p[4] + t_min) - zpr[0]) * THfunc[0] +
                  (-p[8] + p[9]) * ((p[4] + t_min) - zpr[2]) * THfunc[2] +
                  (-p[10] + p[11]) * ((p[4] + t_min) - zpr[4]) * THfunc[4] +
                  (-p[12] + p[13]) * ((p[4] + t_min) - zpr[6]) * THfunc[6];
   double funcY = p[2] + p[3] * t_min +
                  (-p[14] + p[15]) * ((p[4] + t_min) - zpr[0]) * THfunc[0] +
                  (-p[16] + p[17]) * ((p[4] + t_min) - zpr[2]) * THfunc[2] +
                  (-p[18] + p[19]) * ((p[4] + t_min) - zpr[4]) * THfunc[4] +
                  (-p[20] + p[21]) * ((p[4] + t_min) - zpr[6]) * THfunc[6];
 
   double fdx = TMath::Power((x - funcX) / errx, 2);
   double fdy = TMath::Power((y - funcY) / erry, 2);
   double fdz = 0;
   double fchi2 = fdx + fdy + fdz;
   return fchi2;
 }

void PndLmdLinFitTask::Exec ( Option_t * opt )

virtual

Definition at line 163 of file PndLmdLinFitTask.cxx.

References Double_t, fCharge, fGeoH, fPbeam, fPDGCode, fRecoArray, fTCandArray, fTrackArray, fVerbose, PndTrackCandHit::GetDetId(), PndTrackCandHit::GetHitId(), PndSdsMergedHit::GetIsMerged(), PndTrackCand::GetNHits(), PndGeoHandling::GetOUVShortId(), PndSdsHit::GetPosition(), PndSdsHit::GetSensorID(), PndTrackCand::GetSortedHit(), hitMergedfl, line3DfitMS(), pid(), and track.

                                      {
   // Reset output Array
   if (fTrackArray == 0) Fatal("PndLmdLinFitTask::Exec", "No TrackArray");
   fTrackArray->Delete();
 
   Int_t ntcand = fTCandArray->GetEntriesFast();
 
   // Detailed output
   if (fVerbose > 1)
     std::cout << " -I- PndLmdLinFitTask: contains " << ntcand
               << " RhoCandidates" << std::endl;
   if (fVerbose > 2) {
     std::cout << " Detailed Debug info on the candidates:" << std::endl;
     unsigned int index = 12345;  // detid=12345,  //[R.K.03/2017] unused
                                  // variable
     for (Int_t itr = 0; itr < ntcand; ++itr) {
       PndTrackCand* trcnd = (PndTrackCand*)fTCandArray->At(itr);
       std::cout << "TrackCand no. " << itr << " has " << trcnd->GetNHits()
                 << " hits." << std::endl;
       std::cout << "Point: \t Index: " << std::endl;
       for (unsigned int ihit = 0; ihit < trcnd->GetNHits();
            ihit++) {                                         // fill Graph
         PndTrackCandHit theHit = trcnd->GetSortedHit(ihit);  // get hit
         index = theHit.GetHitId();
         // detid = theHit.GetDetId(); //[R.K.03/2017] unused variable
         std::cout << ihit << "\t" << index << std::endl;
       }
     }
   }
 
   // Fitting
   // ----------------------------------------------------------------------------------
   //  if(fVerbose>1)
   // std::cout<<" -I- PndLmdLinFitTask: start Fitting "<<std::endl;
   int rec_tkr = 0;
   for (Int_t track = 0; track < ntcand; track++) {
     PndTrackCand* trcnd = (PndTrackCand*)fTCandArray->At(track);
     const int numPts = trcnd->GetNHits();  // read how many points in this track
 
     TGraph2DErrors fitme(numPts);  // new graph for fitting
     // Int_t firstHit=-1, lastHit=-1; //[R.K.03/2017] unused variable
     TVector3 hit0, hit1;
     for (int ihit = 0; ihit < numPts; ihit++) {            // fill Graph
       PndTrackCandHit theHit = trcnd->GetSortedHit(ihit);  // get hit
       Int_t index = theHit.GetHitId();
       // Int_t detId = theHit.GetDetId(); //[R.K. 01/2017] unused variable
       // if(fVerbose>2) std::cout << "Point: "<< ihit<< " index: "<< index
       // <<std::endl;
 
       //      PndSdsHit* addHit = (PndSdsHit*) fRecoArray->At(index);
       PndSdsMergedHit* addHit = (PndSdsMergedHit*)fRecoArray->At(index);
       //      cout<<"IsMerged??? "<<addHit->GetIsMerged()<<endl;
       hitMergedfl[ihit] = addHit->GetIsMerged();
 
       TVector3 addPos = addHit->GetPosition();
       double xhit = addPos.X();
       double yhit = addPos.Y();
       double zhit = addPos.Z();
       //      lmddim->Transform_global_to_lmd_local(xhit,yhit,zhit,false);
       // TVector3 addPos2(xhit,yhit,zhit);
       if (ihit == 0) {
         // firstHit=index; //[R.K.03/2017] unused variable
         hit0 = addPos;
       } else {
         if (ihit == 2) hit1 = addPos;
         //      if(ihit==1)       hit1 = addPos;
         //      if(ihit==numPts-1)        hit1 = addPos;
         if (ihit == numPts - 1) {
           // lastHit=index; //[R.K.03/2017] unused variable
           //      hit1 = addPos2;
         }
       }
       double errxhit = addHit->GetDx();
       double erryhit = addHit->GetDy();
       double errzhit = addHit->GetDz();
       fitme.SetPoint(ihit, xhit, yhit, zhit);
       fitme.SetPointError(ihit, errxhit, erryhit, errzhit);
       // fitme.SetPointError(ihit, 0.5*errxhit,0.5*erryhit,0.5*errzhit);
     }  // end of Hits in TCand
 
     TVector3 dirSeed = hit1 - hit0;
     TVector3 posSeed = hit0;
     // cout<<"DirSeed before fit and norm:"<<endl;
     // dirSeed.Print();
     dirSeed = dirSeed.Unit();
     std::vector<double> fit_parameters(22);  // fit-parameter
     TMatrixDSym* COVmatrix = new TMatrixDSym(6);
     Double_t accuracy = line3DfitMS(&fitme, posSeed, dirSeed, fit_parameters,
                                     COVmatrix);  // with kink angles
     //    Double_t accuracy = line3Dfit(numPts, &fitme, posSeed, dirSeed,
     //    parFit, COVmatrix); //MS taking into account by hits errors
 
     // save as PndTrack
     TVector3 FitPoint(fit_parameters[0], fit_parameters[2], fit_parameters[4]);
     TVector3 FitDir(fit_parameters[1], fit_parameters[3], fit_parameters[5]);
     //  TVector3 FitDir(parFit[1]+parFit[9], parFit[3]+parFit[17],
     //  sqrt(1-(parFit[1]+parFit[9])*(parFit[1]+parFit[9])-(parFit[3]+parFit[17])*(parFit[3]+parFit[17])));//TEST
     TVector3 FitMom = FitDir * fPbeam;
     Double_t COVmatrixPosMom[6][6];
 
     for (int ij = 0; ij < 6; ij++) {
       for (int ji = 0; ji < 6; ji++) COVmatrixPosMom[ij][ji] = 0;
     }
 
     int iconver[6] = {3, 0, 4, 1, 5, 2};
     for (int ij = 0; ij < 6; ij++) {
       for (int ji = 0; ji < 6; ji++) {
         if (ij == 1 || ij == 3 || ij == 5) (*COVmatrix)(ij, ji) *= fPbeam;
         if (ji == 1 || ji == 3 || ji == 5) (*COVmatrix)(ij, ji) *= fPbeam;
         int km = iconver[ij];
         int mk = iconver[ji];
         COVmatrixPosMom[km][mk] = (*COVmatrix)(ij, ji);
       }
     }
     delete COVmatrix;
 
     // Read info about 1st plane(sensor)
     PndTrackCandHit theHit = trcnd->GetSortedHit(0);  // get 1st hit
     Int_t hitID = theHit.GetHitId();
     int sysID = theHit.GetDetId();
     PndSdsHit* myHit = (PndSdsHit*)(fRecoArray->At(hitID));
     Int_t id = myHit->GetSensorID();
 
     TVector3 oo, uu, vv;
     fGeoH->GetOUVShortId(id, oo, uu, vv);
     // TVector3 o = oo;
     TVector3 o = FitPoint;
     TVector3 dj = uu;
     TVector3 dk = vv;
 
     FairTrackParP* trkfitp = new FairTrackParP(
         FitPoint, FitMom, COVmatrixPosMom, fCharge, o, dj, dk);
     int flagpndtrk = 0;
     if (accuracy > 1e3)
       flagpndtrk = -1;  // quality of track: bad if chi2 is too large
     int ndf = 2 * (trcnd->GetNHits()) - 4;  // number d.o.f
     int pid = fPDGCode;
 
     PndTrack* trackfit = new PndTrack(
         *trkfitp, *trkfitp, *trcnd, flagpndtrk, accuracy, ndf, pid, track,
         sysID);  // TODO: FairTrackParP at 1st and last point???
     new ((*fTrackArray)[rec_tkr]) PndTrack(*(trackfit));  // save Track
 
     delete trackfit;  // TEST
     rec_tkr++;
     //}
   }  // end of TCand's
 
   // Done--------------------------------------------------------------------------------------
   if (fVerbose > 2) std::cout << "Fitting done" << std::endl;
 
   return;
 }

double PndLmdLinFitTask::FCN_MS ( const double * vars )

protected

Definition at line 351 of file PndLmdLinFitTask.cxx.

References distance_MS(), fGraph2D, GetSigmaMS(), i, npoints, sqrt(), x, y, and z.

Referenced by line3DfitMS().

                                                   {
   double* x = fGraph2D->GetX();
   double* y = fGraph2D->GetY();
   double* z = fGraph2D->GetZ();
   double* errx = fGraph2D->GetEX();
   double* erry = fGraph2D->GetEY();
   int npoints = fGraph2D->GetN();
   double sum(0.0);
   double erralA = GetSigmaMS(0);
   double erralB = GetSigmaMS(1);
   double zpr[8];
   int izcur = 0;
   for (int izpl = 0; izpl < 4; izpl++) {
     zpr[izcur] = z[izpl];
     izcur++;
     zpr[izcur] = z[izpl];
     izcur++;
   }
   for (int i = 0; i < npoints; ++i) {
     double chi2 = distance_MS(x[i], y[i], z[i], errx[i], erry[i], vars, zpr);
     sum += chi2;
   }
   double fdal = 0;
   for (int jms = 0; jms < 8; jms = jms + 2) {
     double errTot = sqrt(erralA * erralA + erralB * erralB);
     fdal += TMath::Power(((vars[7 + jms] - vars[6 + jms]) / errTot), 2);
     fdal += TMath::Power(((vars[15 + jms] - vars[14 + jms]) / errTot), 2);
   }
   sum += fdal;
   return sum;
 }

void PndLmdLinFitTask::Finish ( )

virtual

Definition at line 516 of file PndLmdLinFitTask.cxx.

516 {}

double PndLmdLinFitTask::GetSigmaMS ( int side )

inlineprotected

Definition at line 56 of file PndLmdLinFitTask.h.

References fsigmaMSa, and fsigmaMSb.

Referenced by FCN_MS().

                               {
     if (side < 1)
       return fsigmaMSa;
     else
       return fsigmaMSb;
   }

InitStatus PndLmdLinFitTask::Init ( )

virtual

Definition at line 108 of file PndLmdLinFitTask.cxx.

References fCharge, fGeoH, fmin, fOutputBranchName, fOutputFolder, fPbeam, fPDGCode, fRecoArray, fRecoBranchName, fRun, fsigmaMSa, fsigmaMSb, fTCandArray, fTCandBranchName, ftotRadLen, fTrackArray, fVerbose, PndGeoHandling::Instance(), par, rtdb, and ScatteredAngle().

Referenced by runLumi4Fitter(), and runLumiPixel4Fitter().

                                   {
   // Get ROOT Manager
   FairRootManager* ioman = FairRootManager::Instance();
 
   if (ioman == 0) {
     Error("PndLmdLinFitTask::Init", "RootManager not instantiated!");
     return kERROR;
   }
 
   // Get input collection
   fTCandArray = (TClonesArray*)ioman->GetObject(fTCandBranchName);
 
   if (fTCandArray == 0) {
     Error("PndLmdLinFitTask::Init", "trackcand-array not found!");
     return kERROR;
   }
 
   fRecoArray = (TClonesArray*)ioman->GetObject(fRecoBranchName);
 
   if (fRecoArray == 0) {
     Error("PndLmdLinFitTask::Init", "reco-array not found!");
     return kERROR;
   }
 
   fTrackArray = new TClonesArray("PndTrack");
   ioman->Register(fOutputBranchName, fOutputFolder, fTrackArray, kTRUE);
 
   // read beam momentum from base
   FairRun* fRun = FairRun::Instance();
   FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
   FairBaseParSet* par =
       (FairBaseParSet*)(rtdb->findContainer("FairBaseParSet"));
   fPbeam = par->GetBeamMom();
   fPDGCode = -2212;  // barp
   fCharge = -1;      // barp
 
   fGeoH = PndGeoHandling::Instance();
   //  double totRadLen = 0.00306;//rad.length of whole plane
   cout << "RadLeng = " << ftotRadLen << endl;
   double totRadLen = ftotRadLen;
   if (totRadLen < 1e-6)
     totRadLen = 0.00306;  // standart rad.length of whole plane
   // totRadLen -=2*0.000175;// -flex-cable
   // totRadLen -=0.00053;// -sensor
   fsigmaMSb = ScatteredAngle(totRadLen);
   double totRadLenCable = 0.053;  // sensor
   totRadLenCable += 0.05;         // rest
   fsigmaMSa = ScatteredAngle(totRadLenCable);
 
   fmin->SetPrintLevel(fVerbose);
 
   std::cout << "-I- PndLmdLinFitTask: Initialisation successfull" << std::endl;
   return kSUCCESS;
 }

double PndLmdLinFitTask::line3DfitMS	(	TGraph2DErrors *	gr,
		const TVector3 &	posSeed,
		const TVector3 &	dirSeed,
		std::vector< double > &	fitpar,
		TMatrixDSym *	covmatrix
	)

protected

!!!!!!!!!!!!!!!!! z0, dz weren't fitted

!!!!!!!!!!!!!!!!!

Definition at line 398 of file PndLmdLinFitTask.cxx.

References Double_t, f, fabs(), FCN_MS(), fGraph2D, fmin, fsigmaMSa, fsigmaMSb, fVerbose, i, and sqrt().

Referenced by Exec().

                                                              {
   const int nparams = fitpar.size();
   ROOT::Math::Functor f(this, &PndLmdLinFitTask::FCN_MS, nparams);
   fmin->SetFunction(f);
 
   if (fVerbose > 2)
     cout << "PndLmdLinFitTask::line3Dfit with SEED is used (multiple "
             "scattering taking into account with kinks)"
          << endl;
   Int_t Npoint = gr->GetN();
   Double_t ErrX1 = gr->GetErrorX(0);
   Double_t ErrY1 = gr->GetErrorY(0);
   Double_t ErrZ1 = gr->GetErrorZ(0);
 
   Double_t ErrX2 = gr->GetErrorX(2);
   Double_t ErrY2 = gr->GetErrorY(2);
 
   Double_t errRx = TMath::Hypot(ErrX1, ErrX2);
   Double_t errRy = TMath::Hypot(ErrY1, ErrY2);
 
   fGraph2D = gr;
 
   if (fVerbose > 5) {
     cout << "Number of hits = " << Npoint << endl;
     cout << "posSeed:" << endl;
     posSeed.Print();
     TVector3 ErrPosSeed(ErrX1, ErrY1, ErrZ1);
     cout << "ErrposSeed:" << endl;
     ErrPosSeed.Print();
     cout << "dirSeed:" << endl;
     dirSeed.Print();
   }
 
   fmin->SetVariable(0, "x0", posSeed.X(), ErrX1);
   fmin->SetVariable(1, "Ax", dirSeed.X(), errRx);
   fmin->SetVariable(2, "y0", posSeed.Y(), ErrY1);
   fmin->SetVariable(3, "Ay", dirSeed.Y(), errRy);
 
   fmin->SetFixedVariable(4, "z0", posSeed.Z());
   fmin->SetFixedVariable(5, "Az", dirSeed.Z());
 
   double scaling_factor(1e-4);
   fmin->SetVariable(6, "al0x_a", 0.0, scaling_factor * fsigmaMSa);
   fmin->SetVariable(7, "al0x_b", 0.0, scaling_factor * fsigmaMSb);
   fmin->SetFixedVariable(8, "al1x_a", 0.0);
   fmin->SetVariable(9, "al1x_b", 0.0, scaling_factor * fsigmaMSb);
   fmin->SetFixedVariable(10, "al2x_a", 0.0);
   fmin->SetVariable(11, "al2x_b", 0.0, scaling_factor * fsigmaMSb);
   fmin->SetFixedVariable(12, "al3x_a", 0.0);
   fmin->SetFixedVariable(13, "al3x_b", 0.0);
 
   fmin->SetVariable(14, "al0y_a", 0.0, scaling_factor * fsigmaMSa);
   fmin->SetVariable(15, "al0y_b", 0.0, scaling_factor * fsigmaMSb);
   fmin->SetFixedVariable(16, "al1y_a", 0.0);
   fmin->SetVariable(17, "al1y_b", 0.0, scaling_factor * fsigmaMSb);
   fmin->SetFixedVariable(18, "al2y_a", 0.0);
   fmin->SetVariable(19, "al2y_b", 0.0, scaling_factor * fsigmaMSb);
   fmin->SetFixedVariable(20, "al3y_a", 0.0);
   fmin->SetFixedVariable(21, "al3y_b", 0.0);
 
   if (Npoint < 4) {
     fmin->FixVariable(11);
     fmin->FixVariable(19);
   }
 
   fmin->Minimize();
 
   double recpres = 1e-7;
 
   if (fmin->Edm() > 1e2 * recpres) {
     return 1e6;
   }
 
   // get fit parameters
   for (int i = 0; i < nparams; ++i) {
     fitpar[i] = fmin->X()[i];
   }
 
   for (size_t i = 0; i < 6; i++) {
     for (size_t j = 0; j < 6; j++) {
       (*covmatrix)(i, j) = fmin->CovMatrix(i, j);
     }
   }
 
   (*covmatrix)(4, 4) = (gr->GetErrorZ(0) * gr->GetErrorZ(0)) / 12.;
 
   fitpar[5] = sqrt(1 - fitpar[1] * fitpar[1] - fitpar[3] * fitpar[3]);
   double dp5_dp1 = fitpar[1] / fitpar[5];
   double dp5_dp3 = fitpar[3] / fitpar[5];
   double errdz2 = pow(dp5_dp1, 2) * (*covmatrix)(1, 1) +
                   pow(dp5_dp3, 2) * (*covmatrix)(3, 3) +
                   2 * fabs(dp5_dp1 * dp5_dp3 * (*covmatrix)(1, 3));
   (*covmatrix)(5, 5) = errdz2;
 
   Double_t chi2 = fmin->MinValue() / (2. * Npoint - 4);
 
   if (fVerbose > 2) {
     cout << " chi2 = " << chi2 << endl;
     cout << "[AFTER FIT] start.point: (" << fitpar[0] << ", " << fitpar[2]
          << ", " << fitpar[4] << ")" << endl;
     cout << "[AFTER FIT] dir: (" << fitpar[1] << ", " << fitpar[3] << ", "
          << fitpar[5] << ")"
          << " sqrt(p1^2+p3^2+p5^2) = "
          << sqrt(fitpar[1] * fitpar[1] + fitpar[3] * fitpar[3] +
                  fitpar[5] * fitpar[5])
          << endl;
   }
 
   return chi2;
 }

double PndLmdLinFitTask::ScatteredAngle ( double radLen )

protected

Definition at line 383 of file PndLmdLinFitTask.cxx.

References Double_t, fPbeam, CAMath::Log(), and CAMath::Sqrt().

Referenced by Init().

                                                      {
   // Calculation of ThetaMS -------------------------------------
   // Charge & mass of particle
   // Int_t PDGCode = -2212; //barp //[R.K. 01/2017] unused variable
   Double_t fMass = 0.938272046;
   Double_t Ebeam = TMath::Hypot(fPbeam, fMass);
   TLorentzVector LorMom(0, 0, fPbeam, Ebeam);
   Double_t beta = LorMom.Beta();
   Double_t X_to_X0 = radLen;
   Double_t thetaMS = 13.6 * 1e-3 * TMath::Sqrt(X_to_X0) *
                      (1 + 0.038 * TMath::Log(X_to_X0)) / (beta * fPbeam);
   return thetaMS;
 }