Millepede Class Reference

#include <AlignmentTools/Millepede.h>

Public Member Functions
	Millepede ()
	Standard constructor.
	~Millepede ()
	Destructor.
bool	initialize ()
	Initialization.
bool	InitMille (bool DOF[], double Sigm[], int nlay, int nglo_on_lay, int nloc, double startfact, int nstd, double res_cut, double res_cut_init)
bool	MakeGlobalFit (double par[], double error[], double pull[])
bool	ParGlo (int index, double param)
bool	ParSig (int index, double sigma)
bool	ConstF (double dercs[], double rhs)
bool	EquLoc (double dergb[], double derlc[], double dernl[], double rmeas, double sigma)
bool	ZerLoc (double dergb[], double derlc[], double dernl[])
bool	FitLoc (int n, double track_params[], int single_fit)
int	GetTrackNumber ()
void	SetTrackNumber (int value)
	Millepede ()
	Standard constructor.
	~Millepede ()
	Destructor.
bool	initialize ()
	Initialization.
bool	InitMille (bool DOF[], double Sigm[], int nglo, int nloc, double startfact, int nstd, double res_cut, double res_cut_init)
bool	MakeGlobalFit (double par[], double error[], double pull[])
bool	ParGlo (int index, double param)
bool	ParSig (int index, double sigma)
bool	ConstF (double dercs[], double rhs)
bool	EquLoc (double dergb[], double derlc[], double dernl[], double rmeas, double sigma)
bool	ZerLoc (double dergb[], double derlc[], double dernl[])
bool	FitLoc (int n, double track_params[], int single_fit)
int	GetTrackNumber ()
void	SetTrackNumber (int value)

Detailed Description

Author

Sebastien Viret

Date

2005-07-29

Definition at line 19 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/CgemAlignAlg/Millepede.h.

Constructor & Destructor Documentation

Millepede() [1/2]

Millepede::Millepede

(

)

Standard constructor.

Definition at line 25 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{}

~Millepede() [1/2]

Millepede::~Millepede

(

)

Destructor.

Definition at line 30 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{} 

Millepede() [2/2]

Millepede::Millepede

(

)

Standard constructor.

~Millepede() [2/2]

Millepede::~Millepede

(

)

Destructor.

Member Function Documentation

ConstF() [1/2]

bool Millepede::ConstF	(	double	dercs[],
		double	rhs )

Definition at line 275 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{  
     if (ncs>=mcs) // mcs is defined in Millepede.h
     {
      cout << "Too many constraints !!!" << endl;
      return false;
     }
 
     for (int i=0; i<nagb; i++) {adercs[ncs][i] = dercs[i];}
 
     arhs[ncs] = rhs;
     ncs++ ;
     cout << "Number of constraints increased to " << ncs << endl;
     return true;
}

ConstF() [2/2]

bool Millepede::ConstF	(	double	dercs[],
		double	rhs )

EquLoc() [1/2]

bool Millepede::EquLoc	(	double	dergb[],
		double	derlc[],
		double	dernl[],
		double	rmeas,
		double	sigma )

Definition at line 305 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{   
     //if (verbose_mode) cout << " sigma " << sigma << endl;
 
     if (sigma<=0.0) // If parameter is fixed, then no equation
     {
      for (int i=0; i<nalc; i++)
      {
           derlc[i] = 0.0;
      }
      for (int i=0; i<nagb; i++)
      {
           dergb[i] = 0.0;
      }
      return true;
     }
 
     // Serious equation, initialize parameters
 
     double wght =  1.0/(sigma*sigma);
     int nonzer  =  0;
     int ialc    = -1;
     int iblc    = -1;
     int iagb    = -1;
     int ibgb    = -1;
     
     if (verbose_mode) cout<<"meas "<<rmeas<<endl;
     for (int i=0; i<nalc; i++) // Retrieve local param interesting indices
     {
          if (verbose_mode) cout<<"i = "<<i<<", derlc = "<<derlc[i]<<endl;
      if (derlc[i]!=0.0)
      {
           nonzer++;
           if (ialc == -1) ialc=i;  // first index
           iblc = i;                // last index
      }
     }
 
     if (verbose_mode) cout << ialc << " / " << iblc << endl;
 
     if (verbose_mode) cout<<"weight "<<wght<<endl;
     for (int i=0; i<nagb; i++)  // Idem for global parameters
     {
          if (verbose_mode) cout<<"i = "<<i<<", dergb = "<<dergb[i]<<endl;
      if (dergb[i]!=0.0)
      {
           nonzer++;
           if (iagb == -1) iagb=i;  // first index
           ibgb = i;            // last index
      }
     }
 
     if (verbose_mode) cout << iagb << " / " << ibgb << endl;
 
     indst.push_back(-1);
     arest.push_back(rmeas);
     arenl.push_back(0.);
 
     for (int i=ialc; i<=iblc; i++)
     {
      if (derlc[i]!=0.0)
      {
           indst.push_back(i);
           arest.push_back(derlc[i]);
           arenl.push_back(0.0);
           derlc[i]   = 0.0;
      }
     }
 
     indst.push_back(-1);
     arest.push_back(wght);
     arenl.push_back(0.0);
 
     for (int i=iagb; i<=ibgb; i++)
     {
      if (dergb[i]!=0.0)
      {
           indst.push_back(i);
           arest.push_back(dergb[i]);
           arenl.push_back(dernl[i]);
           dergb[i]   = 0.0;
      }
     }  
 
     if (verbose_mode) cout << "Out Equloc --  NST = " << arest.size() << endl;
 
     return true;   
}

EquLoc() [2/2]

bool Millepede::EquLoc	(	double	dergb[],
		double	derlc[],
		double	dernl[],
		double	rmeas,
		double	sigma )

FitLoc() [1/2]

bool Millepede::FitLoc	(	int	n,
		double	track_params[],
		int	single_fit )

Definition at line 435 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{
     // Few initializations
 
     int i, j, k, ik, ij, ist, nderlc, ndergl, ndf;
     int ja      = -1;
     int jb      = 0;
     int nagbn   = 0;
 
     double rmeas, wght, rms, cutval;
 
     double summ  = 0.0;
     int    nsum  = 0;
     nst   = 0; 
     nst   = arest.size();
     if(debug_mode) cout<<"Start FitLoc, nst = "<<nst<<endl;
 
 
     // Fill the track store at first pass
 
     if (itert < 2 && single_fit != 1)  // Do it only once 
     {
      if (debug_mode) cout << "Store equation no: " << n << endl; 
 
      for (i=0; i<nst; i++)    // Store the track parameters
      {
           storeind.push_back(indst[i]);
           storeare.push_back(arest[i]);
           storenl.push_back(arenl[i]);
 
           if (arenl[i] != 0.) arest[i] = 0.0; // Reset global derivatives if non linear and first iteration
      }
 
      arenl.clear();
 
      storeplace.push_back(storeind.size());
 
      if (verbose_mode) cout << "StorePlace size = " << storeplace[n] << endl; 
      if (verbose_mode) cout << "StoreInd size   = " << storeind.size() << endl; 
     }  
 
 
     for (i=0; i<nalc; i++) // reset local params
     {
      blvec[i] = 0.0;
 
      for (j=0; j<nalc; j++) {clmat[i][j] = 0.0;}
     }
 
     for (i=0; i<nagb; i++) {indnz[i] = -1;} // reset mixed params
 
 
     /*
 
     LOOPS : HOW DOES IT WORKS ?    
 
     Now start by reading the informations stored with EquLoc.
     Those informations are in vector INDST and AREST.
     Each -1 in INDST delimits the equation parameters:
 
     First -1  ---> rmeas in AREST 
     Then we have indices of local eq in INDST, and derivatives in AREST
     Second -1 ---> weight in AREST
     Then follows indices and derivatives of global eq.
     ....
 
     We took them and store them into matrices.
 
     As we want ONLY local params, we substract the part of the estimated value
     due to global params. Indeed we could have already an idea of these params,
     with previous alignment constants for example (set with PARGLO). Also if there
     are more than one iteration (FITLOC could be called by FITGLO)
 
     */
 
 
     //
     // FIRST LOOP : local track fit
     //
 
     ist = 0;
     indst.push_back(-1);
 
     while (ist <= nst)
     {
      if (indst[ist] == -1)
      {
           if (ja == -1)     {ja = ist;}  // First  0 : rmeas
           else if (jb == 0) {jb = ist;}  // Second 0 : weight 
           else                           // Third  0 : end of equation  
           {
            rmeas   = arest[ja];
            wght    = arest[jb];
            if (verbose_mode) cout << "rmeas = " << rmeas << endl ;
            if (verbose_mode) cout << "wght = " << wght << endl ;
 
            for (i=(jb+1); i<ist; i++)   // Now suppress the global part   
             // (only relevant with iterations)
            {
             j = indst[i];              // Global param indice
             if (verbose_mode) cout << "dparm[" << j << "] = " << dparm[j] << endl;        
             if (verbose_mode) cout << "Starting misalignment = " << pparm[j] << endl;        
             rmeas -= arest[i]*(pparm[j]+dparm[j]);
            }
 
            if (verbose_mode) cout << "rmeas after global stuff removal = " << rmeas << endl ;
 
            for (i=(ja+1); i<jb; i++)    // Finally fill local matrix and vector
            {
             j = indst[i];   // Local param indice (the matrix line) 
             // Beta   1/err2 residual derivative  -->   Beta_j = d R^2/ d p_j * R/(sigma)^2
             blvec[j] += wght*rmeas*arest[i];  // See note for precisions 
 
             if (verbose_mode) cout << "blvec[" << j << "] = " << blvec[j] << endl ;
 
             for (k=(ja+1); k<=i ; k++) // Symmetric matrix, don't bother k>j coeffs
             {
                  ik = indst[k];                        
                  // Gamma_jk     1/err2 derivative    --> Gamma_jk = dR^2/dp_j  * dR^2/dp_k * 1/(sigma)^2 
                  clmat[j][ik] += wght*arest[i]*arest[k];
 
                  if (verbose_mode) cout << "clmat[" << j << "][" << ik << "] = " << clmat[j][ik] << endl;
             } 
            }  
            ja = -1;
            jb = 0;
            ist--;
           } // End of "end of equation" operations
      } // End of loop on equation
      ist++;
     } // End of loop on all equations used in the fit
 
 
     //
     // Local params matrix is completed, now invert to solve...
     //
 
     nrank = 0;  // Rank is the number of nonzero diagonal elements, what are the last two variables? 
     nrank = Millepede::SpmInv(clmat, blvec, nalc, scdiag, scflag);
 
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << " __________________________________________________" << endl;
     if (debug_mode) cout << " Printout of local fit  (FITLOC)  with rank= "<< nrank << endl;
     if (debug_mode) cout << " Result of local fit :      (index/parameter/error)" << endl;
 
     for (i=0; i<nalc; i++)
     {
      if (debug_mode) cout << std::setprecision(4) << std::fixed;
      if (debug_mode) cout << std::setw(20) << i << "   /   " << std::setw(10) 
                   << blvec[i] << "   /   " << sqrt(clmat[i][i]) << endl;// Why the error is sqrt(clmat[i][i]) ? By guoaq-19-09-26  
     }
 
 
     // Store the track params and errors
 
     for (i=0; i<nalc; i++)
     {
      track_params[2*i] = blvec[i];
      track_params[2*i+1] = sqrt(fabs(clmat[i][i]));
     }
 
 
     //
     // SECOND LOOP : residual calculation
     //
 
     ist = 0;
     ja = -1;
     jb = 0;
 
     while (ist <= nst)
     {
      if (verbose_mode)    cout<<"ist , nst = "<<ist<<", "<<nst<<endl;
      if (indst[ist] == -1)
      {
           if (ja == -1)     {ja = ist;}  // First  0 : rmeas
           else if (jb == 0) {jb = ist;}  // Second 0 : weight 
           else                           // Third  0 : end of equation  
           {    
            rmeas   = arest[ja];
            wght    = arest[jb];
 
            nderlc = jb-ja-1;    // Number of local derivatives involved
            ndergl = ist-jb-1;   // Number of global derivatives involved
 
            // Print all (for debugging purposes)
 
            if (verbose_mode) cout << "" << endl;
            if (verbose_mode) cout << std::setprecision(4) << std::fixed;
            if (verbose_mode) cout << ". equation:  measured value " << std::setw(13) 
                       << rmeas << " +/- " << std::setw(13) << 1.0/sqrt(wght) << endl;
            if (verbose_mode) cout << "Number of derivatives (global, local): " 
                       << ndergl << ", " << nderlc << endl;
            if (verbose_mode) cout << "Global derivatives are: (index/derivative/parvalue) " << endl;
 
            for (i=(jb+1); i<ist; i++)
            {if (verbose_mode) cout << indst[i] << " / " << arest[i] 
                        << " / " << pparm[indst[i]] << endl;} 
 
                    if (verbose_mode) cout << "Local derivatives are: (index/derivative) " << endl;
 
                    for (i=(ja+1); i<jb; i++) {if (verbose_mode) cout << indst[i] << " / " << arest[i] << endl;}      
 
                    // Now suppress local and global parts to RMEAS;
 
                    for (i=(ja+1); i<jb; i++) // First the local part 
                    {
             j = indst[i];
             rmeas -= arest[i]*blvec[j];
                    }
 
                    for (i=(jb+1); i<ist; i++) // Then the global part
                    {
             j = indst[i];
             rmeas -= arest[i]*(pparm[j]+dparm[j]);
                    }
 
                    // rmeas contains now the residual value
                    //  if (verbose_mode) cout << "Residual value : "<< rmeas << endl;
                    if (verbose_reject) cout << "Residual value : "<< rmeas << endl;
                    if (verbose_reject) cout << "Residual cut init value : "<< m_residual_cut_init << endl;
                    if (verbose_reject) cout << "Residual cut value : "<< m_residual_cut << endl;
 
                    // reject the track if rmeas is too important (outlier)
                    if (fabs(rmeas) >= m_residual_cut_init && itert <= 1)  
                    {
             //       if (verbose_mode) cout << "Rejected track !!!!!" << endl;
             if (verbose_reject) cout << "Rejected track !!!!!" << endl;
             locrej++;      
             indst.clear(); // reset stores and go to the next track 
             arest.clear();   
             return false;
                    }
 
                    if (fabs(rmeas) >= m_residual_cut && itert > 1)   
                    {
             //       if (verbose_mode) cout << "Rejected track !!!!!" << endl;
             if (verbose_reject) cout << "Rejected track !!!!!" << endl;
             locrej++;      
             indst.clear(); // reset stores and go to the next track 
             arest.clear();   
             return false;
                    }
                     
            if(debug_mode) cout<<"chisq , wght, rmeas = "<<summ<<", "<<wght<<", "<<rmeas<<endl; 
                    summ += wght*rmeas*rmeas ; // total chi^2
                    nsum++;                    // number of equations           
                    ja = -1;
                    jb = 0;
                    ist--;
           } // End of "end of equation" operations
      }   // End of loop on equation
      ist++;
     } // End of loop on all equations used in the fit
 
     ndf = nsum-nrank;  
     if (debug_mode) cout<<"nsum, nrank = "<<nsum<<", "<<nrank<<endl;
     rms = 0.0;
 
     if (debug_mode) cout << "Final chi square / degrees of freedom "<< summ << " / " << ndf << endl;
     //cout << "Final chi square / degrees of freedom "<< summ << " / " << ndf << endl;
 
 
 
 
     if (ndf > 0) rms = summ/float(ndf);  // Chi^2/dof
 
     if (single_fit == 0) loctot++;
 
     if (nstdev != 0 && ndf > 0 && single_fit != 1) // Chisquare cut
     {
      cutval = Millepede::chindl(nstdev, ndf)*cfactr;
          if (debug_mode) cout<<"nstdev, ndf, cfactr, cutval = "<<nstdev<<", "<<ndf<<", "<<cfactr<<", "<<cutval<<endl;
 
      if (debug_mode) cout << "Reject if Chisq/Ndf = " << rms << "  >  " << cutval << endl;
 
      if (rms > cutval) // Reject the track if too much...
      {
           if (verbose_mode) cout << "Rejected track !!!!!" << endl;
           if (single_fit == 0) locrej++;      
           indst.clear(); // reset stores and go to the next track 
           arest.clear();
           return false;
      }
     }
 
     if (single_fit == 1) // Stop here if just updating the track parameters
     {
      indst.clear(); // Reset store for the next track 
      arest.clear();
      return true;
     }
 
     //  
     // THIRD LOOP: local operations are finished, track is accepted 
     // We now update the global parameters (other matrices)
     //
     if(record_mill)
     {
         ofstream fmill("mill.csv",ios::app);
         //     itert        track        chisq           ndf         nHit      cut_chi_p_dnf     drho           phi0            dz               tanl           dx1             dy1             dz1             rz1
         fmill<<itert<<", "<< n << ", "<< summ << ", " << ndf <<", "<<nsum<<", "<<cutval<<", "<< blvec[0]<<", "<<blvec[1]<<", "<<blvec[3]<<", "<<blvec[4]<<", "<<dparm[0]<<", "<<dparm[3]<<", "<<dparm[6]<<", "<<dparm[9]<< endl;
         fmill.close();
     }
     
     ist = 0;
     ja = -1;
     jb = 0;
 
     while (ist <= nst)
     {
      if (indst[ist] == -1)
      {
           if (ja == -1)     {ja = ist;}  // First  0 : rmeas
           else if (jb == 0) {jb = ist;}  // Second 0 : weight 
           else                           // Third  0 : end of equation  
           {    
            rmeas   = arest[ja];
            wght    = arest[jb];
 
            for (i=(jb+1); i<ist; i++) // Now suppress the global part
            {
             j = indst[i];   // Global param indice
             rmeas -= arest[i]*(pparm[j]+dparm[j]);
            }
 
            // First of all, the global/global terms (exactly like local matrix)
 
            for (i=(jb+1); i<ist; i++)  
            {
             j = indst[i];   // Global param indice (the matrix line)          
 
             bgvec[j] += wght*rmeas*arest[i];  // See note for precisions
             if (verbose_mode) cout << "bgvec[" << j << "] = " << bgvec[j] << endl ;
 
             for (k=(jb+1); k<ist ; k++)
             {
                  ik = indst[k];                        
                  cgmat[j][ik] += wght*arest[i]*arest[k];
                  if (verbose_mode) cout << "cgmat[" << j << "][" << ik << "] = " << cgmat[j][ik] << endl;
             } 
            }
 
            // Now we have also rectangular matrices containing global/local terms.
 
            for (i=(jb+1); i<ist; i++) 
            {
             j  = indst[i];  // Global param indice (the matrix line)
             ik = indnz[j];  // Index of index
 
             if (ik == -1)    // New global variable
             {
                  for (k=0; k<nalc; k++) {clcmat[nagbn][k] = 0.0;} // Initialize the row
 
                  indnz[j] = nagbn;
                  indbk[nagbn] = j;
                  ik = nagbn;
                  nagbn++;
             }
 
             for (k=(ja+1); k<jb ; k++) // Now fill the rectangular matrix
             {
                  ij = indst[k];                        
                  clcmat[ik][ij] += wght*arest[i]*arest[k];
                  if (verbose_mode) cout << "clcmat[" << ik << "][" << ij << "] = " << clcmat[ik][ij] << endl;
             } 
            }
            ja = -1;
            jb = 0;
            ist--;
           } // End of "end of equation" operations
      }   // End of loop on equation
      ist++;
     } // End of loop on all equations used in the fit
 
     // Third loop is finished, now we update the correction matrices (see notes)
 
     Millepede::SpAVAt(clmat, clcmat, corrm, nalc, nagbn);
     Millepede::SpAX(clcmat, blvec, corrv, nalc, nagbn);
 
     for (i=0; i<nagbn; i++)
     {
      j = indbk[i];
      bgvec[j] -= corrv[i];
 
      for (k=0; k<nagbn; k++)
      {
           ik = indbk[k];
           cgmat[j][ik] -= corrm[i][k];
      }
     }
 
     indst.clear(); // Reset store for the next track 
     arest.clear();
 
     return true;
}

Referenced by MakeGlobalFit().

FitLoc() [2/2]

bool Millepede::FitLoc	(	int	n,
		double	track_params[],
		int	single_fit )

GetTrackNumber() [1/2]

int Millepede::GetTrackNumber

(

)

Definition at line 1591 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{return m_track_number;}

Referenced by CgemMilleAlign::fillHist(), MilleAlign::fillHist(), and MakeGlobalFit().

GetTrackNumber() [2/2]

int Millepede::GetTrackNumber

(

)

initialize() [1/2]

bool Millepede::initialize

(

)

Initialization.

initialize() [2/2]

bool Millepede::initialize

(

)

Initialization.

InitMille() [1/2]

bool Millepede::InitMille	(	bool	DOF[],
		double	Sigm[],
		int	nglo,
		int	nloc,
		double	startfact,
		int	nstd,
		double	res_cut,
		double	res_cut_init )

Definition at line 46 of file Mdc/MdcAlignAlg/MdcAlignAlg-00-01-04/src/Millepede.cxx.

{
 
     cout << "                                           " << endl;
     cout << "            * o o                   o      " << endl;
     cout << "              o o                   o      " << endl;
     cout << "   o ooooo  o o o  oo  ooo   oo   ooo  oo  " << endl;
     cout << "    o  o  o o o o o  o o  o o  o o  o o  o " << endl;
     cout << "    o  o  o o o o oooo o  o oooo o  o oooo " << endl;
     cout << "    o  o  o o o o o    ooo  o    o  o o    " << endl;
     cout << "    o  o  o o o o  oo  o     oo   ooo  oo  ++ starts" << endl;       
     cout << "                                           " << endl;
 
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "    Entering InitMille" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "-----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
 
     ncs = 0;
     loctot  = 0;                        // Total number of local fits
     locrej  = 0;                        // Total number of local fits rejected
     cfactref  = 1.0;                    // Reference value for Chi^2/ndof cut
 
     Millepede::SetTrackNumber(0);       // Number of local fits (starts at 0)
 
     m_residual_cut = res_cut;
     m_residual_cut_init = res_cut_init; 
 
     //   nagb    = 6*nglo;    // Number of global derivatives
     nagb    = 3*nglo;      // modified by wulh
     nalc     = nloc;       // Number of local derivatives
     nstdev  = nstd;     // Number of StDev for local fit chisquare cut
 
     if (debug_mode) cout << "Number of global parameters   : " << nagb << endl;
     if (debug_mode) cout << "Number of local parameters    : " << nalc << endl;
     if (debug_mode) cout << "Number of standard deviations : " << nstdev << endl;
 
     if (nagb>mglobl || nalc>mlocal)
     {
      if (debug_mode) cout << "Two many parameters !!!!!" << endl;
      return false;
     }
 
     // Global parameters initializations
 
     for (int i=0; i<nagb; i++)
     {
      bgvec[i]=0.;
      pparm[i]=0.;
      dparm[i]=0.;
      psigm[i]=-1.;
      indnz[i]=-1;
      indbk[i]=-1;
      nlnpa[i]=0;
 
      for (int j=0; j<nagb;j++)
      {
           cgmat[i][j]=0.;
      }
     }
 
     // Local parameters initializations
 
     for (int i=0; i<nalc; i++)
     {
      blvec[i]=0.;
 
      for (int j=0; j<nalc;j++)
      {
           clmat[i][j]=0.;
      }
     }
 
     // Then we fix all parameters...
 
     for (int j=0; j<nagb; j++)  {Millepede::ParSig(j,0.0);}
 
     // ...and we allow them to move if requested
 
     //   for (int i=0; i<3; i++)
     for (int i=0; i<3; i++)    // modified by wulh on 06/08/27
     {
      if (verbose_mode) cout << "GetDOF(" << i << ")= " << DOF[i] << endl;
 
      if (DOF[i]) 
      {
           for (int j=i*nglo; j<(i+1)*nglo; j++) 
           {Millepede::ParSig(j,Sigm[i]);}
      }
     }
 
     // Activate iterations (if requested)
 
     itert   = 0;   // By default iterations are turned off
     cfactr = 500.0;
     if (m_iteration) Millepede::InitUn(startfact);          
 
     arest.clear();  // Number of stored parameters when doing local fit
     arenl.clear();  // Is it linear or not?
     indst.clear(); 
 
     storeind.clear();
     storeare.clear();
     storenl.clear();
     storeplace.clear();
 
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "    InitMille has been successfully called!" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "-----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
 
     return true;
}

InitMille() [2/2]

bool Millepede::InitMille	(	bool	DOF[],
		double	Sigm[],
		int	nlay,
		int	nglo_on_lay,
		int	nloc,
		double	startfact,
		int	nstd,
		double	res_cut,
		double	res_cut_init )

Definition at line 47 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{
 
     cout << "                                           " << endl;
     cout << "            * o o                   o      " << endl;
     cout << "              o o                   o      " << endl;
     cout << "   o ooooo  o o o  oo  ooo   oo   ooo  oo  " << endl;
     cout << "    o  o  o o o o o  o o  o o  o o  o o  o " << endl;
     cout << "    o  o  o o o o oooo o  o oooo o  o oooo " << endl;
     cout << "    o  o  o o o o o    ooo  o    o  o o    " << endl;
     cout << "    o  o  o o o o  oo  o     oo   ooo  oo  ++ starts" << endl;       
     cout << "                                           " << endl;
 
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "    Entering InitMille" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "-----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
 
     ncs = 0;
     loctot  = 0;                        // Total number of local fits
     locrej  = 0;                        // Total number of local fits rejected
     //cfactref  = 1.0;                    // Reference value for Chi^2/ndof cut
     cfactref  = 3.0;                    // Reference value for Chi^2/ndof cut
     //cfactref  = 6.0;                    // Reference value for Chi^2/ndof cut
 
     record_mill = true;
     record_err = true;
 
     Millepede::SetTrackNumber(0);       // Number of local fits (starts at 0)
 
     m_residual_cut = res_cut;
     m_residual_cut_init = res_cut_init; 
     nglolay = nglo_on_lay;
     nagb    = nglo_on_lay*nlay;  // nglo stands for the number of layer, and 4 stands for the number of alignment parameters by guoaq - 2019/09/13
     nalc    = nloc;              // Number of local derivatives
     nstdev  = nstd;              // Number of StDev for local fit chisquare cut
 
     if (debug_mode) cout << "Number of global parameters   : " << nagb << endl;
     if (debug_mode) cout << "Number of local parameters    : " << nalc << endl;
     if (debug_mode) cout << "Number of standard deviations : " << nstdev << endl;
 
     if (nagb>mglobl || nalc>mlocal)
     {
      if (debug_mode) cout << "Two many parameters !!!!!" << endl;
      return false;
     }
 
     // Global parameters initializations
 
     for (int i=0; i<nagb; i++)
     {
      bgvec[i]=0.;
      pparm[i]=0.;
      dparm[i]=0.;
      psigm[i]=-1.;
      indnz[i]=-1;
      indbk[i]=-1;
      nlnpa[i]=0;
 
      for (int j=0; j<nagb;j++)
      {
           cgmat[i][j]=0.;
      }
     }
 
     // Local parameters initializations
 
     for (int i=0; i<nalc; i++)
     {
      blvec[i]=0.;
 
      for (int j=0; j<nalc;j++)
      {
           clmat[i][j]=0.;
      }
     }
 
     // Then we fix all parameters...
 
     for (int j=0; j<nagb; j++)  {Millepede::ParSig(j,0.0);}
 
     // ...and we allow them to move if requested
 
     //   for (int i=0; i<3; i++)
     for (int i=0; i<nglo_on_lay; i++)  // modified by wulh on 06/08/27 // 4 is the number of alginment parameters -- By guoaq 2019-09-14 
     {
      if (verbose_mode) cout << "GetDOF(" << i << ")= " << DOF[i] << endl;
 
      if (DOF[i]) 
      {
           for (int j=i*nlay; j<(i+1)*nlay; j++) 
           {
               Millepede::ParSig(j,Sigm[j]);
           } // Use layer dependent resolution, is it correct ? -- By guoaq 2019-10-11 
      }
     }
 
     // Activate iterations (if requested)
 
     itert   = 0;   // By default iterations are turned off
     cfactr = 500;    // cf stands for Chisq/ndof, ct stands for cut, what does a mean ?
                        // or cfact is the abbreviation of "cut factor" ? God knows ! -- guoaq-19-09-14
            // initial chisq / nof cut  =  cfactr / n_local_par
     
     cout<<"initial cfactr = "<<cfactr<<endl;
     if (m_iteration) Millepede::InitUn(startfact); // startfact is the initial chisq cut          
 
     arest.clear();  // Number of stored parameters when doing local fit
     arenl.clear();  // Is it linear or not?
     indst.clear(); 
 
     storeind.clear();
     storeare.clear();
     storenl.clear();
     storeplace.clear();
     
 
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "    InitMille has been successfully called!" << endl;
     if (debug_mode) cout << "" << endl;
     if (debug_mode) cout << "-----------------------------------------------------" << endl;
     if (debug_mode) cout << "" << endl;
 
     return true;
}

MakeGlobalFit() [1/2]

bool Millepede::MakeGlobalFit	(	double	par[],
		double	error[],
		double	pull[] )

Definition at line 851 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{
     int i, j, nf, nvar;
     int itelim = 0;
 
     int nstillgood;
 
     double sum;
 
     double step[150];
 
     double trackpars[2*mlocal];
 
     int ntotal_start, ntotal;
 
     cout << "..... Making global fit ....." << endl;
 
     ntotal_start = Millepede::GetTrackNumber();
 
     std::vector<double> track_slopes; // What is track_slopes? - By Guoaq-19-10-10
 
     track_slopes.resize(2*ntotal_start);
 
     for (i=0; i<2*ntotal_start; i++) track_slopes[i] = 0.;
 
     if (itert <= 1) itelim=10;    // Max number of iterations how to optimize this value - By Guoaq-19-09-26
 
     cout<<"itert, itelim = "<<itert<<", "<<itelim<<endl;
 
     while (itert < itelim)  // Iteration for the final loop
     {
      if (debug_mode) cout << "ITERATION --> " << itert << endl;
 
      ntotal = Millepede::GetTrackNumber();
      cout << "...using " << ntotal << " tracks..." << endl;
 
      // Start by saving the diagonal elements
 
      for (i=0; i<nagb; i++) {diag[i] = cgmat[i][i];}
 
      //  Then we retrieve the different constraints: fixed parameter or global equation
 
      nf = 0; // First look at the fixed global params
 
      for (i=0; i<nagb; i++)
      {
           if (psigm[i] <= 0.0)   // fixed global param
           {
            nf++;
 
            for (j=0; j<nagb; j++)
            {
             cgmat[i][j] = 0.0;  // Reset row and column
             cgmat[j][i] = 0.0;
            }           
           }
           else if (psigm[i] > 0.0) 
           {
            cgmat[i][i] += 1.0/(psigm[i]*psigm[i]);
           }
      }
 
      nvar = nagb;  // Current number of equations  
      if (debug_mode) cout << "Number of constraint equations : " << ncs << endl;
 
      if (ncs > 0) // Then the constraint equation
      {
           for (i=0; i<ncs; i++)
           {
            sum = arhs[i];
            for (j=0; j<nagb; j++)
            {
             cgmat[nvar][j] = float(ntotal)*adercs[i][j];
             cgmat[j][nvar] = float(ntotal)*adercs[i][j];          
             sum -= adercs[i][j]*(pparm[j]+dparm[j]);
 
            }
 
            cgmat[nvar][nvar] = 0.0;
            bgvec[nvar] = float(ntotal)*sum;
            nvar++;
           }
      }
 
 
      // Intended to compute the final global chisquare
 
      double final_cor = 0.0;
 
      if (itert > 1)
      {
           for (j=0; j<nagb; j++)
           {
            for (i=0; i<nagb; i++)
            {
             if (psigm[i] > 0.0)
             {
                  final_cor += step[j]*cgmat[j][i]*step[i]; 
                  if (i == j) final_cor -= step[i]*step[i]/(psigm[i]*psigm[i]);
             }
            }
           }
      }
 
      cout << " Final coeff is " << final_cor << endl;      
      cout << " Final NDOFs =  " << nagb << endl;
 
      //  The final matrix inversion
 
      nrank = Millepede::SpmInv(cgmat, bgvec, nvar, scdiag, scflag);
 
      for (i=0; i<nagb; i++)
      {
           dparm[i] += bgvec[i];    // Update global parameters values (for iterations)
           if (verbose_mode) cout << "dparm[" << i << "] = " << dparm[i] << endl;
           if (verbose_mode) cout << "cgmat[" << i << "][" << i << "] = " << cgmat[i][i] << endl;
           if (verbose_mode) cout << "err = " << sqrt(fabs(cgmat[i][i])) << endl;
 
           step[i] = bgvec[i];
 
           if (itert == 1) error[i] = cgmat[i][i]; // Unfitted error
      }
 
      cout << "" << endl;
      cout << "The rank defect of the symmetric " << nvar << " by " << nvar 
           << " matrix is " << nvar-nf-nrank << " (bad if non 0)" << endl;
 
      if (itert == 0)  break;       
      itert++;
 
      cout << "" << endl;
      cout << "Total : " << loctot << " local fits, " 
           << locrej << " rejected." << endl;
 
      // Reinitialize parameters for iteration
 
      loctot = 0;
      locrej = 0;
 
      if (cfactr != cfactref && sqrt(cfactr) > 1.2*cfactref)
      {
           cfactr = sqrt(cfactr);
      }
      else
      {
           cfactr = cfactref;
           //      itert = itelim;
      }
 
      if (itert == itelim)  break;  // End of story         
 
      cout << "Iteration " << itert << " with cut factor " << cfactr << endl;
 
      // Reset global variables
 
      for (i=0; i<nvar; i++)
      {
           bgvec[i] = 0.0;
           for (j=0; j<nvar; j++)
           {
            cgmat[i][j] = 0.0;
           }
      }
 
      //
      // We start a new iteration
      //
 
      // First we read the stores for retrieving the local params
 
      nstillgood = 0;   
 
      for (i=0; i<ntotal_start; i++)
      {
           int rank_i = 0;
           int rank_f = 0;
 
           (i>0) ? rank_i = abs(storeplace[i-1]) : rank_i = 0;
           rank_f = storeplace[i];
 
           if (verbose_mode) cout << "Track " << i << " : " << endl;
           if (verbose_mode) cout << "Starts at " << rank_i << endl;
           if (verbose_mode) cout << "Ends at " << rank_f << endl;
 
           if (rank_f >= 0) // Fit is still OK
           {
 
            if (itert > 3)
            {
             indst.clear();
             arest.clear();
 
             for (j=rank_i; j<rank_f; j++)
             {
                  indst.push_back(storeind[j]);
 
                  if (storenl[j] == 0) arest.push_back(storeare[j]);
                  if (storenl[j] == 1) arest.push_back(track_slopes[2*i]*storeare[j]);
                  if (storenl[j] == 2) arest.push_back(track_slopes[2*i+1]*storeare[j]);
             }
 
             for (j=0; j<2*nalc; j++) {trackpars[j] = 0.;}  
 
             Millepede::FitLoc(i,trackpars,1);
 
             track_slopes[2*i] = trackpars[2];
             track_slopes[2*i+1] = trackpars[6];
            }
            
            if (verbose_mode) cout << "after inters>3 " << endl;
 
            indst.clear();
            arest.clear();
 
 
            if (verbose_mode) cout << "rank_f is "<< rank_f << endl;
            for (j=rank_i; j<rank_f; j++)
            {
             if (verbose_mode) cout << "storeind  "<<j<<" is "<<storeind[j] << endl;
             if (verbose_mode) cout << "storearea  "<<j<<" is "<<storeare[j] << endl;
             if (verbose_mode) cout << "storenl  "<<j<<" is "<<storenl[j] << endl;
             
             indst.push_back(storeind[j]);
             if (verbose_mode) cout << "after 1st push_back  " << endl;
             if (storenl[j] == 0) arest.push_back(storeare[j]);
             if (verbose_mode) cout << "after 2nd push_back  " << endl;
             if (storenl[j] == 1) arest.push_back(track_slopes[2*i]*storeare[j]);
             if (verbose_mode) cout << "after 3rd push_back  " << endl;
             if (storenl[j] == 2) arest.push_back(track_slopes[2*i+1]*storeare[j]);
            }
            
            if (verbose_mode) cout << "after push_back arest " << endl;
 
            for (j=0; j<2*nalc; j++) {trackpars[j] = 0.;}   
 
            bool sc = Millepede::FitLoc(i,trackpars,0);
            
            if (verbose_mode) cout << "after invork FitLoc " << endl;
 
            (sc) 
             ? nstillgood++
             : storeplace[i] = -rank_f;     
           }
      } // End of loop on fits
 
      Millepede::SetTrackNumber(nstillgood);
 
     } // End of iteration loop
 
     // record the track number fit parameters and error
 
     Millepede::PrtGlo(); // Print the final results
 
     for (j=0; j<nagb; j++)
     {
      par[j]   = dparm[j];
      dparm[j] = 0.;
      pull[j]  = par[j]/sqrt(psigm[j]*psigm[j]-cgmat[j][j]);
      error[j] = sqrt(fabs(cgmat[j][j]));
     }
    
     if(record_err)
     {
         ofstream ferr("merr.csv",ios::app);
         ferr<<ntotal<<", "<< par[0]<<", "<< error[0]<<", "<<par[1]<<", "<< error[1]<<", "<<par[2]<<", "<< error[2]<<", "<<par[3]<<", "<< error[3]<<", "<<par[4]<<", "<< error[4]<<", "<<par[5]<<", "<< error[5]<<", "<<par[6]<<", "<< error[6]<<", "<<par[7] <<", "<< error[7]<<", "<<par[8]<<", "<< error[8]<<", "<<par[9]<<", "<< error[9]<<", "<<par[10]<<", "<< error[10]<<", "<<par[11]<<", "<< error[11]<< endl;
         ferr.close();
     }
 
     cout << std::setw(10) << " " << endl;
     cout << std::setw(10) << "            * o o                   o      " << endl;
     cout << std::setw(10) << "              o o                   o      " << endl;
     cout << std::setw(10) << "   o ooooo  o o o  oo  ooo   oo   ooo  oo  " << endl;
     cout << std::setw(10) << "    o  o  o o o o o  o o  o o  o o  o o  o " << endl;
     cout << std::setw(10) << "    o  o  o o o o oooo o  o oooo o  o oooo " << endl;
     cout << std::setw(10) << "    o  o  o o o o o    ooo  o    o  o o    " << endl;
     cout << std::setw(10) << "    o  o  o o o o  oo  o     oo   ooo  oo ++ ends." << endl;
     cout << std::setw(10) << "                       o                   " << endl;      
 
     return true;
}

MakeGlobalFit() [2/2]

bool Millepede::MakeGlobalFit	(	double	par[],
		double	error[],
		double	pull[] )

ParGlo() [1/2]

bool Millepede::ParGlo	(	int	index,
		double	param )

Definition at line 194 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{
     if (index<0 || index>=nagb)
     {return false;}
     else
     {pparm[index] = param;}
 
     return true;
}

ParGlo() [2/2]

bool Millepede::ParGlo	(	int	index,
		double	param )

ParSig() [1/2]

bool Millepede::ParSig	(	int	index,
		double	sigma )

Definition at line 221 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{
     if (index>=nagb) 
     {return false;}
     else
     {psigm[index] = sigma;}
 
     return true;
}

Referenced by InitMille(), and InitMille().

ParSig() [2/2]

bool Millepede::ParSig	(	int	index,
		double	sigma )

SetTrackNumber() [1/2]

void Millepede::SetTrackNumber

(

int

value

)

Definition at line 1592 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{m_track_number = value;}

Referenced by InitMille(), InitMille(), and MakeGlobalFit().

SetTrackNumber() [2/2]

void Millepede::SetTrackNumber

(

int

value

)

ZerLoc() [1/2]

bool Millepede::ZerLoc	(	double	dergb[],
		double	derlc[],
		double	dernl[] )

Definition at line 405 of file Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx.

{
     for(int i=0; i<nalc; i++) {derlc[i] = 0.0;}
     for(int i=0; i<nagb; i++) {dergb[i] = 0.0;}
     for(int i=0; i<nagb; i++) {dernl[i] = 0.0;}
 
     return true;
}

ZerLoc() [2/2]

bool Millepede::ZerLoc	(	double	dergb[],
		double	derlc[],
		double	dernl[] )

---

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

• Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/CgemAlignAlg/Millepede.h
• Mdc/MdcAlignAlg/MdcAlignAlg-00-01-04/MdcAlignAlg/Millepede.h
• Cgem/CgemAlignAlg/CgemAlignAlg-00-00-08/src/Millepede.cxx
• Mdc/MdcAlignAlg/MdcAlignAlg-00-01-04/src/Millepede.cxx