CgemLineFit Class Reference

#include <CgemLineFit.h>

Inheritance diagram for CgemLineFit:

Public Member Functions
	CgemLineFit (const std::string &name, ISvcLocator *pSvcLocator)
	~CgemLineFit ()
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
bool	Data_Max ()
bool	Data_Closest ()
bool	Loop_All ()
bool	Loop_MaxQ ()
bool	ToyMC ()
void	OrderClusters ()
void	OrderClusterSizeQ ()
void	FilterClusters ()
bool	Erase_outer ()
void	erase (int i)
int	GetVirLay (int geolay, double phi)
void	Discard (int layer)
StraightLine *	IniPar (double phi1, double z1, int i, double phi2, double z2, int j)
void	Filter (int layerid, StraightLine *l1)
void	Store_Trk (TMinuit *fit, int trkid)
TMinuit *	Fit (StraightLine *l_outer, int sheetflag, int typ)
bool	Get_MCInf ()
void	GetIntersection (StraightLine *line)
vector< double >	MC_truth ()
vector< double >	Get4Par (HepLorentzVector p4, Hep3Vector bp)
bool	HitPos (HepLorentzVector p4, Hep3Vector bp)
void	Get_OtherIniPar (int clst_0, int clst_1, int clst_2)
vector< double >	Get_Clusters (vector< double >Vec_truth)
void	Rec_Clusters ()
void	Rec_Clusters_mTPC ()
void	Rec_ClusterSize ()
void	Rec_ClusterQ ()
void	fireStripQ ()
void	fireStripT ()
void	fireStripID ()
void	Fit_Clusters (double par[])
void	InAngle (StraightLine sl)
double	RealPhi (double SimuPhi)
StatusCode	registerTrack (RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
void	Erase (int _layerid)
double	Min_dis (int R_id, double x, double z, StraightLine *l1)
double	Min_dis_up (int R_id, double x, double z, StraightLine *l1)
double	Min_dis_down (int R_id, double x, double z, StraightLine *l1)
bool	Layer_cross (int R_id, StraightLine *l1)
void	Swap (int i, int j)

Static Public Member Functions
static void	fcn (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static void	fcn2 (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static void	fcn3 (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static double	dx (double layerid, double R, double dr, double phi0, double z0, double tanl, double x)
static double	dV (int layer, double R, double dr, double phi0, double z0, double tglam, double x, double V)
static double	CalSigma2 (double layerid, double flag, double x, double V, double z, StraightLine *cosmic_ray)
static double	CalLineSAtY (StraightLine &l, double y)
static void	IncidentAngle (StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[])
static vector< double >	Trans (double par0, double par1, double par2, double par3)
static vector< int >	GetNMaxQ (vector< double > Q_11, vector< int > L_11, int Nmax)
static double	func (double layer, double R, double dr, double phi0, double z0, double tglam, double V, double x)
static void	to_0_2pi (double &arg)
static double	Min_diff (double arg1, double arg2)

Public Attributes
bool	debug
bool	MC
int	MAX_COMB
int	TEST_N
int	Nmax
int	Switch_CCmTPC
int	NXComb
int	NVComb
double	MinQ_Clus2D

Detailed Description

Definition at line 49 of file CgemLineFit.h.

Constructor & Destructor Documentation

CgemLineFit()

CgemLineFit::CgemLineFit	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 127 of file CgemLineFit.cxx.

                                                                        : Algorithm(name, pSvcLocator) {
    declareProperty("sigma", sigma = 0.13);              // unit:mm
    declareProperty("Run10_Flag", Run10_Flag = false);   // true:4 clusters, false: 3 clusters
    declareProperty("Align_Flag", Align_Flag = false);   // Alignment
    declareProperty("Switch_CCmTPC", Switch_CCmTPC = 0); // 0:Center charge 1:mTPC  2: merged position
    declareProperty("Method", Method = 0);      // 0: ToyMC 1: max_charge 2:closest to intersection  3: Loop_all  4: Loop_MaxQ
    declareProperty("MAX_COMB", MAX_COMB = 50); // The maximum number of combinations for method 3
    declareProperty("MAX_Clusters", MAX_Clusters = 500); // The maximum number of combinations for method 3
    declareProperty("MaxClu_Sheet", Nmax = 3);           // Max number of clusters on each sheet for method 4
    declareProperty("Chisq_cut", Chisq_cut = 2000);      // Max chisq for the selected set of clusters
    declareProperty("TEST_N", TEST_N = 0);               // set the sheet you want to study the efficiency / resolution. 1
    // to 6 from top to bottom. 0 use all the layers
    declareProperty("Debug", debug = 0);                 // Switch of debug
    declareProperty("MinQ_Clus2D", MinQ_Clus2D = 0);     // Add charge cut on 2D cluster to reduce the noise
    declareProperty("MC", MC = 0);                       // remove the hits on 3rd layer for MC
    declareProperty("LUTfile", m_LUTfile = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_17.root");
}

~CgemLineFit()

CgemLineFit::~CgemLineFit

(

)

Definition at line 144 of file CgemLineFit.cxx.

{ delete lutreader; }

Member Function Documentation

CalLineSAtY()

double CgemLineFit::CalLineSAtY ( StraightLine & l, double y )

static

Definition at line 4496 of file CgemLineFit.cxx.

                                                         {
    double drho = l.dr();
    double phi0 = l.phi0();
    double s = (y - drho * sin(phi0)) / sin(phi0 - TMath::Pi() / 2);
    return s;
}

CalSigma2()

double CgemLineFit::CalSigma2 ( double layerid, double flag, double x, double V, double z, StraightLine * cosmic_ray )

static

Definition at line 4401 of file CgemLineFit.cxx.

                                                                                                                 {
    double sigma = 0;
    double Ang[3];
    if (layerid == 0) {
        HepPoint3D l0(cos(x) * R_layer[0], sin(x) * R_layer[0], z);
        IncidentAngle(cosmic_ray, l0, pl_00->getStereoAngle(), Ang);
        if (flag == 0) {
            const int n = 28;
            double angle[n] = {-1.35, -1.25, -1.15, -1.05, -0.95, -0.85, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
                0.05,  0.15,  0.25,  0.35,  0.45,  0.55,  0.65,  0.75,  0.85,  0.95,  1.05,  1.15,  1.25,  1.35};
            double y[n] = {2.3576,   1.44448,  1.137,    0.998306,  0.819508,  0.695739, 0.583851, 0.487318, 0.408487, 0.317606,
                0.245717, 0.178443, 0.109751, 0.0469862, 0.0463951, 0.109757, 0.180819, 0.252827, 0.326331, 0.409922,
                0.482766, 0.594053, 0.690917, 0.835596,  0.93621,   1.25398,  1.56933,  1.93546};
            double *ex = nullptr;
            double ey[n] = {0.179035,    0.0641824,  0.0347424,  0.0256942,  0.015515,   0.0123214,  0.00893213,
                0.00584735,  0.00470714, 0.00334024, 0.00247131, 0.00171306, 0.00101272, 0.000571988,
                0.000525056, 0.00106782, 0.00170087, 0.00260604, 0.00347387, 0.00486819, 0.00638214,
                0.00910626,  0.0117511,  0.0172709,  0.0236933,  0.040096,   0.0607565,  0.171064};
            TGraphErrors gr(n, angle, y, ex, ey);
            sigma = gr.Eval(Ang[0]);
        }
        if (flag == 1) {
            const int n = 26;
            double angle[n] = {-1.25, -1.15, -1.05, -0.95, -0.85, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
                0.05,  0.15,  0.25,  0.35,  0.45,  0.55,  0.65,  0.75,  0.85,  0.95,  1.05,  1.15,  1.25};
            double y[n] = {1.32155,  1.00793,  0.844239, 0.745985,  0.607887,  0.525712, 0.449061, 0.372476, 0.305783,
                0.231835, 0.167037, 0.105903, 0.0526632, 0.0518158, 0.106611, 0.170067, 0.235356, 0.301801,
                0.374196, 0.46125,  0.526767, 0.61776,   0.749122,  0.848473, 1.01892,  1.50534};
            double *ex = nullptr;
            double ey[n] = {0.108111,    0.0523342,  0.0275105,  0.0227304, 0.0135576,   0.0095431,   0.00675217,
                0.00439062,  0.00328685, 0.00225977, 0.0013645, 0.000902477, 0.000515666, 0.000505206,
                0.000917112, 0.00148768, 0.00221462, 0.0031313, 0.0046925,   0.00745832,  0.00983847,
                0.0151778,   0.0227792,  0.033388,   0.0477401, 0.159787};
            TGraphErrors gr(n, angle, y, ex, ey);
            sigma = gr.Eval(Ang[1]);
        }
    }
    if (layerid == 1) {
        HepPoint3D l1(cos(x) * R_layer[1], sin(x) * R_layer[1], z);
        IncidentAngle(cosmic_ray, l1, pl_11->getStereoAngle(), Ang);
        if (flag == 0) {
            const int n = 14;
            double angle[n] = {-0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65};
            double y[n] = {0.459771,  0.381421, 0.310736, 0.238213, 0.169008, 0.103486, 0.0461908,
                0.0453169, 0.103585, 0.166698, 0.240332, 0.3125,   0.37973,  0.462382};
            double *ex = nullptr;
            double ey[n] = {0.00685022,  0.00437997,  0.00320362, 0.00214694, 0.00136347, 0.00081222, 0.000440626,
                0.000435363, 0.000818715, 0.00128686, 0.00214114, 0.00313311, 0.00436737, 0.0071195};
 
            TGraphErrors gr(n, angle, y, ex, ey);
            sigma = gr.Eval(Ang[0]);
        }
        if (flag == 1) {
            const int n = 16;
            double angle[n] = {-0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75};
            double y[n] = {0.543351,  0.449597, 0.376926, 0.299783, 0.234983, 0.165406, 0.105945, 0.0544188,
                0.0540844, 0.10723,  0.167146, 0.234566, 0.300303, 0.37238,  0.461751, 0.570197};
            double *ex = nullptr;
            double ey[n] = {0.0255211,   0.00933694,  0.00529557, 0.00317545, 0.00209458, 0.0013353,  0.000773631, 0.00046874,
                0.000453118, 0.000828085, 0.00125195, 0.00202037, 0.00315857, 0.00489179, 0.0100004,   0.032045};
            TGraphErrors gr(n, angle, y, ex, ey);
            sigma = gr.Eval(Ang[1]);
        }
    }
    /////
    if (layerid == 2) {
        HepPoint3D l2(cos(x) * R_layer[2], sin(x) * R_layer[2], z);
        IncidentAngle(cosmic_ray, l2, pl_21->getStereoAngle(), Ang);
        if (flag == 0) {
            const int n = 14;
            double angle[n] = {-0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65};
            double y[n] = {0.459771,  0.381421, 0.310736, 0.238213, 0.169008, 0.103486, 0.0461908,
                0.0453169, 0.103585, 0.166698, 0.240332, 0.3125,   0.37973,  0.462382};
            double *ex = nullptr;
            double ey[n] = {0.00685022,  0.00437997,  0.00320362, 0.00214694, 0.00136347, 0.00081222, 0.000440626,
                0.000435363, 0.000818715, 0.00128686, 0.00214114, 0.00313311, 0.00436737, 0.0071195};
 
            TGraphErrors gr(n, angle, y, ex, ey);
            sigma = gr.Eval(Ang[0]);
        }
        if (flag == 1) {
            const int n = 16;
            double angle[n] = {-0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75};
            double y[n] = {0.543351,  0.449597, 0.376926, 0.299783, 0.234983, 0.165406, 0.105945, 0.0544188,
                0.0540844, 0.10723,  0.167146, 0.234566, 0.300303, 0.37238,  0.461751, 0.570197};
            double *ex = nullptr;
            double ey[n] = {0.0255211,   0.00933694,  0.00529557, 0.00317545, 0.00209458, 0.0013353,  0.000773631, 0.00046874,
                0.000453118, 0.000828085, 0.00125195, 0.00202037, 0.00315857, 0.00489179, 0.0100004,   0.032045};
            TGraphErrors gr(n, angle, y, ex, ey);
            sigma = gr.Eval(Ang[1]);
        }
    }
    /////
    return sigma * sigma;
}

Referenced by fcn(), fcn2(), and fcn3().

Data_Closest()

bool CgemLineFit::Data_Closest

(

)

Definition at line 1180 of file CgemLineFit.cxx.

                               {
    double cl_00(0), cl_01(0), cl_11(0), cl_10(0);
    int _loop(0);
    double maxQ_11(0), maxQ_10(0), maxQ_00(0);
    int max_11(-1), max_10(-1), max_00(-1);
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
    RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
 
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        int layerid = cluster->getlayerid();
        double _Q = cluster->getenergydeposit();
        int sheetid = cluster->getsheetid();
        double _phi = cluster->getrecphi();
        if (flag != 2) continue;
        if (layerid == 1 && sheetid == 1 && maxQ_11 < _Q) {
            maxQ_11 = _Q;
            max_11 = _loop;
        }
        if (layerid == 1 && sheetid == 0 && maxQ_10 < _Q) {
            maxQ_10 = _Q;
            max_10 = _loop;
        }
        if (layerid == 1 && sheetid == 1) cl_11++;
        if (layerid == 0 && sheetid == 0 && _phi > 0) cl_01++;
        if (layerid == 0 && sheetid == 0 && _phi < 0) cl_00++;
        if (layerid == 1 && sheetid == 0) cl_10++;
    }
    clst_11 = cl_11;
    clst_01 = cl_01;
    clst_10 = cl_10;
    clst_00 = cl_00;
    _loop = 0;
    iter = recCgemClusterCol->begin();
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        if (flag != 2) continue;
 
        int layerid = cluster->getlayerid();
        int sheetid = cluster->getsheetid();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        double _Z = cluster->getRecZ();
        double _Q = cluster->getenergydeposit();
        double t_phi = cluster->getrecphi_mTPC();
        double t_v = cluster->getrecv_mTPC();
        double t_Z = cluster->getRecZ_mTPC();
        int vir_layerid = GetVirLay(layerid, _phi); // Added by Farinelli Riccardo in 2023.12.29
 
        RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
        RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
        double merged_phi = cluster_x->get_merge_phi();
        RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
        double merged_v = cluster_v->get_merge_v();
 
        if (_loop != max_11 && _loop != max_10 && layerid != 0) continue;
        Vec_flag.push_back(flag);
        Vec_layerid.push_back(layerid);
        Vec_sheetid.push_back(sheetid);
        Vec_layer.push_back(R_layer[layerid]);
        Vec_Virlayerid.push_back(vir_layerid); // Added by Farinelli Riccardo in 2023.12.29
 
        Vec_m_phi.push_back(t_phi);
        Vec_m_v.push_back(t_v);
        Vec_m_Z.push_back(t_Z);
 
        if (Switch_CCmTPC == 0) {
            Vec_phi.push_back(_phi);
            Vec_v.push_back(_v);
            Vec_Z.push_back(_Z);
        } else if (Switch_CCmTPC == 1) {
            Vec_phi.push_back(t_phi);
            Vec_v.push_back(t_v);
            Vec_Z.push_back(t_Z);
        } else if (Switch_CCmTPC == 2) {
 
            CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
            double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
            Vec_phi.push_back(merged_phi);
            Vec_v.push_back(merged_v);
            Vec_Z.push_back(merged_Z);
        }
 
        Vec_Q.push_back(_Q);
        cluster->setTrkId(0);
        vecclusters.push_back(cluster);
    }
    NC = vecclusters.size();
 
    if (Vec_layer.size() < 3) { return false; }
    for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
        if (Vec_layerid[i_cl] == 2) _clst_2++;
        if (Vec_layerid[i_cl] == 1) _clst_1++;
        if (Vec_layerid[i_cl] == 0) _clst_0++;
    }
 
    OrderClusters();
 
    if (_clst_1 == 2) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 1, Vec_phi[1], Vec_Z[1], 1);
    else return false;
 
    if (!Layer_cross(0, l_outer)) return false;
    if (_clst_0 > 1 || (Run10_Flag && _clst_0 > 2)) Filter(0, l_outer);
 
    return true;
}

Referenced by execute().

Data_Max()

bool CgemLineFit::Data_Max

(

)

Definition at line 857 of file CgemLineFit.cxx.

                           {
    if (debug) cout << "CgemLineFit::DataMax - START" << endl;
 
    double cl_00(0), cl_01(0), cl_11(0), cl_10(0), cl_20(0), cl_21(0);
    int _loop(0);
    double maxQ_11(0), maxQ_10(0), maxQ_00(0), maxQ_01(0), maxQ_21(0), maxQ_20(0);
    int max_11(-1), max_10(-1), max_00(-1), max_01(-1), max_21(-1), max_20(-1);
    double maxv_00(0), maxv_01(0);
 
    double phi_11(-1), phi_10(-1), phi_00(-1), phi_01(-1), phi_21(-1), phi_20(-1);
    double z_11(0), z_10(0), z_00(0), z_01(0), z_21(0), z_20(0);
    double cid_11(-1), cid_10(-1), cid_00(-1), cid_01(-1), cid_21(-1), cid_20(-1);
    double Xid_11(-1), Xid_10(-1), Xid_00(-1), Xid_01(-1), Xid_21(-1), Xid_20(-1);
    double Vid_11(-1), Vid_10(-1), Vid_00(-1), Vid_01(-1), Vid_21(-1), Vid_20(-1);
 
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
    RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
 
    int ic(0);
    int ic_tot(0);
    if (debug) cout << "CgemLineFit::DataMax - First loop on cluster" << endl;
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        if (debug) cout << _loop << endl;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        int layerid = cluster->getlayerid();
        double _Q = cluster->getenergydeposit();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        int sheetid = cluster->getsheetid();
 
        int clusterid = cluster->getclusterid();
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        double x = R_layer[layerid] * cos(_phi);
        double y = R_layer[layerid] * sin(_phi);
        double z = cluster->getRecZ();
        if (debug) cout << " flag: " << flag << endl;
        if (debug)
            cout << "CgemLineFit::DataMax - cluster info - ID: " << clusterid << " layer: " << layerid << " sheet: " << sheetid
                << " Xid: " << clusterflagb << " Vid: " << clusterflage << " Phi: " << _phi << " Z: " << z << " X: " << x
                << " Y: " << y << " Q: " << _Q << endl;
 
        if (_loop > MAX_Clusters) break;
        if (flag != 2) continue;        // only cluster 2D
        if (_Q < MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
 
        tr_id_cluster[ic] = clusterid;
        tr_x_cluster[ic] = x;
        tr_y_cluster[ic] = y;
        tr_z_cluster[ic] = z;
        tr_Q_cluster[ic] = _Q;
        // tr_phi_cluster[ic] = _phi;
        // tr_layer_cluster[ic] = layerid;
        // tr_sheet_cluster[ic] = sheetid;
 
        ic_tot++;
        ic++;
 
        if (layerid == 2 && sheetid == 1 && maxQ_21 < _Q) {
            maxQ_21 = _Q;
            max_21 = _loop;
            phi_21 = _phi;
            z_21 = z;
            cid_21 = clusterid;
            Xid_21 = clusterflagb;
            Vid_21 = clusterflage;
        }
        if (layerid == 2 && sheetid == 0 && maxQ_20 < _Q) {
            maxQ_20 = _Q;
            max_20 = _loop;
            phi_20 = _phi;
            z_20 = z;
            cid_20 = clusterid;
            Xid_20 = clusterflagb;
            Vid_20 = clusterflage;
        }
        if (layerid == 1 && sheetid == 1 && maxQ_11 < _Q) {
            maxQ_11 = _Q;
            max_11 = _loop;
            phi_11 = _phi;
            z_11 = z;
            cid_11 = clusterid;
            Xid_11 = clusterflagb;
            Vid_11 = clusterflage;
        }
        if (layerid == 1 && sheetid == 0 && maxQ_10 < _Q) {
            maxQ_10 = _Q;
            max_10 = _loop;
            phi_10 = _phi;
            z_10 = z;
            cid_10 = clusterid;
            Xid_10 = clusterflagb;
            Vid_10 = clusterflage;
        }
        if (layerid == 0 && sheetid == 0 && _phi > 0 && maxQ_01 < _Q && maxv_00 != _v) {
            maxQ_01 = _Q;
            max_01 = _loop;
            phi_01 = _phi;
            z_01 = z;
            cid_01 = clusterid;
            Xid_01 = clusterflagb;
            Vid_01 = clusterflage;
        }
        if (layerid == 0 && sheetid == 0 && _phi < 0 && maxQ_00 < _Q && maxv_00 != _v) {
            maxQ_00 = _Q;
            max_00 = _loop;
            phi_00 = _phi;
            z_00 = z;
            cid_00 = clusterid;
            Xid_00 = clusterflagb;
            Vid_00 = clusterflage;
        }
 
        if (layerid == 2 && sheetid == 1) cl_21++;
        if (layerid == 2 && sheetid == 0) cl_20++;
        if (layerid == 1 && sheetid == 1) cl_11++;
        if (layerid == 1 && sheetid == 0) cl_10++;
        if (layerid == 0 && sheetid == 0 && _phi > 0) cl_01++;
        if (layerid == 0 && sheetid == 0 && _phi < 0) cl_00++;
    }
 
    if (debug) cout << "CgemLineFit::DataMax - First loop on cluster ended - number of cluster: " << ic << endl;
 
    tr_ncluster = ic;
    for (int i = 0; i < ic; i++) {
        tr_Is_selected_cluster[i] = 0;
        if (i == (max_00 - 1) && maxQ_00 > 0) tr_Is_selected_cluster[i] = 1;
        if (i == (max_01 - 1) && maxQ_01 > 0) tr_Is_selected_cluster[i] = 1;
        if (i == (max_10 - 1) && maxQ_10 > 0) tr_Is_selected_cluster[i] = 1;
        if (i == (max_11 - 1) && maxQ_11 > 0) tr_Is_selected_cluster[i] = 1;
        if (i == (max_20 - 1) && maxQ_20 > 0) tr_Is_selected_cluster[i] = 1;
        if (i == (max_21 - 1) && maxQ_21 > 0) tr_Is_selected_cluster[i] = 1;
    }
 
    clst_21 = cl_21;
    clst_20 = cl_20;
    clst_11 = cl_11;
    clst_01 = cl_01;
    clst_10 = cl_10;
    clst_00 = cl_00;
    _loop = 0;
 
    if (debug)
        cout << "CgemLineFit::DataMax - number of cluster per sheet/layer: " << clst_00 << " " << clst_10 << " " << clst_01 << " "
            << clst_11 << " " << clst_20 << " " << clst_21 << endl;
 
    iter = recCgemClusterCol->begin();
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        if (flag != 2) continue;
 
        int layerid = cluster->getlayerid();
        int sheetid = cluster->getsheetid();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        double _Z = cluster->getRecZ();
        double _Q = cluster->getenergydeposit();
        double t_phi = cluster->getrecphi_mTPC();
        double t_v = cluster->getrecv_mTPC();
        double t_Z = cluster->getRecZ_mTPC();
        int vir_layerid = GetVirLay(layerid, _phi);
 
        RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
        RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
        double merged_phi = cluster_x->get_merge_phi();
        RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
        double merged_v = cluster_v->get_merge_v();
 
        // Added to be checked - farinelli 20240106
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        RecCgemClusterCol::iterator id_XCluster = recCgemClusterCol->begin() + clusterflagb;
        RecCgemCluster *xcluster = (*id_XCluster);
        int xclustersize = abs(xcluster->getclusterflage() - xcluster->getclusterflagb()) + 1;
        double X_Q = xcluster->getenergydeposit();
        RecCgemClusterCol::iterator id_VCluster = recCgemClusterCol->begin() + clusterflage;
        RecCgemCluster *vcluster = (*id_VCluster);
        int vclustersize = abs(vcluster->getclusterflage() - vcluster->getclusterflagb()) + 1;
        double V_Q = vcluster->getenergydeposit();
        int xclusterb = xcluster->getclusterflagb();
        int xclustere = xcluster->getclusterflage();
        int vclusterb = vcluster->getclusterflagb();
        int vclustere = vcluster->getclusterflage();
        int layerid_str, sheetid_str, stripid_str, time_chnnel;
        int nxstrips = 0;
        int nvstrips = 0;
        double energydeposit;
        double Q_fC, T_ns, Tf_ns;
        SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(), "/Event/Digi/CgemDigiCol");
        CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
        Vec_x_stripsQ.clear();
        Vec_v_stripsQ.clear();
        Vec_x_stripsT.clear();
        Vec_v_stripsT.clear();
        Vec_x_stripsTf.clear();
        Vec_v_stripsTf.clear();
        Vec_x_stripsid.clear();
        Vec_v_stripsid.clear();
        for (; iter_digi != cgemDigiCol->end(); iter_digi++) {
            layerid_str = CgemID::layer((*iter_digi)->identify());
            sheetid_str = CgemID::sheet((*iter_digi)->identify());
            stripid_str = CgemID::strip((*iter_digi)->identify());
            bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
            int flagxv;
            if (striptype == true) flagxv = 0;
            else flagxv = 1;
            if (layerid == layerid_str && sheetid == sheetid_str) {
                Q_fC = (*iter_digi)->getCharge_fc();
                T_ns = (*iter_digi)->getTime_ns();
                Tf_ns = get_Time(iter_digi);
                if (flagxv == 0) {
                    if (stripid_str >= xclusterb && stripid_str <= xclustere) {
                        if (xclustersize < 100) {
                            Vec_x_stripsid.push_back(stripid_str);
                            Vec_x_stripsQ.push_back(Q_fC);
                            Vec_x_stripsT.push_back(T_ns);
                            Vec_x_stripsTf.push_back(Tf_ns);
                            nxstrips++;
                        }
                    }
                }
                if (flagxv == 1) {
                    if (stripid_str >= vclusterb && stripid_str <= vclustere) {
                        if (vclustersize < 100) {
                            Vec_v_stripsid.push_back(stripid_str);
                            Vec_v_stripsQ.push_back(Q_fC);
                            Vec_v_stripsT.push_back(T_ns);
                            Vec_v_stripsTf.push_back(Tf_ns);
                            nvstrips++;
                        }
                    }
                }
            }
        }
        // End to be checked
 
        if (debug) cout << "cl ID: " << _loop << endl;
        if (!(_loop == max_00 || _loop == max_10 || _loop == max_01 || _loop == max_11 || _loop == max_20 || _loop == max_21))
            continue;
        if (debug) cout << "is selected" << endl;
 
        Vec_flag.push_back(flag);
        Vec_layerid.push_back(layerid);
        Vec_sheetid.push_back(sheetid);
        Vec_layer.push_back(R_layer[layerid]);
        Vec_Virlayerid.push_back(vir_layerid);
 
        Vec_m_phi.push_back(t_phi);
        Vec_m_v.push_back(t_v);
        Vec_m_Z.push_back(t_Z);
 
        Vec_xCluSize.push_back(xclustersize);
        Vec_vCluSize.push_back(vclustersize);
        Vec_XQ.push_back(X_Q);
        Vec_VQ.push_back(V_Q);
        Vec_XstripQ.push_back(Vec_x_stripsQ);
        Vec_VstripQ.push_back(Vec_v_stripsQ);
        Vec_XstripT.push_back(Vec_x_stripsT);
        Vec_VstripT.push_back(Vec_v_stripsT);
        Vec_XstripTf.push_back(Vec_x_stripsTf);
        Vec_VstripTf.push_back(Vec_v_stripsTf);
        Vec_XstripID.push_back(Vec_x_stripsid);
        Vec_VstripID.push_back(Vec_v_stripsid);
 
        if (Switch_CCmTPC == 0) {
            Vec_phi.push_back(_phi);
            Vec_v.push_back(_v);
            Vec_Z.push_back(_Z);
 
        } else if (Switch_CCmTPC == 1) {
            Vec_phi.push_back(t_phi);
            Vec_v.push_back(t_v);
            Vec_Z.push_back(t_Z);
 
        } else if (Switch_CCmTPC == 2) {
 
            CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
            double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
            Vec_phi.push_back(merged_phi);
            Vec_v.push_back(merged_v);
            Vec_Z.push_back(merged_Z);
        }
 
        Vec_Q.push_back(_Q);
        cluster->setTrkId(0);
        vecclusters.push_back(cluster);
    }
    NC = vecclusters.size();
 
    if (debug) cout << "CgemLineFit::DataMax - Number of points: " << Vec_layer.size() << endl;
 
    // if(Vec_layer.size()!=6) return false;
    if (Vec_layer.size() < 2) return false;
    NXComb = 1;
    NVComb = 1;
 
    for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
        if (Vec_layerid[i_cl] == 2) _clst_2++;
        if (Vec_layerid[i_cl] == 1) _clst_1++;
        if (Vec_layerid[i_cl] == 0) _clst_0++;
    }
 
    OrderClusters();
 
    // if(debug) cout<<"CgemLineFit::DataMax - _clst_2: "<<_clst_2<<endl;
    // if(_clst_1==2) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
    // changed to be checked - farinelli 20240106
    if (debug) cout << "CgemLineFit::DataMax - Vec_layerid.size(): " << Vec_layerid.size() << endl;
    if (Vec_layerid.size() >= 2)
        l_outer = IniPar(Vec_phi[0], Vec_Z[0], 2 * Vec_layerid[0], Vec_phi[1], Vec_Z[1], 2 * Vec_layerid[1]);
    else return false;
 
    return true;
}

Referenced by execute().

Discard()

void CgemLineFit::Discard

(

int

layer

)

Definition at line 2985 of file CgemLineFit.cxx.

                                   {
    for (int i = 0; i < Vec_layer.size(); i++) {
        if (Vec_layerid[i] == layer) erase(i);
    }
}

Referenced by Filter().

dV()

double CgemLineFit::dV ( int layer, double R, double dr, double phi0, double z0, double tglam, double x, double V )

static

Definition at line 2720 of file CgemLineFit.cxx.

                                                                                                                     {
 
    HepPoint3D Up, Down;
    StraightLine l1(dr, phi0, z0, tglam);
    double phiV_up[2], phiV_down[2];
    bool findinter = false;
    if (!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
    else findinter = Mp->getPointAligned(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
 
    if (!findinter) return 9999999999999;
 
    double dx1 = Min_diff(x, phiV_down[0]);
    double dx2 = Min_diff(x, phiV_up[0]);
 
    if (dx1 < dx2) return pow(fabs(phiV_down[1] - V), 2.0);
    else return pow(fabs(phiV_up[1] - V), 2.0);
}

Referenced by fcn2(), and fcn3().

dx()

double CgemLineFit::dx ( double layerid, double R, double dr, double phi0, double z0, double tanl, double x )

static

Definition at line 2695 of file CgemLineFit.cxx.

                                                                                                             {
    if (dr > R) return 9999999999999;
    StraightLine l1(dr, phi0, dz, tanl);
    double phiV_up[2], phiV_down[2];
    HepPoint3D Up, Down;
    double phiV_up1_old[2], phiV_down1_old[2];
    HepPoint3D Up1_old, Down1_old;
    bool findinter = false;
    if (!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
    else findinter = Mp->getPointAligned(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
 
    if (!findinter) return 9999999999999;
 
    phiV_up[0] = fmod((phiV_up[0]), (2.0 * TMath::Pi()));
    phiV_down[0] = fmod((phiV_down[0]), (2.0 * TMath::Pi()));
 
    if (phiV_up[0] < 0) phiV_up[0] += 2.0 * TMath::Pi();
    if (phiV_down[0] < 0) phiV_down[0] += 2.0 * TMath::Pi();
 
    double dx1 = Min_diff(x, phiV_down[0]);
    double dx2 = Min_diff(x, phiV_up[0]);
    double mdv = min(pow(dx1, 2.0), pow(dx2, 2.0));
    return mdv * R * R;
}

Referenced by fcn(), and fcn3().

Erase()

void CgemLineFit::Erase

(

int

_layerid

)

erase()

void CgemLineFit::erase

(

int

i

)

Definition at line 2955 of file CgemLineFit.cxx.

                             {
    vector<double>::iterator _iter_v;
    _iter_v = Vec_v.begin();
    Vec_v.erase(_iter_v + i);
 
    vector<double>::iterator _iter_phi;
    _iter_phi = Vec_phi.begin();
    Vec_phi.erase(_iter_phi + i);
 
    vector<double>::iterator _iter_Z;
    _iter_Z = Vec_Z.begin();
    Vec_Z.erase(_iter_Z + i);
 
    vector<double>::iterator _iter_layer;
    _iter_layer = Vec_layer.begin();
    Vec_layer.erase(_iter_layer + i);
 
    vector<double>::iterator _iter_layerid;
    _iter_layerid = Vec_layerid.begin();
    Vec_layerid.erase(_iter_layerid + i);
 
    vector<double>::iterator _iter_Virlayerid;
    _iter_Virlayerid = Vec_Virlayerid.begin();
    Vec_Virlayerid.erase(_iter_Virlayerid + i);
 
    vector<double>::iterator _iter_flag;
    _iter_flag = Vec_flag.begin();
    Vec_flag.erase(_iter_flag + i);
}

Referenced by Discard(), Erase_outer(), and FilterClusters().

Erase_outer()

bool CgemLineFit::Erase_outer

(

)

Definition at line 2934 of file CgemLineFit.cxx.

                              {
 
    vector<double> _x, _v;
    vector<int> _iter;
    int wrong(-1);
    for (int i = 0; i < 3; i++) {
        _x.push_back(Vec_phi[i]);
        _v.push_back(Vec_v[i]);
        _iter.push_back(i);
    }
    if (_x[0] == _x[1] && _v[0] == _v[2]) wrong = 0;
    else if (_x[0] == _x[1] && _v[1] == _v[2]) wrong = 1;
    else if (_x[0] == _x[2] && _v[0] == _v[1]) wrong = 0;
    else if (_x[0] == _x[2] && _v[2] == _v[1]) wrong = 2;
    else if (_x[1] == _x[2] && _v[1] == _v[0]) wrong = 1;
    else if (_x[1] == _x[2] && _v[2] == _v[0]) wrong = 2;
    else { return false; }
    erase(_iter[wrong]);
    return true;
}

execute()

StatusCode CgemLineFit::execute

(

)

Definition at line 561 of file CgemLineFit.cxx.

                                {
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in execute()" << endreq;
    if (debug) cout << "CgemLineFit::execute START" << endl;
 
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(), "/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::FATAL << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
 
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
    if (!recCgemClusterCol) {
        log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
        return StatusCode::FAILURE;
    }
    SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(), "/Event/Digi/CgemDigiCol");
    if (!cgemDigiCol) {
        log << MSG::WARNING << "Could not retrieve Cgem digi list" << endreq;
        return StatusCode::FAILURE;
    }
 
    if (debug) cout << "CgemLineFit::execute - read the collection" << endl;
 
    m_event = eventHeader->eventNumber();
    m_run = eventHeader->runNumber();
    event = m_event;
    run = m_run;
    tr_run = run;
    tr_event = event;
    if ((m_event % 500) == 0) cout << " begin evt : " << event << endl;
    //cout<<__FILE__<<": evt = "<<event<<endl;
    //////////////////Collecting Clusters///////////
 
    // release the memory
    Vec_layer.clear();
    Vec_phi.clear();
    Vec_Z.clear();
    Vec_v.clear();
    Vec_layerid.clear();
    Vec_Virlayerid.clear();
    Vec_flag.clear();
    Vec_Q.clear();
    Vec_sheetid.clear();
    Vec_m_layer.clear();
    Vec_m_phi.clear();
    Vec_m_Z.clear();
    Vec_m_v.clear();
    Vec_m_layerid.clear();
    Vec_m_flag.clear();
    Vec_m_Q.clear();
    Vec_m_sheetid.clear();
    vecclusters.clear();
    SmartRefVector<RecCgemCluster>().swap(vecclusters);
 
    _clst_0 = 0;
    _clst_1 = 0;
    _clst_2 = 0;
 
    NXComb = -1;
    NVComb = -1;
    // int is_cluster_exsited[6] = {0,0,0,0,0,0};
    // int nEvent_with_6cluster = 0;
    // int nEvent_with_4cluster = 0;
    // for(int i=0;i<6;i++){is_cluster_exsited[i]=0;}
 
    if (debug) cout << "CgemLineFit::execute - cleaned the variables" << endl;
 
    TStopwatch gtimer;
    gtimer.Start();
    if (debug) cout << "start Straight line fit" << endl;
    if (Method == 0) {
        if (!ToyMC()) {
            tr_Is_fitted = 0;
            tr_chisq = -1;
            m_tuple_track->write();
            if (debug) cout << "CgemLineFit::execute - method 0 failed" << endl;
            return StatusCode::SUCCESS;
        }
    } else if (Method == 1) {
        if (!Data_Max()) {
            tr_Is_fitted = 0;
            tr_chisq = -1;
            m_tuple_track->write();
            if (debug) cout << "CgemLineFit::execute - method 1 failed" << endl;
            return StatusCode::SUCCESS;
        }
    } else if (Method == 2) {
        if (!Data_Closest()) {
            tr_Is_fitted = 0;
            tr_chisq = -1;
            m_tuple_track->write();
            if (debug) cout << "CgemLineFit::execute - method 2 failed" << endl;
            return StatusCode::SUCCESS;
        }
    } else if (Method == 3) {
        if (!Loop_All()) {
            tr_Is_fitted = 0;
            tr_chisq = -1;
            m_tuple_track->write();
            if (debug) cout << "CgemLineFit::execute - method 3 failed" << endl;
            return StatusCode::SUCCESS;
        }
    } else if (Method == 4) {
        if (!Loop_MaxQ()) {
            tr_Is_fitted = 0;
            tr_chisq = -1;
            m_tuple_track->write();
            // cout << "Nevent with 6 clusters: " << nEvent_with_6cluster << endl;
            // cout << "Nevent with 4 clusters: " << nEvent_with_4cluster << endl;
            if (debug) cout << "CgemLineFit::execute - method 4 failed" << endl;
            return StatusCode::SUCCESS;
        }
        // cout << "Nevent with 6 clusters: " << nEvent_with_6cluster << endl;
        // cout << "Nevent with 4 clusters: " << nEvent_with_4cluster << endl;
    }
 
    if (Vec_layerid.size() < 4 && Run10_Flag) return StatusCode::SUCCESS;
 
    /////Set initial parameter according to the clusters on the outtermost layer and perform the final fit:
    /////1, fit phi0, drho; 2, fit tanl z0; 3, fit 4 par. at the same time/////
 
    if (debug) cout << "CgemLineFit::execute - set initial parameter" << endl;
 
    TMinuit *fit = Fit(l_outer, fa, 0);
    double trkpar[4] = {-999, -999, -999, -999};
    double trkparErr[4] = {-999, -999, -999, -999};
    Int_t ierflg, npari, nparx, istat3;
    Double_t fmin, edm, errdef;
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat3);
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    double emat[16];
    fit->mnemat(emat, 4);
    chisq_3 = fit->fAmin;
    registerTrack(m_trackList_tds, m_hitList_tds);
    Store_Trk(fit, 0);
    delete fit;
    double chisq_last = chisq_3;
 
    gtimer.Stop();
    vector<double> _trkpar = Trans(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
 
    par0 = _trkpar[0];
    par1 = _trkpar[1];
    par2 = _trkpar[2];
    par3 = _trkpar[3];
    Err_par0 = trkparErr[0];
    Err_par1 = trkparErr[1];
    Err_par2 = trkparErr[2];
    Err_par3 = trkparErr[3];
 
    status3 = istat3;
 
    clst_0 = _clst_0;
    clst_1 = _clst_1;
    clst_2 = _clst_2;
 
    R0 = R_layer[0];
    R1 = R_layer[1];
    R2 = R_layer[2];
 
    // the initial track parameters based on the clusters on the outtermost layer
    ini_0 = l_outer->dr();
    ini_1 = l_outer->phi0();
    ini_2 = l_outer->dz();
    ini_3 = l_outer->tanl();
 
    if (debug) cout << "Get_otherIniPar" << endl;
    // Get the initial track parameters based on the clusters on middle and innermost layers
    if (TEST_N == 0) Get_OtherIniPar(_clst_0, _clst_1, _clst_2);
 
    // Store the x,v,z information of clusters on the avialable layers,_clst_0: N clusters on middle layer, _clst_1: N clusters on
    // innermost layer _clst_0=2; _clst_1=2; _clst_2=0;
 
    if (debug) cout << "Rec_Clusters" << endl;
    Rec_Clusters();
    Rec_ClusterQ();
    Rec_ClusterSize();
    fireStripQ();
    fireStripT();
    fireStripID();
    // Store the x,v,z information of clusters by mTPC method on the avialable layers
    if (Switch_CCmTPC == 1) Rec_Clusters_mTPC();
 
    if (debug) cout << "Fit_Clusters" << endl;
    // Store the intersections between fitted line and Cgem
    Fit_Clusters(trkpar);
 
    if (debug) cout << "GetIntersection" << endl;
    // Store the intersections according to the extroplation of the clusters on other layers
    StraightLine interline(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    GetIntersection(&interline);
 
    if (debug) cout << "GetRes" << endl;
    // Get the resolution
    StraightLine l_fit(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    GetRes(&l_fit, TEST_N);
    m_tuple->write();
 
    tr_drho = trkpar[0];
    tr_phi0 = trkpar[1];
    tr_dz0 = trkpar[2];
    tr_tanLam = trkpar[3];
    tr_Is_fitted = 1;
    tr_chisq = chisq_last;
    m_tuple_track->write();
 
    delete l_outer;
    if (debug) cout << "CgemLineFit::execute END" << endl;
 
    ////////
    //cout << "NC = " << NC << endl;
    //SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
    //if (!recMdcTrackCol) {
    //  log << MSG::WARNING << "Could not retrieve Rec Track list" << endreq;
    //  return StatusCode::FAILURE;
    //}
 
    //RecMdcTrackCol::iterator itera = recMdcTrackCol->begin();
    //cout << "Ncluster = " << (*itera)->getNcluster() << endl;
    /////////
    return StatusCode::SUCCESS;
}

fcn()

void CgemLineFit::fcn ( Int_t & npar, Double_t * gin, Double_t & f, Double_t * par, Int_t iflag )

static

Definition at line 2623 of file CgemLineFit.cxx.

                                                                                         {
    double chisq = 0;
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
 
        if (Vec_Virlayerid[i] == 5 && bool(sheet_flag & f21)) continue;
        if (Vec_Virlayerid[i] == 4 && bool(sheet_flag & f20)) continue;
        if (Vec_Virlayerid[i] == 3 && bool(sheet_flag & f11)) continue;
        if (Vec_Virlayerid[i] == 2 && bool(sheet_flag & f10)) continue;
        if (Vec_Virlayerid[i] == 1 && bool(sheet_flag & f01)) continue;
        if (Vec_Virlayerid[i] == 0 && bool(sheet_flag & f00)) continue;
        double dx1 = dx(Vec_Virlayerid[i], Vec_layer[i], par[0], par[1], dz_set, tanl_set, Vec_phi[i]);
        StraightLine *l = new StraightLine(par[0], par[1], dz_set, tanl_set);
        double flagg = Vec_flag[i];
        if (flagg == 2) flagg = 0;
        double _sigma2 = CalSigma2(Vec_layerid[i], flagg, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
        chisq += dx1 / _sigma2;
        delete l;
    }
    f = chisq;
}

Referenced by Fit().

fcn2()

void CgemLineFit::fcn2 ( Int_t & npar, Double_t * gin, Double_t & f, Double_t * par, Int_t iflag )

static

Definition at line 2644 of file CgemLineFit.cxx.

                                                                                          {
    double chisq = 0;
 
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5 && bool(sheet_flag & f21)) continue;
        if (Vec_Virlayerid[i] == 4 && bool(sheet_flag & f20)) continue;
        if (Vec_Virlayerid[i] == 3 && bool(sheet_flag & f11)) continue;
        if (Vec_Virlayerid[i] == 2 && bool(sheet_flag & f10)) continue;
        if (Vec_Virlayerid[i] == 1 && bool(sheet_flag & f01)) continue;
        if (Vec_Virlayerid[i] == 0 && bool(sheet_flag & f00)) continue;
 
        if (Vec_flag[i] != 2) continue;
        double flagg = Vec_flag[i];
        if (flagg == 2) flagg = 1;
        double dv1 = dV(Vec_Virlayerid[i], Vec_layer[i], dr_set, phi0_set, par[0], par[1], Vec_phi[i], Vec_v[i]);
        StraightLine *l = new StraightLine(dr_set, phi0_set, par[0], par[1]);
        // double _sigma2 = CalSigma2(Vec_layerid[i], Vec_flag[i], Vec_phi[i], Vec_v[i], Vec_Z[i], l);
        double _sigma2 = CalSigma2(Vec_layerid[i], flagg, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
        chisq += dv1 / _sigma2;
        delete l;
    }
    f = chisq;
}

Referenced by Fit().

fcn3()

void CgemLineFit::fcn3 ( Int_t & npar, Double_t * gin, Double_t & f, Double_t * par, Int_t iflag )

static

Definition at line 2668 of file CgemLineFit.cxx.

                                                                                          {
 
    double chisq = 0;
    double dv1(0), dx1(0);
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_flag[i] == 2) {
 
            if (Vec_Virlayerid[i] == 5 && bool(sheet_flag & f21)) continue;
            if (Vec_Virlayerid[i] == 4 && bool(sheet_flag & f20)) continue;
            if (Vec_Virlayerid[i] == 3 && bool(sheet_flag & f11)) continue;
            if (Vec_Virlayerid[i] == 2 && bool(sheet_flag & f10)) continue;
            if (Vec_Virlayerid[i] == 1 && bool(sheet_flag & f01)) continue;
            if (Vec_Virlayerid[i] == 0 && bool(sheet_flag & f00)) continue;
 
            dv1 = dV(Vec_Virlayerid[i], Vec_layer[i], par[0], par[1], par[2], par[3], Vec_phi[i], Vec_v[i]);
            StraightLine *l = new StraightLine(par[0], par[1], par[2], par[3]);
            double _sigma2 = CalSigma2(Vec_layerid[i], 1, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
            chisq += (dv1 / _sigma2);
            _sigma2 = CalSigma2(Vec_layerid[i], 0, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
            dx1 = dx(Vec_Virlayerid[i], Vec_layer[i], par[0], par[1], par[2], par[3], Vec_phi[i]);
            chisq += dx1 / _sigma2;
            delete l;
        }
    }
    f = chisq;
}

Referenced by Fit().

Filter()

void CgemLineFit::Filter	(	int	layerid,
		StraightLine *	l1 )

Definition at line 3008 of file CgemLineFit.cxx.

                                                      {
    if (!Layer_cross(layerid, l1)) { Discard(layerid); }
 
    vector<int> number;
    vector<double> dst;
    vector<double> dst_up;
    vector<double> dst_down;
    vector<double> v_x;
    vector<double> v_v;
 
    if (Run10_Flag) {
        for (int i = 0; i < Vec_layer.size(); i++) {
            if (Vec_layerid[i] != layerid) continue;
            number.push_back(i);
            // obtain the distance between the up/down side of crosssections to the given position (phi, V)
            dst_up.push_back(Min_dis_up(layerid, Vec_phi[i], Vec_v[i], l1));
            dst_down.push_back(Min_dis_down(layerid, Vec_phi[i], Vec_v[i], l1));
        }
 
        for (int j = 0; j < number.size(); j++) { int i = number[j]; }
 
        // sort the clusters according to the up side distance
        int Loop;
        Loop = 0;
        while (Loop != number.size() - 1) {
            Loop = 0;
            for (int j = 0; j < number.size() - 1; j++) {
                if (dst_up[j] > dst_up[j + 1]) {
                    swap(dst_up[j], dst_up[j + 1]);
                    swap(dst_down[j], dst_down[j + 1]);
                    swap(number[j], number[j + 1]);
                } else Loop++;
            }
        }
 
        // set the index of the nearest cluster
        for (int j = 0; j < number.size(); j++) { int i = number[j]; }
 
        // remove the first member
        number.erase(number.begin());
        dst_up.erase(dst_up.begin());
        dst_down.erase(dst_down.begin());
 
        // what's the difference to previous operation?
        for (int j = 0; j < number.size(); j++) { int i = number[j]; }
 
        // sort the clusters according to the down side distance
        Loop = 0;
        while (Loop != number.size() - 1) {
            Loop = 0;
            for (int j = 0; j < number.size() - 1; j++) {
                if (dst_down[j] > dst_down[j + 1]) {
                    swap(dst_down[j], dst_down[j + 1]);
                    swap(number[j], number[j + 1]);
                } else Loop++;
            }
        }
 
        for (int j = 0; j < number.size(); j++) { int i = number[j]; }
 
        number.erase(number.begin());
 
    } else // if run10_flag is false
    {
        for (int i = 0; i < Vec_layer.size(); i++) {
            if (Vec_layerid[i] != layerid) continue;
            number.push_back(i);
 
            dst.push_back(Min_dis(layerid, Vec_phi[i], Vec_v[i], l1));
        }
 
        int Loop(0);
 
        while (Loop != number.size() - 1) {
            Loop = 0;
            for (int j = 0; j < number.size() - 1; j++) {
                if (dst[j] > dst[j + 1]) {
                    swap(dst[j], dst[j + 1]);
                    swap(number[j], number[j + 1]);
                } else Loop++;
            }
        }
 
        number.erase(number.begin());
 
    } // end of if run10_flag
 
    for (int j = 0; j < number.size(); j++) { int i = number[j]; }
 
    sort(number.begin(), number.end(), less<int>());
    for (int i = number.size() - 1; i != -1; i--) //{ erase(number[i]); }//removed by Farinelli Riccardo in 2023.12.29
 
        return;
}

Referenced by Data_Closest(), and ToyMC().

FilterClusters()

void CgemLineFit::FilterClusters

(

)

Definition at line 4065 of file CgemLineFit.cxx.

                                 {
    int max_1_0(0), max_0_0(0);
    int _size = Vec_flag.size();
    max_1_0 = _size;
    max_0_0 = _size;
 
    for (int i = 0; i < Vec_flag.size(); i++) {
        if (1 == Vec_layerid[i] && 0 == Vec_sheetid[i]) {
            max_1_0 == i;
            break;
        }
    }
 
    for (int j = 0; j < Vec_flag.size(); j++) {
        if (0 == Vec_layerid[j] && 0 == Vec_sheetid[j]) {
            max_0_0 == j;
            break;
        }
    }
    for (int k = Vec_flag.size() - 1; k > max_0_0; k--) { erase(k); }
    for (int m = Vec_flag.size() - 2; m > max_1_0; m--) { erase(m); }
}

finalize()

StatusCode CgemLineFit::finalize

(

)

Definition at line 785 of file CgemLineFit.cxx.

                                 {
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in finalize()" << endreq;
    return StatusCode::SUCCESS;
}

fireStripID()

void CgemLineFit::fireStripID

(

)

Definition at line 3651 of file CgemLineFit.cxx.

                              {
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5) {
            for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_2_up_stripID[j] = Vec_XstripID[i][j]; }
            for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_2_up_stripID[j] = Vec_VstripID[i][j]; }
        }
        if (Vec_Virlayerid[i] == 4) {
            for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_2_down_stripID[j] = Vec_XstripID[i][j]; }
            for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_2_down_stripID[j] = Vec_VstripID[i][j]; }
        }
        if (Vec_Virlayerid[i] == 3) {
            for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_1_up_stripID[j] = Vec_XstripID[i][j]; }
            for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_1_up_stripID[j] = Vec_VstripID[i][j]; }
        }
        if (Vec_Virlayerid[i] == 2) {
            for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_1_down_stripID[j] = Vec_XstripID[i][j]; }
            for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_1_down_stripID[j] = Vec_VstripID[i][j]; }
        }
        if (Vec_Virlayerid[i] == 1) {
            for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_0_up_stripID[j] = Vec_XstripID[i][j]; }
            for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_0_up_stripID[j] = Vec_VstripID[i][j]; }
        }
        if (Vec_Virlayerid[i] == 0) {
            for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_0_down_stripID[j] = Vec_XstripID[i][j]; }
            for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_0_down_stripID[j] = Vec_VstripID[i][j]; }
        }
    }
}

Referenced by execute().

fireStripQ()

void CgemLineFit::fireStripQ

(

)

Definition at line 3557 of file CgemLineFit.cxx.

                             {
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5) {
            for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_2_up_stripQ[j] = Vec_XstripQ[i][j]; }
            for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_2_up_stripQ[j] = Vec_VstripQ[i][j]; }
        }
        if (Vec_Virlayerid[i] == 4) {
            for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_2_down_stripQ[j] = Vec_XstripQ[i][j]; }
            for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_2_down_stripQ[j] = Vec_VstripQ[i][j]; }
        }
        if (Vec_Virlayerid[i] == 3) {
            for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_1_up_stripQ[j] = Vec_XstripQ[i][j]; }
            for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_1_up_stripQ[j] = Vec_VstripQ[i][j]; }
        }
        if (Vec_Virlayerid[i] == 2) {
            for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_1_down_stripQ[j] = Vec_XstripQ[i][j]; }
            for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_1_down_stripQ[j] = Vec_VstripQ[i][j]; }
        }
        if (Vec_Virlayerid[i] == 1) {
            for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_0_up_stripQ[j] = Vec_XstripQ[i][j]; }
            for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_0_up_stripQ[j] = Vec_VstripQ[i][j]; }
        }
        if (Vec_Virlayerid[i] == 0) {
            for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_0_down_stripQ[j] = Vec_XstripQ[i][j]; }
            for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_0_down_stripQ[j] = Vec_VstripQ[i][j]; }
        }
    }
}

Referenced by execute().

fireStripT()

void CgemLineFit::fireStripT

(

)

Definition at line 3586 of file CgemLineFit.cxx.

                             {
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5) {
            for (int j = 0; j < Vec_XstripT[i].size(); j++) {
                x_2_up_stripT[j] = Vec_XstripT[i][j];
                x_2_up_stripTf[j] = Vec_XstripTf[i][j];
            }
            for (int j = 0; j < Vec_VstripT[i].size(); j++) {
                v_2_up_stripT[j] = Vec_VstripT[i][j];
                v_2_up_stripTf[j] = Vec_VstripTf[i][j];
            }
        }
        if (Vec_Virlayerid[i] == 4) {
            for (int j = 0; j < Vec_XstripT[i].size(); j++) {
                x_2_down_stripT[j] = Vec_XstripT[i][j];
                x_2_down_stripTf[j] = Vec_XstripTf[i][j];
            }
            for (int j = 0; j < Vec_VstripT[i].size(); j++) {
                v_2_down_stripT[j] = Vec_VstripT[i][j];
                v_2_down_stripTf[j] = Vec_VstripTf[i][j];
            }
        }
        if (Vec_Virlayerid[i] == 3) {
            for (int j = 0; j < Vec_XstripT[i].size(); j++) {
                x_1_up_stripT[j] = Vec_XstripT[i][j];
                x_1_up_stripTf[j] = Vec_XstripTf[i][j];
            }
            for (int j = 0; j < Vec_VstripT[i].size(); j++) {
                v_1_up_stripT[j] = Vec_VstripT[i][j];
                v_1_up_stripTf[j] = Vec_VstripTf[i][j];
            }
        }
        if (Vec_Virlayerid[i] == 2) {
            for (int j = 0; j < Vec_XstripT[i].size(); j++) {
                x_1_down_stripT[j] = Vec_XstripT[i][j];
                x_1_down_stripTf[j] = Vec_XstripTf[i][j];
            }
            for (int j = 0; j < Vec_VstripT[i].size(); j++) {
                v_1_down_stripT[j] = Vec_VstripT[i][j];
                v_1_down_stripTf[j] = Vec_VstripTf[i][j];
            }
        }
        if (Vec_Virlayerid[i] == 1) {
            for (int j = 0; j < Vec_XstripT[i].size(); j++) {
                x_0_up_stripT[j] = Vec_XstripT[i][j];
                x_0_up_stripTf[j] = Vec_XstripTf[i][j];
            }
            for (int j = 0; j < Vec_VstripT[i].size(); j++) {
                v_0_up_stripT[j] = Vec_VstripT[i][j];
                v_0_up_stripTf[j] = Vec_VstripTf[i][j];
            }
        }
        if (Vec_Virlayerid[i] == 0) {
            for (int j = 0; j < Vec_XstripT[i].size(); j++) {
                x_0_down_stripT[j] = Vec_XstripT[i][j];
                x_0_down_stripTf[j] = Vec_XstripTf[i][j];
            }
            for (int j = 0; j < Vec_VstripT[i].size(); j++) {
                v_0_down_stripT[j] = Vec_VstripT[i][j];
                v_0_down_stripTf[j] = Vec_VstripTf[i][j];
            }
        }
    }
}

Referenced by execute().

Fit()

TMinuit * CgemLineFit::Fit	(	StraightLine *	l_outer,
		int	sheetflag,
		int	typ )

Definition at line 2738 of file CgemLineFit.cxx.

                                                                         {
    // typ == 0 : 2D fit on XY plane and SZ plane indipendently, then use the fit results as initial values to perform  3D fit!
    // typ == 1 : 2D fit on XY plane
    // typ == 2 : 3D fit on XY and SZ plane
 
    double trkpar[4] = {-999, -999, -999, -999};
    double trkparErr[4] = {-999, -999, -999, -999};
 
    sheet_flag = ~_sheet_flag;
    dz_set = l_outer->dz();
    tanl_set = l_outer->tanl();
    Int_t ierflg, npari, nparx, istat1, istat2, istat3;
    Double_t fmin, edm, errdef;
    Double_t arglist[10];
 
    arglist[0] = 2000000;
    arglist[1] = 0.0001;
    if (typ == 2) {
        TMinuit *Min3 = new TMinuit(4);
        Min3->SetFCN(fcn3);
        Min3->SetPrintLevel(-1);
        Min3->mnparm(0, "dr", l_outer->dr(), 0.1, -1 * R_layer[2], R_layer[2], ierflg);
        Min3->mnparm(1, "phi0", l_outer->phi0(), 0.01, -100 * acos(-1), 100 * acos(-1),
                ierflg); // Why the range of phi0 is so big -- Guoaq--2020/11/01
        Min3->mnparm(2, "dz", l_outer->dz(), 0.1, -0.5 * length, 0.5 * length, ierflg);
        Min3->mnparm(3, "tanL", l_outer->tanl(), 0.01, -9999, 9999, ierflg);
        Min3->mnexcm("MIGRAD", arglist, 2, ierflg);
 
        return Min3;
    }
 
    if (typ == 0) {
 
        arglist[0] = 2000000;
        arglist[1] = 0.0001;
    }
    TMinuit *Min = new TMinuit(2);
    Min->SetPrintLevel(-1);
    Min->SetFCN(fcn);
    Min->SetErrorDef(1.0); // For Chisq
    Min->mnparm(0, "dr", l_outer->dr(), 0.001, -1 * R_layer[2], R_layer[2], ierflg);
    Min->mnparm(1, "phi0", l_outer->phi0(), 0.001, -100 * acos(-1), 100 * acos(-1), ierflg);
    arglist[0] = 2000000;
    arglist[1] = 0.001;
    Min->mnexcm("MIGRAD", arglist, 2, ierflg);
    Min->mnstat(fmin, edm, errdef, npari, nparx, istat1);
    Min->GetParameter(0, trkpar[0], trkparErr[0]);
    // cout <<"typ1 : trkparErr[0] ----> " <<trkparErr[0] << endl;
    Min->GetParameter(1, trkpar[1], trkparErr[1]);
 
    to_0_2pi(trkpar[1]);
    dr_set = trkpar[0];
    phi0_set = trkpar[1];
 
    // cout << "trkpar[0] = " << trkpar[0] << "  trkpar[1] = " << trkpar[1] << endl;
    chisq_1 = Min->fAmin;
    // cout << "chisq_1 : " << chisq_1 << endl;
    if (typ == 1) return Min;
    else delete Min;
 
    /////////////////////////Fit to z0 & tanlam////////////
    // arglist[0] = 2000000;
    // arglist[1] = 0.0001;
    TMinuit *Min2 = new TMinuit(2);
    Min2->SetFCN(fcn2);
    Min2->SetPrintLevel(-1);
    Min2->mnparm(0, "dz", l_outer->dz(), 0.001, -0.5 * length, 0.5 * length, ierflg);
    Min2->mnparm(1, "tanL", l_outer->tanl(), 0.001, -9999, 9999, ierflg);
    Min2->mnexcm("MIGRAD", arglist, 2, ierflg);
    Min2->mnstat(fmin, edm, errdef, npari, nparx, istat2);
    for (int i = 0; i < 2; i++) { Min2->GetParameter(i, trkpar[i + 2], trkparErr[i + 2]); }
    chisq_2 = Min2->fAmin;
 
    delete Min2;
    //////////////////re-Fit 4 parameters to get correct EDM /////////////
 
    // arglist[0] = 2000000;
    // arglist[1] = 0.0001;
    TMinuit *Min3 = new TMinuit(4);
    Min3->SetFCN(fcn3);
    Min3->SetPrintLevel(-1);
    Min3->mnparm(0, "dr", trkpar[0], 0.001, -1 * R_layer[2], R_layer[2], ierflg);
    Min3->mnparm(1, "phi0", trkpar[1], 0.001, -100 * acos(-1), 100 * acos(-1), ierflg);
    Min3->mnparm(2, "dz", trkpar[2], 0.001, -0.5 * length, 0.5 * length, ierflg);
    Min3->mnparm(3, "tanL", trkpar[3], 0.001, -9999, 9999, ierflg);
    Min3->mnexcm("MIGRAD", arglist, 2, ierflg);
    /////
    // Min3->mnstat(fmin, edm, errdef, npari, nparx, istat3);
    // for (int i = 0; i < 4; i++) { Min3->GetParameter(i, trkpar[i], trkparErr[i]); }
    // cout <<"typ0 : trkparErr[0] ----> " <<trkparErr[0] << endl;
    /////
 
    if (typ == 0) return Min3;
    else delete Min3;
}

Referenced by execute(), GetIntersection(), Loop_All(), and Loop_MaxQ().

Fit_Clusters()

void CgemLineFit::Fit_Clusters

(

double

par[]

)

Definition at line 3680 of file CgemLineFit.cxx.

                                           {
    vector<double> Vec_truth;
    Vec_truth.push_back(par[0]);
    Vec_truth.push_back(par[1]);
    Vec_truth.push_back(par[2]);
    Vec_truth.push_back(par[3]);
 
    vector<double> clusters = Get_Clusters(Vec_truth);
 
    x_0_up_f = clusters[0];
    v_0_up_f = clusters[1];
    x_0_down_f = clusters[2];
    v_0_down_f = clusters[3];
 
    x_1_up_f = clusters[4];
    v_1_up_f = clusters[5];
    x_1_down_f = clusters[6];
    v_1_down_f = clusters[7];
 
    x_2_up_f = clusters[8];
    v_2_up_f = clusters[9];
    x_2_down_f = clusters[10];
    v_2_down_f = clusters[11];
}

Referenced by execute().

func()

static double CgemLineFit::func ( double layer, double R, double dr, double phi0, double z0, double tglam, double V, double x )

static

Get4Par()

vector< double > CgemLineFit::Get4Par	(	HepLorentzVector	p4,
		Hep3Vector	bp )

Definition at line 3319 of file CgemLineFit.cxx.

                                                                      {
    vector<double> _V_par;
    _V_par.clear();
    double xb = p4.px();
    double yb = p4.py();
    double zb = p4.pz();
    double xa = bp[0];
    double ya = bp[1];
    double za = bp[2];
    double k = -1 * (xa * xb + ya * yb) / (xb * xb + yb * yb);
    double xc = k * xb + xa;
    double yc = k * yb + ya;
    double zc = k * zb + za;
    HepPoint3D a(xa, ya, za);
    HepPoint3D c(xc, yc, zc);
    StraightLine l(a, c);
    _V_par.push_back(l.dr());
    _V_par.push_back(l.phi0());
    _V_par.push_back(l.dz());
    _V_par.push_back(l.tanl());
    return _V_par;
}

Referenced by MC_truth().

Get_Clusters()

vector< double > CgemLineFit::Get_Clusters

(

vector< double >

Vec_truth

)

Definition at line 3893 of file CgemLineFit.cxx.

                                                                 {
    vector<double> clusters;
    StraightLine l(Vec_truth[0], Vec_truth[1], Vec_truth[2], Vec_truth[3]);
    double phiV_up[2], phiV_down[2];
    HepPoint3D Up, Down;
 
    double phiV_up_Tmp[2], phiV_down_Tmp[2];
    HepPoint3D Up_Tmp, Down_Tmp;
 
    if (!Align_Flag) Mp->getPointIdealGeom(0 * 2, l, Up, Down, phiV_up, phiV_down);
    else {
        Mp->getPointAligned(0, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(1, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    }
    to_0_2pi(phiV_up[0]);
    to_0_2pi(phiV_down[0]);
 
    clusters.push_back(phiV_up[0]);
    clusters.push_back(phiV_up[1]);
    clusters.push_back(phiV_down[0]);
    clusters.push_back(phiV_down[1]);
 
    if (!Align_Flag) Mp->getPointIdealGeom(1 * 2, l, Up, Down, phiV_up, phiV_down);
    else {
        Mp->getPointAligned(2, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(3, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    }
    to_0_2pi(phiV_up[0]);
    to_0_2pi(phiV_down[0]);
 
    clusters.push_back(phiV_up[0]);
    clusters.push_back(phiV_up[1]);
    clusters.push_back(phiV_down[0]);
    clusters.push_back(phiV_down[1]);
 
    if (!Align_Flag) Mp->getPointIdealGeom(2 * 2, l, Up, Down, phiV_up, phiV_down);
    else {
        Mp->getPointAligned(4, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(5, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    }
    to_0_2pi(phiV_up[0]);
    to_0_2pi(phiV_down[0]);
 
    clusters.push_back(phiV_up[0]);
    clusters.push_back(phiV_up[1]);
    clusters.push_back(phiV_down[0]);
    clusters.push_back(phiV_down[1]);
 
    return clusters;
}

Referenced by Fit_Clusters(), and Get_MCInf().

Get_MCInf()

bool CgemLineFit::Get_MCInf

(

)

Definition at line 3103 of file CgemLineFit.cxx.

                            {
    vector<double> Vec_truth;
    Vec_truth.clear();
    Vec_truth = MC_truth();
 
    vector<double> _Vec_truth = Trans(Vec_truth[0], Vec_truth[1], Vec_truth[2], Vec_truth[3]);
 
    MC_par0 = _Vec_truth[0];
    MC_par1 = _Vec_truth[1];
    MC_par2 = _Vec_truth[2];
    MC_par3 = _Vec_truth[3];
    MC_px = Vec_truth[4];
    MC_py = Vec_truth[5];
    MC_pz = Vec_truth[6];
 
    vector<double> V_MC_clusters = Get_Clusters(Vec_truth);
    x_0_up_mc = V_MC_clusters[0];
    v_0_up_mc = V_MC_clusters[1];
    x_0_down_mc = V_MC_clusters[2];
    v_0_down_mc = V_MC_clusters[3];
 
    x_1_up_mc = V_MC_clusters[4];
    v_1_up_mc = V_MC_clusters[5];
    x_1_down_mc = V_MC_clusters[6];
    v_1_down_mc = V_MC_clusters[7];
 
    x_2_up_mc = V_MC_clusters[8];
    v_2_up_mc = V_MC_clusters[9];
    x_2_down_mc = V_MC_clusters[10];
    v_2_down_mc = V_MC_clusters[11];
    return true;
}

Referenced by Loop_MaxQ(), and ToyMC().

Get_OtherIniPar()

void CgemLineFit::Get_OtherIniPar	(	int	clst_0,
		int	clst_1,
		int	clst_2 )

Definition at line 3367 of file CgemLineFit.cxx.

                                                                    {
 
    if (clst_2 == 0) { // only 2 layers situation
 
        if (clst_0 >= 2 && clst_1 >= 2) {
            StraightLine *l_ini_i = IniPar(Vec_phi[2], Vec_Z[2], 1, Vec_phi[3], Vec_Z[3],
                    0); // why use index of 4 and 5 because we only have 2 layer. I guess hongpeng want to get
            // the initial parameters by the innermost layer
 
            vector<double> inii_par = Trans(l_ini_i->dr(), l_ini_i->phi0(), l_ini_i->dz(), l_ini_i->tanl());
            inii_0 = inii_par[0];
            inii_1 = inii_par[1];
            inii_2 = inii_par[2];
            inii_3 = inii_par[3];
            delete l_ini_i;
        }
        if (clst_1 >= 2) {
            StraightLine *l_ini_m = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1],
                    2); // I guess Hongpeng want to get the initial parameters by the middle layer
            vector<double> inim_par = Trans(l_ini_m->dr(), l_ini_m->phi0(), l_ini_m->dz(), l_ini_m->tanl());
 
            inim_0 = inim_par[0];
            inim_1 = inim_par[1];
            inim_2 = inim_par[2];
            inim_3 = inim_par[3];
            delete l_ini_m;
        }
 
    } else {
        if (clst_0 >= 2 && clst_1 >= 2) {
            StraightLine *l_ini_i = IniPar(Vec_phi[4], Vec_Z[4], 1, Vec_phi[5], Vec_Z[5],
                    0); // why use index of 4 and 5 because we only have 2 layer. I guess hongpeng want to get
            // the initial parameters by the innermost layer
 
            vector<double> inii_par = Trans(l_ini_i->dr(), l_ini_i->phi0(), l_ini_i->dz(), l_ini_i->tanl());
            inii_0 = inii_par[0];
            inii_1 = inii_par[1];
            inii_2 = inii_par[2];
            inii_3 = inii_par[3];
            delete l_ini_i;
        }
        if (clst_1 >= 2) {
            StraightLine *l_ini_m = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3],
                    2); // I guess Hongpeng want to get the initial parameters by the middle layer
            vector<double> inim_par = Trans(l_ini_m->dr(), l_ini_m->phi0(), l_ini_m->dz(), l_ini_m->tanl());
 
            inim_0 = inim_par[0];
            inim_1 = inim_par[1];
            inim_2 = inim_par[2];
            inim_3 = inim_par[3];
            delete l_ini_m;
        }
    }
}

Referenced by execute().

GetIntersection()

void CgemLineFit::GetIntersection

(

StraightLine *

line

)

Definition at line 3705 of file CgemLineFit.cxx.

                                                    {
 
    if (debug) cout << "Start Getintersection" << endl;
    Double_t fmin, edm, errdef;
    Int_t ierflg, npari, nparx, istat;
    double trkpar[4] = {-999, -999, -999, -999};
    double trkparErr[4] = {-999, -999, -999, -999};
 
    double phiV_up[2], phiV_down[2];
    HepPoint3D Up, Down;
    double phiV_up_Tmp[2], phiV_down_Tmp[2];
    HepPoint3D Up_Tmp, Down_Tmp;
 
    double Ang[3];
 
    //StraightLine *l1a = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
    TMinuit *fit = Fit(line, fa - f21, 0);
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
    inter_1_flag = istat;
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    StraightLine *l1b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    Store_Trk(fit, 1000);
    delete fit;
 
    if (Align_FLAG) {
        Mp->getPointAligned(4, l1b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(5, l1b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    } else Mp->getPointIdealGeom(2 * 2, l1b, Up, Down, phiV_up, phiV_down);
 
    to_0_2pi(phiV_up[0]);
    to_0_2pi(phiV_down[0]);
    inter_1_x = phiV_up[0];
    inter_1_V = phiV_up[1];
    inter_1_z = Up[2];
 
    /*inter_4_x=phiV_down[0];
      inter_4_V=phiV_down[1];
      inter_4_z=Down[2];*/
 
    HepPoint3D l2_up(cos(inter_1_x) * R_layer[2], sin(inter_1_x) * R_layer[2], inter_1_z);
    // HepPoint3D l1_down(cos(Vec_phi[1])*R_layer[1],sin(Vec_phi[1])*R_layer[1],Vec_Z[1]);
 
    IncidentAngle(l1b, l2_up, pl_21->getStereoAngle(), Ang);
    ang_1_x = Ang[0];
    ang_1_v = Ang[1];
    ang_1 = Ang[2];
 
    /* IncidentAngle(l1b,l1_down,pl_11->getStereoAngle(),Ang);
       ang_4=Ang[0];
       ang_4_x=Ang[1];
       ang_4_v=Ang[2];*/
 
    //delete l1a;
    delete l1b;
 
    //StraightLine *l2a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
    fit = Fit(line, fa - f11, 0);
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
    inter_2_flag = istat;
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    Store_Trk(fit, 2000);
    delete fit;
    StraightLine *l2b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    if (Align_FLAG) {
        Mp->getPointAligned(2, l2b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(3, l2b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    } else Mp->getPointIdealGeom(1 * 2, l2b, Up, Down, phiV_up, phiV_down);
    to_0_2pi(phiV_up[0]);
    inter_2_x = phiV_up[0];
    inter_2_V = phiV_up[1];
    inter_2_z = Up[2];
    HepPoint3D l1_up(cos(inter_2_x) * R_layer[1], sin(inter_2_x) * R_layer[1], inter_2_z);
    IncidentAngle(l2b, l1_up, pl_11->getStereoAngle(), Ang);
    ang_2_x = Ang[0];
    ang_2_v = Ang[1];
    ang_2 = Ang[2];
    //delete l2a;
    delete l2b;
 
    // if(Run10_Flag)
    // {
    //StraightLine *l3a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
    fit = Fit(line, fa - f01, 0);
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
    inter_3_flag = istat;
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    Store_Trk(fit, 3000);
    delete fit;
    StraightLine *l3b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
 
    if (Align_FLAG) {
        Mp->getPointAligned(0, l3b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(1, l3b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    } else Mp->getPointIdealGeom(0 * 2, l3b, Up, Down, phiV_up, phiV_down);
    to_0_2pi(phiV_up[0]);
    inter_3_x = phiV_up[0];
    inter_3_V = phiV_up[1];
    inter_3_z = Up[2];
 
    HepPoint3D l0_up(cos(inter_3_x) * R_layer[0], sin(inter_3_x) * R_layer[0], inter_3_z);
    IncidentAngle(l3b, l0_up, pl_00->getStereoAngle(), Ang);
    ang_3_x = Ang[0];
    ang_3_v = Ang[1];
    ang_3 = Ang[2];
    //delete l3a;
    delete l3b;
    //}
 
    //StraightLine *l4a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
    fit = Fit(line, fa - f00, 0);
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
    inter_4_flag = istat;
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    Store_Trk(fit, 4000);
    delete fit;
    StraightLine *l4b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    if (Align_FLAG) {
        Mp->getPointAligned(0, l4b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(1, l4b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    } else Mp->getPointIdealGeom(0 * 2, l4b, Up, Down, phiV_up, phiV_down);
 
    // to_0_2pi(phiV_up[0]);
    to_0_2pi(phiV_down[0]);
    inter_4_x = phiV_down[0];
    inter_4_V = phiV_down[1];
    inter_4_z = Down[2];
 
    HepPoint3D l0_down(cos(inter_4_x) * R_layer[0], sin(inter_4_x) * R_layer[0], inter_4_z);
    IncidentAngle(l4b, l0_down, pl_00->getStereoAngle(), Ang);
    ang_4_x = Ang[0];
    ang_4_v = Ang[1];
    ang_4 = Ang[2];
    //delete l4a;
    delete l4b;
 
    //StraightLine *l5a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
    fit = Fit(line, fa - f10, 0);
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
    inter_5_flag = istat;
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    Store_Trk(fit, 5000);
    delete fit;
    StraightLine *l5b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    if (Align_FLAG) {
        Mp->getPointAligned(2, l5b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(3, l5b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    } else Mp->getPointIdealGeom(1 * 2, l5b, Up, Down, phiV_up, phiV_down);
    to_0_2pi(phiV_down[0]);
    inter_5_x = phiV_down[0];
    inter_5_V = phiV_down[1];
    inter_5_z = Down[2];
    HepPoint3D l1_down(cos(inter_5_x) * R_layer[1], sin(inter_5_x) * R_layer[1], inter_5_z);
    IncidentAngle(l5b, l1_down, pl_10->getStereoAngle(), Ang);
    ang_5_x = Ang[0];
    ang_5_v = Ang[1];
    ang_5 = Ang[2];
    //delete l5a;
    delete l5b;
    //}
 
    //StraightLine *l6a = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
    fit = Fit(line, fa - f20, 0);
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
    inter_6_flag = istat;
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    Store_Trk(fit, 6000);
    delete fit;
    StraightLine *l6b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
    if (Align_FLAG) {
        Mp->getPointAligned(4, l6b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
        Mp->getPointAligned(5, l6b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
    } else Mp->getPointIdealGeom(2 * 2, l6b, Up, Down, phiV_up, phiV_down);
 
// to_0_2pi(phiV_up[0]);
to_0_2pi(phiV_down[0]);
inter_6_x = phiV_down[0];
inter_6_V = phiV_down[1];
inter_6_z = Down[2];
 
HepPoint3D l2_down(cos(inter_6_x) * R_layer[2], sin(inter_6_x) * R_layer[2], inter_6_z);
IncidentAngle(l6b, l2_down, pl_20->getStereoAngle(), Ang);
ang_6_x = Ang[0];
ang_6_v = Ang[1];
ang_6 = Ang[2];
//delete l6a;
delete l6b;
}

Referenced by execute().

GetNMaxQ()

vector< int > CgemLineFit::GetNMaxQ ( vector< double > Q_11, vector< int > L_11, int Nmax = 3 )

static

Definition at line 4036 of file CgemLineFit.cxx.

                                                                                     {
    // chose the N X clusters with the largest charge  and save corresponding index
    vector<int> L_11_s;
    vector<double> Q_11_s;
    double t, tt;
    if (L_11.size() > Nmax) {
        for (int i = 0; i < Nmax; i++) {
            for (int j = i + 1; j < L_11.size(); j++) {
                if (Q_11[i] < Q_11[j]) {
                    t = Q_11[i];
                    Q_11[i] = Q_11[j];
                    Q_11[j] = t;
 
                    tt = L_11[i];
                    L_11[i] = L_11[j];
                    L_11[j] = tt;
                }
            }
            Q_11_s.push_back(Q_11[i]);
            L_11_s.push_back(L_11[i]);
        }
 
    } else {
        L_11_s.assign(L_11.begin(), L_11.end());
        Q_11_s.assign(Q_11.begin(), Q_11.end());
    }
    return L_11_s;
}

Referenced by Loop_MaxQ().

GetVirLay()

int CgemLineFit::GetVirLay	(	int	geolay,
		double	phi )

Definition at line 4030 of file CgemLineFit.cxx.

                                                 {
    int virsheet = (phi < 0) ? 0 : 1;
    int virlay = geolay * 2 + virsheet;
    return virlay;
}

Referenced by Data_Closest(), Data_Max(), Loop_All(), Loop_MaxQ(), Min_dis(), Min_dis_down(), Min_dis_up(), and ToyMC().

HitPos()

bool CgemLineFit::HitPos	(	HepLorentzVector	p4,
		Hep3Vector	bp )

Definition at line 3342 of file CgemLineFit.cxx.

                                                           {
    double L(500), W(160), H(600);
 
    vector<double> _V_par;
    _V_par.clear();
    double xb = p4.px();
    double yb = p4.py();
    double zb = p4.pz();
    double xa = bp[0];
    double ya = bp[1];
    double za = bp[2];
    double k = H / yb * -1;
    double xc = k * xb + xa;
    double yc = k * yb + ya;
    double zc = k * zb + za;
    hit_x = xc;
    hit_y = yc;
    hit_z = zc;
 
    if (fabs(xc) > W / 2) return false;
    if (fabs(zc) > L / 2) return false;
 
    else return true;
}

Referenced by MC_truth().

InAngle()

void CgemLineFit::InAngle

(

StraightLine

sl

)

IncidentAngle()

void CgemLineFit::IncidentAngle ( StraightLine * cosmic_ray, HepPoint3D cluster, double theta, double ang[] )

static

Definition at line 4149 of file CgemLineFit.cxx.

                                                                                                        {
 
    const Hep3Vector Normals(cluster[0], cluster[1], 0);
 
    HepPoint3D x1 = cosmic_ray->x(1);
    HepPoint3D x2 = cosmic_ray->x(-1);
    const Hep3Vector CosmicRay(x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2]);
    const Hep3Vector z(0, 0, 1);
    const Hep3Vector phi = z.cross(Normals);
    Hep3Vector _phi = phi;
 
    const Hep3Vector V = _phi.rotate(-1 * theta, Normals);
 
    const Hep3Vector N_V = Normals.cross(V);
    const Hep3Vector N_phi = Normals.cross(phi);
 
    Hep3Vector C_phi = CosmicRay - CosmicRay.project(N_phi);
    double s_C_phi = Normals.dot(C_phi);
    Hep3Vector v_phi = Normals.cross(C_phi);
    // ang[0]=Normals.angle(C_phi);
    if (s_C_phi <= 0) {
        double s_phi = v_phi.dot(z);
        if (s_phi >= 0)
            ang[0] = -(Normals.angle(-C_phi)); // 宇宙线在x-y平面上的投影和法线的夹角,空间角不分正负，所以加几个判断条件区分正负方向
        if (s_phi < 0) ang[0] = Normals.angle(-C_phi);
        /* if(s_phi==0)
           {
           ang[0]=Normals.angle(C_phi);
           cout<<"s_C_phi<=0  s_phi==0     "<<ang[0]<<endl;
           }
           if(ang[0]==0)
           cout<<"s_C_phi<=0      "<<ang[0]<<endl;*/
    } else if (s_C_phi > 0) {
        double s_phi = v_phi.dot(z);
        if (s_phi <= 0)
            ang[0] = -(Normals.angle(
                        C_phi)); // 宇宙线在x-y平面上的投影和法线的夹角,C_phi,Normals都是向量,算出的间角不分正负，所以加几个判断条件区分正负方向
        if (s_phi > 0) ang[0] = Normals.angle(C_phi);
        /* if(s_phi==0)
           {
           ang[0]=Normals.angle(C_phi);
           cout<<"s_C_phi>0  s_phi==0     "<<ang[0]<<endl;
           }
           if(ang[0]==0)
           cout<<"s_C_phi>0      "<<ang[0]<<endl;*/
    }
    Hep3Vector C_V = CosmicRay - CosmicRay.project(N_V); // C_V是宇宙线在与V向垂直的平面上的投影
    double s_C_V = Normals.dot(C_V);
    Hep3Vector v_V = Normals.cross(C_V);
    // ang[1]=Normals.angle(C_V);
    if (s_C_V <= 0) {
        double s_V = v_V.dot(z);
        if (s_V >= 0) ang[1] = -(Normals.angle(-C_V)); // 宇宙线在与V向垂直的平面上的投影和法线的夹角
        if (s_V < 0) ang[1] = Normals.angle(-C_V);
        /* if(s_V==0)
           {
           ang[1]=Normals.angle(C_V);
           cout<<"s_C_V<=0  s_V==0     "<<ang[1]<<endl;
           }
           if(ang[1]==0)
           cout<<"s_C_V<=0     "<<ang[1]<<endl;*/
    } else if (s_C_V > 0) {
        double s_V = v_V.dot(z);
        if (s_V <= 0) ang[1] = -(Normals.angle(C_V));
        if (s_V > 0) ang[1] = Normals.angle(C_V);
        /* if(s_V==0)
           {
           ang[1]=Normals.angle(C_V);
           cout<<"s_C_V>0  s_V==0     "<<ang[1]<<endl;
           }
           if(ang[1]==0)
           cout<<"s_C_V>0     "<<ang[1]<<endl;*/
    }
    ang[2] = Normals.angle(CosmicRay); // 宇宙线和法线的夹角即是入射角,但这是个空间上的三维的角，所以我们要用投影平面的夹角。
    /*Hep3Vector C_V=CosmicRay-CosmicRay.project(N_V);
      ang[1]=Normals.angle(C_V);
      ang[2]=Normals.angle(CosmicRay);*/
}

Referenced by CalSigma2(), and GetIntersection().

IniPar()

StraightLine * CgemLineFit::IniPar	(	double	phi1,
		double	z1,
		int	i,
		double	phi2,
		double	z2,
		int	j )

Definition at line 2991 of file CgemLineFit.cxx.

                                                                                              {
    int Gi = int(i / 2);
    int Gj = int(j / 2);
 
    HepPoint3D up(R_layer[Gi] * cos(phi1), R_layer[Gi] * sin(phi1), z1);
    HepPoint3D down(R_layer[Gj] * cos(phi2), R_layer[Gj] * sin(phi2), z2);
 
    if (Align_Flag) up = Al->point_invTransform(i, up); // added by Guoaq-2020/03/13
    if (Align_Flag) down = Al->point_invTransform(j, down);
 
    StraightLine *l = new StraightLine(up, down);
    return l;
}

Referenced by Data_Closest(), Data_Max(), Get_OtherIniPar(), Loop_All(), Loop_MaxQ(), and ToyMC().

initialize()

StatusCode CgemLineFit::initialize

(

)

Definition at line 147 of file CgemLineFit.cxx.

                                   {
    MsgStream log(msgSvc(), name());
    sigma2 = sigma * sigma;
    Align_FLAG = Align_Flag;
    log << MSG::INFO << "in initialize()" << endreq;
    StatusCode sc;
    NTuplePtr cosmic(ntupleSvc(), "LineFit/cosmic");
    if (cosmic) m_tuple = cosmic;
    else {
        m_tuple = ntupleSvc()->book("LineFit/cosmic", CLID_ColumnWiseTuple, "N-Tuple example");
        if (m_tuple) {
 
            sc = m_tuple->addItem("run", run);
            sc = m_tuple->addItem("event", event);
            sc = m_tuple->addItem("R0", R0);
            sc = m_tuple->addItem("R1", R1);
            sc = m_tuple->addItem("R2", R2);
            sc = m_tuple->addItem("clst_0", clst_0);
            sc = m_tuple->addItem("clst_1", clst_1);
            sc = m_tuple->addItem("clst_2", clst_2);
            sc = m_tuple->addItem("clst_01", clst_01);
            sc = m_tuple->addItem("clst_00", clst_00);
            sc = m_tuple->addItem("clst_11", clst_11);
            sc = m_tuple->addItem("clst_10", clst_10);
            sc = m_tuple->addItem("clst_21", clst_21);
            sc = m_tuple->addItem("clst_20", clst_20);
            sc = m_tuple->addItem("status1", status1);
            sc = m_tuple->addItem("status2", status2);
            sc = m_tuple->addItem("status3", status3);
            sc = m_tuple->addItem("ini_0", ini_0);
            sc = m_tuple->addItem("ini_1", ini_1);
            sc = m_tuple->addItem("ini_2", ini_2);
            sc = m_tuple->addItem("ini_3", ini_3);
            sc = m_tuple->addItem("inim_0", inim_0);
            sc = m_tuple->addItem("inim_1", inim_1);
            sc = m_tuple->addItem("inim_2", inim_2);
            sc = m_tuple->addItem("inim_3", inim_3);
            sc = m_tuple->addItem("inii_0", inii_0);
            sc = m_tuple->addItem("inii_1", inii_1);
            sc = m_tuple->addItem("inii_2", inii_2);
            sc = m_tuple->addItem("inii_3", inii_3);
            sc = m_tuple->addItem("par0", par0);
            sc = m_tuple->addItem("par1", par1);
            sc = m_tuple->addItem("par2", par2);
            sc = m_tuple->addItem("par3", par3);
            sc = m_tuple->addItem("MC_par0", MC_par0);
            sc = m_tuple->addItem("MC_par1", MC_par1);
            sc = m_tuple->addItem("MC_par2", MC_par2);
            sc = m_tuple->addItem("MC_par3", MC_par3);
            sc = m_tuple->addItem("MC_px", MC_px);
            sc = m_tuple->addItem("MC_py", MC_py);
            sc = m_tuple->addItem("MC_pz", MC_pz);
            sc = m_tuple->addItem("Err_par0", Err_par0);
            sc = m_tuple->addItem("Err_par1", Err_par1);
            sc = m_tuple->addItem("Err_par2", Err_par2);
            sc = m_tuple->addItem("Err_par3", Err_par3);
 
            sc = m_tuple->addItem("x_2_up", x_2_up);
            sc = m_tuple->addItem("x_1_up", x_1_up);
            sc = m_tuple->addItem("x_0_up", x_0_up);
            sc = m_tuple->addItem("v_2_up", v_2_up);
            sc = m_tuple->addItem("v_1_up", v_1_up);
            sc = m_tuple->addItem("v_0_up", v_0_up);
            sc = m_tuple->addItem("z_2_up", z_2_up);
            sc = m_tuple->addItem("z_1_up", z_1_up);
            sc = m_tuple->addItem("z_0_up", z_0_up);
            sc = m_tuple->addItem("x_2_down", x_2_down);
            sc = m_tuple->addItem("x_1_down", x_1_down);
            sc = m_tuple->addItem("x_0_down", x_0_down);
            sc = m_tuple->addItem("v_2_down", v_2_down);
            sc = m_tuple->addItem("v_1_down", v_1_down);
            sc = m_tuple->addItem("v_0_down", v_0_down);
            sc = m_tuple->addItem("z_2_down", z_2_down);
            sc = m_tuple->addItem("z_1_down", z_1_down);
            sc = m_tuple->addItem("z_0_down", z_0_down);
 
            sc = m_tuple->addItem("x_2_up_m", x_2_up_m);
            sc = m_tuple->addItem("x_1_up_m", x_1_up_m);
            sc = m_tuple->addItem("x_0_up_m", x_0_up_m);
            sc = m_tuple->addItem("v_2_up_m", v_2_up_m);
            sc = m_tuple->addItem("v_1_up_m", v_1_up_m);
            sc = m_tuple->addItem("v_0_up_m", v_0_up_m);
            sc = m_tuple->addItem("z_2_up_m", z_2_up_m);
            sc = m_tuple->addItem("z_1_up_m", z_1_up_m);
            sc = m_tuple->addItem("z_0_up_m", z_0_up_m);
 
            sc = m_tuple->addItem("x_2_down_m", x_2_down_m);
            sc = m_tuple->addItem("x_1_down_m", x_1_down_m);
            sc = m_tuple->addItem("x_0_down_m", x_0_down_m);
            sc = m_tuple->addItem("v_2_down_m", v_2_down_m);
            sc = m_tuple->addItem("v_1_down_m", v_1_down_m);
            sc = m_tuple->addItem("v_0_down_m", v_0_down_m);
            sc = m_tuple->addItem("z_2_down_m", z_2_down_m);
            sc = m_tuple->addItem("z_1_down_m", z_1_down_m);
            sc = m_tuple->addItem("z_0_down_m", z_0_down_m);
 
            sc = m_tuple->addItem("x_2_up_f", x_2_up_f);
            sc = m_tuple->addItem("x_1_up_f", x_1_up_f);
            sc = m_tuple->addItem("x_0_up_f", x_0_up_f);
            sc = m_tuple->addItem("v_2_up_f", v_2_up_f);
            sc = m_tuple->addItem("v_1_up_f", v_1_up_f);
            sc = m_tuple->addItem("v_0_up_f", v_0_up_f);
 
            sc = m_tuple->addItem("x_2_down_f", x_2_down_f);
            sc = m_tuple->addItem("x_1_down_f", x_1_down_f);
            sc = m_tuple->addItem("x_0_down_f", x_0_down_f);
            sc = m_tuple->addItem("v_2_down_f", v_2_down_f);
            sc = m_tuple->addItem("v_1_down_f", v_1_down_f);
            sc = m_tuple->addItem("v_0_down_f", v_0_down_f);
 
            sc = m_tuple->addItem("x_2_up_pre", x_2_up_mc);
            sc = m_tuple->addItem("x_1_up_pre", x_1_up_mc);
            sc = m_tuple->addItem("x_0_up_pre", x_0_up_mc);
            sc = m_tuple->addItem("v_2_up_pre", v_2_up_mc);
            sc = m_tuple->addItem("v_1_up_pre", v_1_up_mc);
            sc = m_tuple->addItem("v_0_up_pre", v_0_up_mc);
 
            sc = m_tuple->addItem("x_2_down_pre", x_2_down_mc);
            sc = m_tuple->addItem("x_1_down_pre", x_1_down_mc);
            sc = m_tuple->addItem("x_0_down_pre", x_0_down_mc);
            sc = m_tuple->addItem("v_2_down_pre", v_2_down_mc);
            sc = m_tuple->addItem("v_1_down_pre", v_1_down_mc);
            sc = m_tuple->addItem("v_0_down_pre", v_0_down_mc);
 
            sc = m_tuple->addItem("x_2_up_hit", x_2_up_mc2);
            sc = m_tuple->addItem("x_1_up_hit", x_1_up_mc2);
            sc = m_tuple->addItem("x_0_up_hit", x_0_up_mc2);
            sc = m_tuple->addItem("v_2_up_hit", v_2_up_mc2);
            sc = m_tuple->addItem("v_1_up_hit", v_1_up_mc2);
            sc = m_tuple->addItem("v_0_up_hit", v_0_up_mc2);
            sc = m_tuple->addItem("z_2_up_hit", z_2_up_mc2);
            sc = m_tuple->addItem("z_1_up_hit", z_1_up_mc2);
            sc = m_tuple->addItem("z_0_up_hit", z_0_up_mc2);
 
            sc = m_tuple->addItem("x_2_down_hit", x_2_down_mc2);
            sc = m_tuple->addItem("x_1_down_hit", x_1_down_mc2);
            sc = m_tuple->addItem("x_0_down_hit", x_0_down_mc2);
            sc = m_tuple->addItem("v_2_down_hit", v_2_down_mc2);
            sc = m_tuple->addItem("v_1_down_hit", v_1_down_mc2);
            sc = m_tuple->addItem("v_0_down_hit", v_0_down_mc2);
            sc = m_tuple->addItem("z_2_down_hit", z_2_down_mc2);
            sc = m_tuple->addItem("z_1_down_hit", z_1_down_mc2);
            sc = m_tuple->addItem("z_0_down_hit", z_0_down_mc2);
 
            sc = m_tuple->addItem("x_2_up_size", x_2_up_size, 0, 100);
            sc = m_tuple->addItem("x_1_up_size", x_1_up_size, 0, 100);
            sc = m_tuple->addItem("x_0_up_size", x_0_up_size, 0, 100);
            sc = m_tuple->addItem("v_2_up_size", v_2_up_size, 0, 100);
            sc = m_tuple->addItem("v_1_up_size", v_1_up_size, 0, 100);
            sc = m_tuple->addItem("v_0_up_size", v_0_up_size, 0, 100);
 
            sc = m_tuple->addItem("x_2_down_size", x_2_down_size, 0, 100);
            sc = m_tuple->addItem("x_1_down_size", x_1_down_size, 0, 100);
            sc = m_tuple->addItem("x_0_down_size", x_0_down_size, 0, 100);
            sc = m_tuple->addItem("v_2_down_size", v_2_down_size, 0, 100);
            sc = m_tuple->addItem("v_1_down_size", v_1_down_size, 0, 100);
            sc = m_tuple->addItem("v_0_down_size", v_0_down_size, 0, 100);
 
            sc = m_tuple->addItem("x_2_up_Q", x_2_up_Q);
            sc = m_tuple->addItem("x_1_up_Q", x_1_up_Q);
            sc = m_tuple->addItem("x_0_up_Q", x_0_up_Q);
            sc = m_tuple->addItem("v_2_up_Q", v_2_up_Q);
            sc = m_tuple->addItem("v_1_up_Q", v_1_up_Q);
            sc = m_tuple->addItem("v_0_up_Q", v_0_up_Q);
 
            sc = m_tuple->addItem("x_2_down_Q", x_2_down_Q);
            sc = m_tuple->addItem("x_1_down_Q", x_1_down_Q);
            sc = m_tuple->addItem("x_0_down_Q", x_0_down_Q);
            sc = m_tuple->addItem("v_2_down_Q", v_2_down_Q);
            sc = m_tuple->addItem("v_1_down_Q", v_1_down_Q);
            sc = m_tuple->addItem("v_0_down_Q", v_0_down_Q);
 
            sc = m_tuple->addIndexedItem("x_2_up_stripQ", x_2_up_size, x_2_up_stripQ);
            sc = m_tuple->addIndexedItem("x_1_up_stripQ", x_1_up_size, x_1_up_stripQ);
            sc = m_tuple->addIndexedItem("x_0_up_stripQ", x_0_up_size, x_0_up_stripQ);
            sc = m_tuple->addIndexedItem("v_2_up_stripQ", v_2_up_size, v_2_up_stripQ);
            sc = m_tuple->addIndexedItem("v_1_up_stripQ", v_1_up_size, v_1_up_stripQ);
            sc = m_tuple->addIndexedItem("v_0_up_stripQ", v_0_up_size, v_0_up_stripQ);
 
            sc = m_tuple->addIndexedItem("x_2_down_stripQ", x_2_down_size, x_2_down_stripQ);
            sc = m_tuple->addIndexedItem("x_1_down_stripQ", x_1_down_size, x_1_down_stripQ);
            sc = m_tuple->addIndexedItem("x_0_down_stripQ", x_0_down_size, x_0_down_stripQ);
            sc = m_tuple->addIndexedItem("v_2_down_stripQ", v_2_down_size, v_2_down_stripQ);
            sc = m_tuple->addIndexedItem("v_1_down_stripQ", v_1_down_size, v_1_down_stripQ);
            sc = m_tuple->addIndexedItem("v_0_down_stripQ", v_0_down_size, v_0_down_stripQ);
 
            sc = m_tuple->addIndexedItem("x_2_up_stripT", x_2_up_size, x_2_up_stripT);
            sc = m_tuple->addIndexedItem("x_1_up_stripT", x_1_up_size, x_1_up_stripT);
            sc = m_tuple->addIndexedItem("x_0_up_stripT", x_0_up_size, x_0_up_stripT);
            sc = m_tuple->addIndexedItem("v_2_up_stripT", v_2_up_size, v_2_up_stripT);
            sc = m_tuple->addIndexedItem("v_1_up_stripT", v_1_up_size, v_1_up_stripT);
            sc = m_tuple->addIndexedItem("v_0_up_stripT", v_0_up_size, v_0_up_stripT);
 
            sc = m_tuple->addIndexedItem("x_2_down_stripT", x_2_down_size, x_2_down_stripT);
            sc = m_tuple->addIndexedItem("x_1_down_stripT", x_1_down_size, x_1_down_stripT);
            sc = m_tuple->addIndexedItem("x_0_down_stripT", x_0_down_size, x_0_down_stripT);
            sc = m_tuple->addIndexedItem("v_2_down_stripT", v_2_down_size, v_2_down_stripT);
            sc = m_tuple->addIndexedItem("v_1_down_stripT", v_1_down_size, v_1_down_stripT);
            sc = m_tuple->addIndexedItem("v_0_down_stripT", v_0_down_size, v_0_down_stripT);
 
            sc = m_tuple->addIndexedItem("x_2_up_stripTf", x_2_up_size, x_2_up_stripTf);
            sc = m_tuple->addIndexedItem("x_1_up_stripTf", x_1_up_size, x_1_up_stripTf);
            sc = m_tuple->addIndexedItem("x_0_up_stripTf", x_0_up_size, x_0_up_stripTf);
            sc = m_tuple->addIndexedItem("v_2_up_stripTf", v_2_up_size, v_2_up_stripTf);
            sc = m_tuple->addIndexedItem("v_1_up_stripTf", v_1_up_size, v_1_up_stripTf);
            sc = m_tuple->addIndexedItem("v_0_up_stripTf", v_0_up_size, v_0_up_stripTf);
 
            sc = m_tuple->addIndexedItem("x_2_down_stripTf", x_2_down_size, x_2_down_stripTf);
            sc = m_tuple->addIndexedItem("x_1_down_stripTf", x_1_down_size, x_1_down_stripTf);
            sc = m_tuple->addIndexedItem("x_0_down_stripTf", x_0_down_size, x_0_down_stripTf);
            sc = m_tuple->addIndexedItem("v_2_down_stripTf", v_2_down_size, v_2_down_stripTf);
            sc = m_tuple->addIndexedItem("v_1_down_stripTf", v_1_down_size, v_1_down_stripTf);
            sc = m_tuple->addIndexedItem("v_0_down_stripTf", v_0_down_size, v_0_down_stripTf);
 
            sc = m_tuple->addIndexedItem("x_2_up_stripID", x_2_up_size, x_2_up_stripID);
            sc = m_tuple->addIndexedItem("x_1_up_stripID", x_1_up_size, x_1_up_stripID);
            sc = m_tuple->addIndexedItem("x_0_up_stripID", x_0_up_size, x_0_up_stripID);
            sc = m_tuple->addIndexedItem("v_2_up_stripID", v_2_up_size, v_2_up_stripID);
            sc = m_tuple->addIndexedItem("v_1_up_stripID", v_1_up_size, v_1_up_stripID);
            sc = m_tuple->addIndexedItem("v_0_up_stripID", v_0_up_size, v_0_up_stripID);
 
            sc = m_tuple->addIndexedItem("x_2_down_stripID", x_2_down_size, x_2_down_stripID);
            sc = m_tuple->addIndexedItem("x_1_down_stripID", x_1_down_size, x_1_down_stripID);
            sc = m_tuple->addIndexedItem("x_0_down_stripID", x_0_down_size, x_0_down_stripID);
            sc = m_tuple->addIndexedItem("v_2_down_stripID", v_2_down_size, v_2_down_stripID);
            sc = m_tuple->addIndexedItem("v_1_down_stripID", v_1_down_size, v_1_down_stripID);
            sc = m_tuple->addIndexedItem("v_0_down_stripID", v_0_down_size, v_0_down_stripID);
 
            sc = m_tuple->addItem("inter_1_x", inter_1_x);
            sc = m_tuple->addItem("inter_2_x", inter_2_x);
            sc = m_tuple->addItem("inter_3_x", inter_3_x);
            sc = m_tuple->addItem("inter_4_x", inter_4_x);
            sc = m_tuple->addItem("inter_5_x", inter_5_x);
            sc = m_tuple->addItem("inter_6_x", inter_6_x);
            sc = m_tuple->addItem("inter_x", inter_x);
 
            sc = m_tuple->addItem("inter_1_V", inter_1_V);
            sc = m_tuple->addItem("inter_2_V", inter_2_V);
            sc = m_tuple->addItem("inter_3_V", inter_3_V);
            sc = m_tuple->addItem("inter_4_V", inter_4_V);
            sc = m_tuple->addItem("inter_5_V", inter_5_V);
            sc = m_tuple->addItem("inter_6_V", inter_6_V);
            sc = m_tuple->addItem("inter_V", inter_V);
 
            sc = m_tuple->addItem("inter_1_flag", inter_1_flag);
            sc = m_tuple->addItem("inter_2_flag", inter_2_flag);
            sc = m_tuple->addItem("inter_3_flag", inter_3_flag);
            sc = m_tuple->addItem("inter_4_flag", inter_4_flag);
            sc = m_tuple->addItem("inter_5_flag", inter_5_flag);
            sc = m_tuple->addItem("inter_6_flag", inter_6_flag);
            sc = m_tuple->addItem("inter_flag", inter_flag);
 
            sc = m_tuple->addItem("inter_1_z", inter_1_z);
            sc = m_tuple->addItem("inter_2_z", inter_2_z);
            sc = m_tuple->addItem("inter_3_z", inter_3_z);
            sc = m_tuple->addItem("inter_4_z", inter_4_z);
            sc = m_tuple->addItem("inter_5_z", inter_5_z);
            sc = m_tuple->addItem("inter_6_z", inter_6_z);
            sc = m_tuple->addItem("inter_z", inter_z);
 
            sc = m_tuple->addItem("chisq_1", chisq_1);
            sc = m_tuple->addItem("chisq_2", chisq_2);
            sc = m_tuple->addItem("chisq_3", chisq_3);
            sc = m_tuple->addItem("pos_x", pos_x);
            sc = m_tuple->addItem("pos_y", pos_y);
            sc = m_tuple->addItem("pos_z", pos_z);
 
            sc = m_tuple->addItem("hit_x", hit_x);
            sc = m_tuple->addItem("hit_y", hit_y);
            sc = m_tuple->addItem("hit_z", hit_z);
 
            sc = m_tuple->addItem("ang_1", ang_1);
            sc = m_tuple->addItem("x_2_up_ang", ang_1_x);
            sc = m_tuple->addItem("v_2_up_ang", ang_1_v);
 
            sc = m_tuple->addItem("ang_2", ang_2);
            sc = m_tuple->addItem("x_1_up_ang", ang_2_x);
            sc = m_tuple->addItem("v_1_up_ang", ang_2_v);
 
            sc = m_tuple->addItem("ang_3", ang_3);
            sc = m_tuple->addItem("x_0_up_ang", ang_3_x);
            sc = m_tuple->addItem("v_0_up_ang", ang_3_v);
 
            sc = m_tuple->addItem("ang_4", ang_4);
            sc = m_tuple->addItem("x_0_down_ang", ang_4_x);
            sc = m_tuple->addItem("v_0_down_ang", ang_4_v);
 
            sc = m_tuple->addItem("ang_5", ang_5);
            sc = m_tuple->addItem("x_1_down_ang", ang_5_x);
            sc = m_tuple->addItem("v_1_down_ang", ang_5_v);
 
            sc = m_tuple->addItem("ang_6", ang_6);
            sc = m_tuple->addItem("x_2_down_ang", ang_6_x);
            sc = m_tuple->addItem("v_2_down_ang", ang_6_v);
 
            sc = m_tuple->addItem("Angle", Angle);
            sc = m_tuple->addItem("ang_x", ang_x);
            sc = m_tuple->addItem("ang_v", ang_v);
 
            sc = m_tuple->addItem("Rec_phi", Rec_phi);
            sc = m_tuple->addItem("Rec_V", Rec_V);
            sc = m_tuple->addItem("Rec_Z", Rec_Z);
            sc = m_tuple->addItem("Test_phi", Test_phi);
            sc = m_tuple->addItem("Test_V", Test_V);
            sc = m_tuple->addItem("Test_Z", Test_Z);
        } else {
            log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple) << endmsg;
            return StatusCode::FAILURE;
        }
    }
 
    NTuplePtr Track(ntupleSvc(), "LineFit/Track");
    if (Track) m_tuple_track = Track;
    else {
        m_tuple_track = ntupleSvc()->book("LineFit/Track", CLID_ColumnWiseTuple, "N-Tuple example");
        if (m_tuple_track) {
            sc = m_tuple_track->addItem("run", tr_run);
            sc = m_tuple_track->addItem("event", tr_event);
            sc = m_tuple_track->addItem("drho", tr_drho);
            sc = m_tuple_track->addItem("phi0", tr_phi0);
            sc = m_tuple_track->addItem("dz0", tr_dz0);
            sc = m_tuple_track->addItem("chisq", tr_chisq);
            sc = m_tuple_track->addItem("Is_fitted", tr_Is_fitted);
            sc = m_tuple_track->addItem("tanLam", tr_tanLam);
            sc = m_tuple_track->addItem("ncluster", tr_ncluster, 0, 5000);
            sc = m_tuple_track->addIndexedItem("id_cluster", tr_ncluster, tr_id_cluster);
            sc = m_tuple_track->addIndexedItem("x_cluster", tr_ncluster, tr_x_cluster);
            sc = m_tuple_track->addIndexedItem("y_cluster", tr_ncluster, tr_y_cluster);
            sc = m_tuple_track->addIndexedItem("z_cluster", tr_ncluster, tr_z_cluster);
            sc = m_tuple_track->addIndexedItem("Q_cluster", tr_ncluster, tr_Q_cluster);
            sc = m_tuple_track->addIndexedItem("phi_cluster", tr_ncluster, tr_phi_cluster);
            sc = m_tuple_track->addIndexedItem("layer_cluster", tr_ncluster, tr_layer_cluster);
            sc = m_tuple_track->addIndexedItem("sheet_cluster", tr_ncluster, tr_sheet_cluster);
            sc = m_tuple_track->addIndexedItem("Is_selected_cluster", tr_ncluster, tr_Is_selected_cluster);
 
        } else {
            log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple_track) << endmsg;
            return StatusCode::FAILURE;
        }
    }
 
    IRawDataProviderSvc *irawDataProviderSvc;
    sc = service("RawDataProviderSvc", irawDataProviderSvc);
    m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc *>(irawDataProviderSvc);
    if (sc.isFailure()) {
        log << MSG::FATAL << name() << " Could not load RawDataProviderSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
    ICgemGeomSvc *icgemGeomSvc;
    sc = service("CgemGeomSvc", icgemGeomSvc);
    m_cgemGeomSvc = dynamic_cast<CgemGeomSvc *>(icgemGeomSvc);
 
    if (sc.isFailure()) {
        log << MSG::FATAL << "Could not load CgemGeomSvc!" << endreq;
        return StatusCode::FAILURE;
    }
    myNCgemLayers = m_cgemGeomSvc->getNumberOfCgemLayer();
    for (int i = 0; i < myNCgemLayers; i++) {
        myCgemLayer[i] = m_cgemGeomSvc->getCgemLayer(i);
        myNSheets[i] = myCgemLayer[i]->getNumberOfSheet();
        myRXstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu2();
        myRVstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu1();
    }
    R_layer[0] = (m_cgemGeomSvc->getCgemLayer(0)->getMiddleROfGapD());
    R_layer[1] = (m_cgemGeomSvc->getCgemLayer(1)->getMiddleROfGapD());
    R_layer[2] = (m_cgemGeomSvc->getCgemLayer(2)->getMiddleROfGapD());
    AngleOfStereo[0] = (m_cgemGeomSvc->getCgemLayer(0)->getAngleOfStereo());
    AngleOfStereo[1] = (m_cgemGeomSvc->getCgemLayer(1)->getAngleOfStereo());
    AngleOfStereo[2] = (m_cgemGeomSvc->getCgemLayer(2)->getAngleOfStereo());
    length = m_cgemGeomSvc->getLengthOfCgem();
    Mp = m_cgemGeomSvc->getMidDriftPtr();
    Al = m_cgemGeomSvc->getAlignPtr();
 
    GeoLayer0 = m_cgemGeomSvc->getCgemLayer(0);
    GeoLayer1 = m_cgemGeomSvc->getCgemLayer(1);
    GeoLayer2 = m_cgemGeomSvc->getCgemLayer(2);
 
    pl_00 = m_cgemGeomSvc->getReadoutPlane(0, 0);
    pl_10 = m_cgemGeomSvc->getReadoutPlane(1, 0);
    pl_11 = m_cgemGeomSvc->getReadoutPlane(1, 1);
    pl_20 = m_cgemGeomSvc->getReadoutPlane(2, 0);
    pl_21 = m_cgemGeomSvc->getReadoutPlane(2, 1);
 
    cout << "CgemLineFit alignment: is it on? " << Align_Flag << endl;
    if (Align_Flag) {
        for (int ilay = 0; ilay < 6; ilay++) {
            cout << "LAYER " << ilay + 1 << endl;
            cout << "shift dx " << Al->getDx(ilay) << " dy " << Al->getDy(ilay) << " dz " << Al->getDz(ilay) << endl;
            cout << "rotation Rx " << Al->getRx(ilay) << " Ry " << Al->getRy(ilay) << " Rz " << Al->getRz(ilay) << endl;
            cout << "midplane radius " << Mp->getR(int(ilay / 2)) << endl;
        }
    }
 
    ICgemCalibFunSvc *icgemCalibSvc;
    sc = service("CgemCalibFunSvc", icgemCalibSvc);
    m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc *>(icgemCalibSvc);
    if (sc.isFailure()) {
        log << MSG::FATAL << "Could not load CgemCalibFunSvc!" << endreq;
        return StatusCode::FAILURE;
    }
    sc = service("CgemCalibFunSvc", myCgemCalibSvc);
    if (sc.isFailure()) {
        cout << name() << "Could not load MdcCalibFunSvc!" << endl;
        return sc;
    }
    // cout<<"sigma from CgemCalibFunSvc: "<<myCgemCalibSvc->getSigma(0,0,0,0,0,0)<<endl;
    lutreader = new CgemLUTReader(m_LUTfile);
    lutreader->ReadLUT();
 
    return StatusCode::SUCCESS;
}

Layer_cross()

bool CgemLineFit::Layer_cross	(	int	R_id,
		StraightLine *	l1 )

Definition at line 4007 of file CgemLineFit.cxx.

                                                        {
 
    // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
    double phiV_up[2], phiV_down[2];
 
    HepPoint3D Up, Down;
    if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
    else Mp->getPointAligned(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
    if (Up[0] > 9999 || Up[1] > 9999 || Down[0] > 9999 || Down[1] > 9999) return false;
    else return true;
}

Referenced by Data_Closest(), and Filter().

Loop_All()

bool CgemLineFit::Loop_All

(

)

Definition at line 1295 of file CgemLineFit.cxx.

                           {
    double cl_00(0), cl_01(0), cl_11(0), cl_10(0);
    int _loop(0);
    double maxQ_11(0), maxQ_10(0), maxQ_00(0);
    int max_11(-1), max_10(-1), max_00(-1);
    double C_00(0), C_10(0), C_11(0), C_01(0);
    vector<double> Vec_00_layer, Vec_00_phi, Vec_00_Z, Vec_00_layerid, Vec_00_v, Vec_00_flag, Vec_00_Q, Vec_00_sheetid; //
    vector<double> Vec_01_layer, Vec_01_phi, Vec_01_Z, Vec_01_layerid, Vec_01_v, Vec_01_flag, Vec_01_Q, Vec_01_sheetid; //
    vector<double> Vec_10_layer, Vec_10_phi, Vec_10_Z, Vec_10_layerid, Vec_10_v, Vec_10_flag, Vec_10_Q, Vec_10_sheetid; //
    vector<double> Vec_11_layer, Vec_11_phi, Vec_11_Z, Vec_11_layerid, Vec_11_v, Vec_11_flag, Vec_11_Q, Vec_11_sheetid; //
 
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
    RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
    vector<int> L_00, L_01, L_10, L_11;
    int ic(0);
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
 
        int flag = cluster->getflag();
        int layerid = cluster->getlayerid();
        double _Q = cluster->getenergydeposit();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        int sheetid = cluster->getsheetid();
        int clusterid = cluster->getclusterid();
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        double x = R_layer[layerid] * cos(_phi);
        double y = R_layer[layerid] * sin(_phi);
        double z = cluster->getRecZ();
        if (_loop > MAX_Clusters) break;
        ic++;
        if (flag != 0) continue;
 
        if (layerid == 1 && sheetid == 1 && TEST_N != 1) { L_11.push_back(_loop); }
        if (layerid == 1 && sheetid == 0 && TEST_N != 4) { L_10.push_back(_loop); }
        if (layerid == 0 && sheetid == 0 && _phi > 0 && TEST_N != 2) { L_01.push_back(_loop); }
        if (layerid == 0 && sheetid == 0 && _phi < 0 && TEST_N != 3) { L_00.push_back(_loop); }
    }
 
    if (TEST_N == 1) L_11.push_back(-1);
    if (TEST_N == 2) L_01.push_back(-1);
    if (TEST_N == 3) L_00.push_back(-1);
    if (TEST_N == 4) L_10.push_back(-1);
 
    double Min_chi(1E10);
 
    if ((L_00.size() * L_01.size() * L_10.size() * L_11.size()) > MAX_COMB) { return false; }
 
    NXComb = L_00.size() * L_01.size() * L_10.size() * L_11.size();
 
    for (int i_11 = 0; i_11 < L_11.size(); i_11++) {
        for (int i_10 = 0; i_10 < L_10.size(); i_10++) {
            for (int i_00 = 0; i_00 < L_00.size(); i_00++) {
                for (int i_01 = 0; i_01 < L_01.size(); i_01++) {
 
                    Vec_layer.clear();
                    Vec_phi.clear();
                    Vec_Z.clear();
                    Vec_v.clear();
                    Vec_layerid.clear();
                    Vec_Virlayerid.clear();
                    Vec_flag.clear();
                    Vec_Q.clear();
                    Vec_sheetid.clear();
                    Vec_m_layer.clear();
                    Vec_m_phi.clear();
                    Vec_m_Z.clear();
                    Vec_m_v.clear();
                    Vec_m_layerid.clear();
                    Vec_m_flag.clear();
                    Vec_m_Q.clear();
                    Vec_m_sheetid.clear();
                    _loop = 0;
                    iter = recCgemClusterCol->begin();
                    for (; iter != recCgemClusterCol->end(); iter++) {
                        _loop++;
                        RecCgemCluster *cluster = (*iter);
                        int flag = cluster->getflag();
                        if (!(_loop == L_11[i_11] || _loop == L_10[i_10] || _loop == L_00[i_00] || _loop == L_01[i_01])) continue;
                        int layerid = cluster->getlayerid();
                        int sheetid = cluster->getsheetid();
                        double _phi = cluster->getrecphi();
                        double _v = cluster->getrecv();
                        double _Z = cluster->getRecZ();
                        double _Q = cluster->getenergydeposit();
                        double t_phi = cluster->getrecphi_mTPC();
                        double t_v = cluster->getrecv_mTPC();
                        double t_Z = cluster->getRecZ_mTPC();
 
                        double merged_phi = cluster->get_merge_phi();
                        double merged_v = cluster->get_merge_v();
 
                        int vir_layerid = GetVirLay(layerid, _phi);
 
                        Vec_flag.push_back(flag);
                        Vec_layerid.push_back(layerid);
                        Vec_Virlayerid.push_back(vir_layerid);
                        Vec_sheetid.push_back(sheetid);
                        Vec_layer.push_back(R_layer[layerid]);
 
                        Vec_m_phi.push_back(t_phi);
                        Vec_m_v.push_back(t_v);
                        Vec_m_Z.push_back(t_Z);
 
                        if (Switch_CCmTPC == 0) {
                            Vec_phi.push_back(_phi);
                            Vec_v.push_back(_v);
                            Vec_Z.push_back(_Z);
 
                        } else if (Switch_CCmTPC == 1) {
                            Vec_phi.push_back(t_phi);
                            Vec_v.push_back(t_v);
                            Vec_Z.push_back(t_Z);
 
                        } else if (Switch_CCmTPC == 2) {
 
                            CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
                            double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
                            Vec_phi.push_back(merged_phi);
                            Vec_v.push_back(merged_v);
                            Vec_Z.push_back(merged_Z);
                        }
 
                    } // one combo
                    if (Vec_layerid.size() != 4 && TEST_N == 0) continue;
                    else if (Vec_layerid.size() != 3 && TEST_N > 0) continue;
 
                    OrderClusters();
 
                    StraightLine *l_outer_tmp;
 
                    if (TEST_N == 1 || TEST_N == 4) l_outer_tmp = IniPar(Vec_phi[1], Vec_Z[1], 1, Vec_phi[2], Vec_Z[2], 0);
                    else if (TEST_N == 2 || TEST_N == 3) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
                    else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
 
                    TMinuit *_fit = Fit(l_outer_tmp, fa, 1);
                    double _chi = _fit->fAmin;
                    if (_chi < Min_chi) {
                        Min_chi = _chi;
                        if (TEST_N == 0) C_00 = Vec_phi[3];
                        else C_00 = -99;
                        C_01 = Vec_phi[2];
                        C_10 = Vec_phi[1];
                        C_11 = Vec_phi[0];
                    }
                    delete _fit;
                    delete l_outer_tmp;
                }
            }
        }
    }
 
    iter = recCgemClusterCol->begin();
    L_00.clear();
    L_01.clear();
    L_10.clear();
    L_11.clear();
    _loop = 0;
    int ic_tot(0);
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        if (flag != 2) continue;
        ic_tot++;
        int layerid = cluster->getlayerid();
        double _Q = cluster->getenergydeposit();
 
        double _phi = 99;
        if (Switch_CCmTPC == 0) _phi = cluster->getrecphi();
        else if (Switch_CCmTPC == 1) _phi = cluster->getrecphi_mTPC();
        else if (Switch_CCmTPC == 2) {
            RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
            RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
            _phi = cluster_x->get_merge_phi();
        }
 
        int sheetid = cluster->getsheetid();
        if (_Q < MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
        if (C_00 == _phi || C_01 == _phi || C_11 == _phi || C_10 == _phi) {
            if (layerid == 1 && sheetid == 1) { L_11.push_back(_loop); }
            if (layerid == 1 && sheetid == 0) { L_10.push_back(_loop); }
            if (layerid == 0 && sheetid == 0 && _phi > 0) { L_01.push_back(_loop); }
            if (layerid == 0 && sheetid == 0 && _phi < 0) { L_00.push_back(_loop); }
        }
    }
    Min_chi = 1E10;
 
    if (TEST_N == 1) L_11.push_back(-1);
    if (TEST_N == 2) L_01.push_back(-1);
    if (TEST_N == 3) L_00.push_back(-1);
    if (TEST_N == 4) L_10.push_back(-1);
 
    clst_11 = L_11.size();
    clst_01 = L_01.size();
    clst_10 = L_10.size();
    clst_00 = L_00.size();
    if ((L_00.size() * L_01.size() * L_10.size() * L_11.size()) > MAX_COMB) {
        cout << "xv size much " << endl;
        return false;
    }
 
    NVComb = L_00.size() * L_01.size() * L_10.size() * L_11.size();
 
    for (int i_11 = 0; i_11 < L_11.size(); i_11++) {
        for (int i_10 = 0; i_10 < L_10.size(); i_10++) {
            for (int i_00 = 0; i_00 < L_00.size(); i_00++) {
                for (int i_01 = 0; i_01 < L_01.size(); i_01++) {
 
                    Vec_layer.clear();
                    Vec_phi.clear();
                    Vec_Z.clear();
                    Vec_v.clear();
                    Vec_layerid.clear();
                    Vec_Virlayerid.clear();
                    Vec_flag.clear();
                    Vec_Q.clear();
                    Vec_sheetid.clear();
                    Vec_m_layer.clear();
                    Vec_m_phi.clear();
                    Vec_m_Z.clear();
                    Vec_m_v.clear();
                    Vec_m_layerid.clear();
                    Vec_m_flag.clear();
                    Vec_m_Q.clear();
                    Vec_m_sheetid.clear();
                    _loop = 0;
                    iter = recCgemClusterCol->begin();
                    for (; iter != recCgemClusterCol->end(); iter++) {
                        _loop++;
                        RecCgemCluster *cluster = (*iter);
                        int flag = cluster->getflag();
                        if (flag != 2) continue;
 
                        int layerid = cluster->getlayerid();
                        int sheetid = cluster->getsheetid();
                        double _phi = cluster->getrecphi();
                        double _v = cluster->getrecv();
                        double _Z = cluster->getRecZ();
                        double _Q = cluster->getenergydeposit();
                        double t_phi = cluster->getrecphi_mTPC();
                        double t_v = cluster->getrecv_mTPC();
                        double t_Z = cluster->getRecZ_mTPC();
 
                        RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
                        RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
                        double merged_phi = cluster_x->get_merge_phi();
                        RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
                        double merged_v = cluster_v->get_merge_v();
 
                        int vir_layerid = GetVirLay(layerid, _phi);
 
                        if (!(_loop == L_11[i_11] || _loop == L_10[i_10] || _loop == L_00[i_00] || _loop == L_01[i_01])) continue;
 
                        Vec_flag.push_back(flag);
                        Vec_layerid.push_back(layerid);
                        Vec_Virlayerid.push_back(vir_layerid);
                        Vec_sheetid.push_back(sheetid);
                        Vec_layer.push_back(R_layer[layerid]);
 
                        Vec_m_phi.push_back(t_phi);
                        Vec_m_v.push_back(t_v);
                        Vec_m_Z.push_back(t_Z);
 
                        if (Switch_CCmTPC == 0) {
                            Vec_phi.push_back(_phi);
                            Vec_v.push_back(_v);
                            Vec_Z.push_back(_Z);
 
                        } else if (Switch_CCmTPC == 1) {
                            Vec_phi.push_back(t_phi);
                            Vec_v.push_back(t_v);
                            Vec_Z.push_back(t_Z);
 
                        } else if (Switch_CCmTPC == 2) {
 
                            CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
                            double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
                            Vec_phi.push_back(merged_phi);
                            Vec_v.push_back(merged_v);
                            Vec_Z.push_back(merged_Z);
                        }
 
                    } // one combo
 
                    if (Vec_layerid.size() != 4 && TEST_N == 0) continue;
                    else if (Vec_layerid.size() != 3 && TEST_N > 0) continue;
 
                    OrderClusters();
 
                    if (Vec_v[0] == Vec_v[1] || Vec_v[0] == Vec_v[2] || Vec_v[0] == Vec_v[3] || Vec_v[1] == Vec_v[2] ||
                            Vec_v[1] == Vec_v[3] || Vec_v[2] == Vec_v[3])
                        continue;
                    StraightLine *l_outer_tmp;
 
                    if (TEST_N == 1 || TEST_N == 4) l_outer_tmp = IniPar(Vec_phi[1], Vec_Z[1], 1, Vec_phi[2], Vec_Z[2], 0);
                    else if (TEST_N == 2 || TEST_N == 3) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
                    else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
 
                    TMinuit *_fit = Fit(l_outer_tmp, fa, 2);
                    double _chi = _fit->fAmin;
                    if (_chi < Min_chi) {
                        Min_chi = _chi;
                        C_00 = L_00[i_00];
                        C_01 = L_01[i_01];
                        C_10 = L_10[i_10];
                        C_11 = L_11[i_11];
                    }
                    delete _fit;
                    delete l_outer_tmp;
                }
            }
        }
    }
 
    Vec_layer.clear();
    Vec_phi.clear();
    Vec_Z.clear();
    Vec_v.clear();
    Vec_layerid.clear();
    Vec_Virlayerid.clear();
    Vec_flag.clear();
    Vec_Q.clear();
    Vec_sheetid.clear();
    Vec_m_layer.clear();
    Vec_m_phi.clear();
    Vec_m_Z.clear();
    Vec_m_v.clear();
    Vec_m_layerid.clear();
    Vec_m_flag.clear();
    Vec_m_Q.clear();
    Vec_m_sheetid.clear();
    _loop = 0;
    ic = 0;
    iter = recCgemClusterCol->begin();
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        int layerid = cluster->getlayerid();
        int sheetid = cluster->getsheetid();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        double _Z = cluster->getRecZ();
        double _Q = cluster->getenergydeposit();
        double t_phi = cluster->getrecphi_mTPC();
        double t_v = cluster->getrecv_mTPC();
        double t_Z = cluster->getRecZ_mTPC();
        int clusterid = cluster->getclusterid();
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        double x = R_layer[layerid] * cos(_phi);
        double y = R_layer[layerid] * sin(_phi);
        double z = cluster->getRecZ();
        int vir_layerid = GetVirLay(layerid, _phi);
 
        if (_loop > MAX_Clusters) break;
        if (flag != 2) continue;
 
        RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
        RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
        double merged_phi = cluster_x->get_merge_phi();
        RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
        double merged_v = cluster_v->get_merge_v();
 
        tr_id_cluster[ic] = clusterid;
        tr_x_cluster[ic] = x;
        tr_y_cluster[ic] = y;
        tr_z_cluster[ic] = z;
        tr_Q_cluster[ic] = _Q;
        // tr_phi_cluster[ic] = _phi;
        // tr_layer_cluster[ic] = layerid;
        // tr_sheet_cluster[ic] = sheetid;
        tr_Is_selected_cluster[ic] = 0;
        ic++;
        if (!(_loop == C_00 || _loop == C_10 || _loop == C_01 || _loop == C_11)) continue;
 
        tr_Is_selected_cluster[ic - 1] = 1;
 
        Vec_flag.push_back(flag);
        Vec_layerid.push_back(layerid);
        Vec_Virlayerid.push_back(vir_layerid);
        Vec_sheetid.push_back(sheetid);
        Vec_layer.push_back(R_layer[layerid]);
 
        Vec_m_phi.push_back(t_phi);
        Vec_m_v.push_back(t_v);
        Vec_m_Z.push_back(t_Z);
 
        if (Switch_CCmTPC == 0) {
            Vec_phi.push_back(_phi);
            Vec_v.push_back(_v);
            Vec_Z.push_back(_Z);
 
        } else if (Switch_CCmTPC == 1) {
            Vec_phi.push_back(t_phi);
            Vec_v.push_back(t_v);
            Vec_Z.push_back(t_Z);
 
        } else if (Switch_CCmTPC == 2) {
 
            CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
            double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
            Vec_phi.push_back(merged_phi);
            Vec_v.push_back(merged_v);
            Vec_Z.push_back(merged_Z);
        }
 
        cluster->setTrkId(0);
        vecclusters.push_back(cluster);
 
    } // one combo
    NC = vecclusters.size();
    tr_ncluster = ic;
 
    for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
        if (Vec_layerid[i_cl] == 2) _clst_2++;
        if (Vec_layerid[i_cl] == 1) _clst_1++;
        if (Vec_layerid[i_cl] == 0) _clst_0++;
    }
 
    if (Vec_layerid.size() != 4 && TEST_N == 0) return false;
    else if (Vec_layerid.size() != 3 && TEST_N > 0) return false;
 
    OrderClusters();
 
    if (TEST_N == 1 || TEST_N == 4) l_outer = IniPar(Vec_phi[1], Vec_Z[1], 1, Vec_phi[2], Vec_Z[2], 0);
    else if (TEST_N == 2 || TEST_N == 3) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
    else if (TEST_N == 0) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
 
    if (Min_chi < Chisq_cut) return true;
    else return false;
}

Referenced by execute().

Loop_MaxQ()

bool CgemLineFit::Loop_MaxQ

(

)

Definition at line 1735 of file CgemLineFit.cxx.

                            {
    // modified by wxn
    double cl_00(0), cl_01(0), cl_11(0), cl_10(0), cl_20(0), cl_21(0);
    int _loop(0);
    double maxQ_21(0), maxQ_20(0), maxQ_11(0), maxQ_10(0), maxQ_00(0);
    int max_21(-1), max_20(-1), max_11(-1), max_10(-1), max_00(-1);
    double C_00(0), C_10(0), C_11(0), C_01(0), C_20(0), C_21(0);
    int iC_00(0), iC_10(0), iC_11(0), iC_01(0), iC_20(0), iC_21(0);
    if (run < 0) Get_MCInf();
    vector<double> Vec_00_layer, Vec_00_phi, Vec_00_Z, Vec_00_layerid, Vec_00_v, Vec_00_flag, Vec_00_Q, Vec_00_sheetid; //
    vector<double> Vec_01_layer, Vec_01_phi, Vec_01_Z, Vec_01_layerid, Vec_01_v, Vec_01_flag, Vec_01_Q, Vec_01_sheetid; //
    vector<double> Vec_10_layer, Vec_10_phi, Vec_10_Z, Vec_10_layerid, Vec_10_v, Vec_10_flag, Vec_10_Q, Vec_10_sheetid; //
    vector<double> Vec_11_layer, Vec_11_phi, Vec_11_Z, Vec_11_layerid, Vec_11_v, Vec_11_flag, Vec_11_Q, Vec_11_sheetid; //
    vector<double> Vec_20_layer, Vec_20_phi, Vec_20_Z, Vec_20_layerid, Vec_20_v, Vec_20_flag, Vec_20_Q, Vec_20_sheetid; //
    vector<double> Vec_21_layer, Vec_21_phi, Vec_21_Z, Vec_21_layerid, Vec_21_v, Vec_21_flag, Vec_21_Q, Vec_21_sheetid; //
 
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
    RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
    vector<int> L_00, L_01, L_10, L_11, L_20, L_21;
    vector<double> Q_00, Q_01, Q_10, Q_11, Q_20, Q_21;
    vector<int> ZL_00, ZL_01, ZL_10, ZL_11, ZL_20, ZL_21;
    set<int> ZL_00_s, ZL_01_s, ZL_10_s, ZL_11_s, ZL_20_s, ZL_21_s;
 
    int ic(0);
    if (debug) cout << "start Loop_MaxQ method" << endl;
 
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
 
        int flag = cluster->getflag();
        if (flag != 2) continue;
        int layerid = cluster->getlayerid();
        double _Q = cluster->getenergydeposit();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        int sheetid = cluster->getsheetid();
        int clusterid = cluster->getclusterid();
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        double x = R_layer[layerid] * cos(_phi);
        double y = R_layer[layerid] * sin(_phi);
        double z = cluster->getRecZ();
        
        // --- cuts Q>10fC, size>1
        //RecCgemCluster* clusterX=*(recCgemClusterCol->begin()+clusterflagb);
        //RecCgemCluster* clusterV=*(recCgemClusterCol->begin()+clusterflage);
        //double QX=clusterX->getenergydeposit();
        //if(QX<10.) continue;
        //double QV=clusterV->getenergydeposit();
        //if(QV<10.) continue;
        //int size_X = 1+abs(clusterX->getclusterflage()-clusterX->getclusterflagb());
        //if(size_X<2) continue;
        //int size_V = 1+abs(clusterV->getclusterflage()-clusterV->getclusterflagb());
        //if(size_V<2) continue;
        
        if (debug) cout << "layerid is " << layerid << endl;
        if (_loop > MAX_Clusters) break;
        ic++;
        if (layerid == 2 && sheetid == 1) {
            ZL_21_s.insert(clusterflagb);
            // is_cluster_exsited[0]++;
        }
        if (layerid == 2 && sheetid == 0) {
            ZL_20_s.insert(clusterflagb);
            // is_cluster_exsited[1]++;
        }
        if (layerid == 1 && sheetid == 1) {
            ZL_11_s.insert(clusterflagb);
            // is_cluster_exsited[2]++;
        }
        if (layerid == 1 && sheetid == 0) {
            ZL_10_s.insert(clusterflagb);
            // is_cluster_exsited[3]++;
        }
        if (layerid == 0 && sheetid == 0 && _phi > 0) {
            ZL_01_s.insert(clusterflagb);
            // is_cluster_exsited[4]++;
        }
        if (layerid == 0 && sheetid == 0 && _phi < 0) {
            ZL_00_s.insert(clusterflagb);
            // is_cluster_exsited[5]++;
        }
    }
    ZL_21.assign(ZL_21_s.begin(), ZL_21_s.end());
    ZL_20.assign(ZL_20_s.begin(), ZL_20_s.end());
    ZL_11.assign(ZL_11_s.begin(), ZL_11_s.end());
    ZL_10.assign(ZL_10_s.begin(), ZL_10_s.end());
    ZL_01.assign(ZL_01_s.begin(), ZL_01_s.end());
    ZL_00.assign(ZL_00_s.begin(), ZL_00_s.end());
    ////
    // if (debug) cout << "ZL_21.size() : " << ZL_21.size() << endl;
    // if (is_cluster_exsited[0]>0&& is_cluster_exsited[1]>0&& is_cluster_exsited[2]>0&& is_cluster_exsited[3]>0&&
    // is_cluster_exsited[4]>0&& is_cluster_exsited[5] > 0)
    //{
    //   nEvent_with_6cluster++;
    // }
 
    // if (is_cluster_exsited[0]>0&& is_cluster_exsited[1]>0&& is_cluster_exsited[2]>0&& is_cluster_exsited[3]> 0)
    //{
    //   nEvent_with_4cluster++;
    // }
 
    ////
    iter = recCgemClusterCol->begin();
    _loop = 0;
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        int layerid = cluster->getlayerid();
        int sheetid = cluster->getsheetid();
        double _Q = cluster->getenergydeposit();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        double _Z = cluster->getRecZ();
        int clusterid = cluster->getclusterid();
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        if (flag != 0) continue;
        if (layerid == 2 && sheetid == 1) {
            for (int i_21 = 0; i_21 < ZL_21.size(); i_21++) {
                if (clusterid == ZL_21[i_21]) {
                    L_21.push_back(_loop);
                    Q_21.push_back(_Q);
                }
            }
        }
        if (layerid == 2 && sheetid == 0) {
            for (int i_20 = 0; i_20 < ZL_20.size(); i_20++) {
                if (clusterid == ZL_20[i_20]) {
                    L_20.push_back(_loop);
                    Q_20.push_back(_Q);
                }
            }
        }
        if (layerid == 1 && sheetid == 1) {
            for (int i_11 = 0; i_11 < ZL_11.size(); i_11++) {
                if (clusterid == ZL_11[i_11]) {
                    L_11.push_back(_loop);
                    Q_11.push_back(_Q);
                }
            }
        }
        if (layerid == 1 && sheetid == 0) {
            for (int i_10 = 0; i_10 < ZL_10.size(); i_10++) {
                if (clusterid == ZL_10[i_10]) {
                    L_10.push_back(_loop);
                    Q_10.push_back(_Q);
                }
            }
        }
        if (layerid == 0 && sheetid == 0 && _phi > 0) {
            for (int i_01 = 0; i_01 < ZL_01.size(); i_01++) {
                if (clusterid == ZL_01[i_01]) {
                    L_01.push_back(_loop);
                    Q_01.push_back(_Q);
                }
            }
        }
        if (layerid == 0 && sheetid == 0 && _phi < 0) {
            for (int i_00 = 0; i_00 < ZL_00.size(); i_00++) {
                if (clusterid == ZL_00[i_00]) {
                    L_00.push_back(_loop);
                    Q_00.push_back(_Q);
                }
            }
        }
    }
    // chose the N X clusters with the largest charge  and save corresponding
    // index
    vector<int> L_00_s, L_01_s, L_10_s, L_11_s, L_20_s, L_21_s;
 
    L_00_s = GetNMaxQ(Q_00, L_00, Nmax);
    L_01_s = GetNMaxQ(Q_01, L_01, Nmax);
    L_10_s = GetNMaxQ(Q_10, L_10, Nmax);
    L_11_s = GetNMaxQ(Q_11, L_11, Nmax);
    L_20_s = GetNMaxQ(Q_20, L_20, Nmax);
    L_21_s = GetNMaxQ(Q_21, L_21, Nmax);
    ////
    if (debug) cout << "L_21_s.size() : " << L_21_s.size() << endl;
    ////
 
    if (TEST_N == 1) {
        L_21_s.clear();
        L_21_s.push_back(-1);
    } else if (TEST_N == 2) {
        L_11_s.clear();
        L_11_s.push_back(-1);
    } else if (TEST_N == 3) {
        L_01_s.clear();
        L_01_s.push_back(-1);
    } else if (TEST_N == 4) {
        L_00_s.clear();
        L_00_s.push_back(-1);
    } else if (TEST_N == 5) {
        L_10_s.clear();
        L_10_s.push_back(-1);
    } else if (TEST_N == 6) {
        L_20_s.clear();
        L_20_s.push_back(-1);
    }
 
    if (debug) cout << "finish the first loop" << endl;
 
    double Min_chi(1E10);
    NXComb = L_00_s.size() * L_01_s.size() * L_10_s.size() * L_11_s.size() * L_20_s.size() * L_21_s.size();
 
    int com0, com1, com2, com3, com4, com5;
    for (int i_21 = 0; i_21 < L_21_s.size(); i_21++) {
        for (int i_20 = 0; i_20 < L_20_s.size(); i_20++) {
            for (int i_11 = 0; i_11 < L_11_s.size(); i_11++) {
                for (int i_10 = 0; i_10 < L_10_s.size(); i_10++) {
                    for (int i_00 = 0; i_00 < L_00_s.size(); i_00++) {
                        for (int i_01 = 0; i_01 < L_01_s.size(); i_01++) {
 
                            Vec_layer.clear();
                            Vec_phi.clear();
                            Vec_Z.clear();
                            Vec_v.clear();
                            Vec_layerid.clear();
                            Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
                            Vec_flag.clear();
                            Vec_Q.clear();
                            Vec_sheetid.clear();
                            Vec_m_layer.clear();
                            Vec_m_phi.clear();
                            Vec_m_Z.clear();
                            Vec_m_v.clear();
                            Vec_m_layerid.clear();
                            Vec_m_flag.clear();
                            Vec_m_Q.clear();
                            Vec_m_sheetid.clear();
                            _loop = 0;
                            iter = recCgemClusterCol->begin();
 
                            vector<int> Com;
                            Com.clear();
                            for (; iter != recCgemClusterCol->end(); iter++) {
                                _loop++;
                                RecCgemCluster *cluster = (*iter);
                                int flag = cluster->getflag();
                                /////
                                // if (debug)  cout << "(wxn)cluster flag is " << flag << endl;
                                /////
                                if (!(_loop == L_21_s[i_21] || _loop == L_20_s[i_20] || _loop == L_11_s[i_11] || _loop == L_10_s[i_10] ||
                                            _loop == L_00_s[i_00] || _loop == L_01_s[i_01]))
                                    continue;
                                /////
                                // if (debug) cout << "(wxn)_loop is " << _loop << endl;
                                /////
 
                                int layerid = cluster->getlayerid();
                                int sheetid = cluster->getsheetid();
                                double _phi = cluster->getrecphi();
                                double _v = cluster->getrecv();
                                double _Z = cluster->getRecZ();
                                double _Q = cluster->getenergydeposit();
                                double t_phi = cluster->getrecphi_mTPC();
                                double t_v = cluster->getrecv_mTPC();
                                double t_Z = cluster->getRecZ_mTPC();
                                double merged_phi = cluster->get_merge_phi();
                                double merged_v = cluster->get_merge_v();
                                /////
                                // cout << "(wxn)mergedv is " << merged_v << endl;
                                /////
 
                                int vir_layerid = GetVirLay(layerid, _phi);
 
                                Vec_flag.push_back(flag);
                                Vec_layerid.push_back(layerid);
                                Vec_Virlayerid.push_back(vir_layerid);
                                Vec_sheetid.push_back(sheetid);
                                Vec_layer.push_back(R_layer[layerid]);
                                Com.push_back(_loop);
 
                                Vec_m_phi.push_back(t_phi);
                                Vec_m_v.push_back(t_v);
                                Vec_m_Z.push_back(t_Z);
 
                                if (Switch_CCmTPC == 0) {
                                    Vec_phi.push_back(_phi);
                                    Vec_v.push_back(_v);
                                    Vec_Z.push_back(_Z);
                                    if (debug) cout << "1st phi, v, by merge: " << merged_phi << ", " << merged_v << endl;
                                    if (debug) cout << "1st phi, v, by CC: " << _phi << ", " << _v << endl;
                                } else if (Switch_CCmTPC == 1) {
                                    Vec_phi.push_back(t_phi);
                                    Vec_v.push_back(t_v);
                                    Vec_Z.push_back(t_Z);
                                } else if (Switch_CCmTPC == 2) {
 
                                    CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
                                    double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
                                    Vec_phi.push_back(merged_phi);
                                    Vec_v.push_back(merged_v);
                                    Vec_Z.push_back(merged_Z);
                                    if (debug) cout << "1st phi, v, by merge: " << merged_phi << ", " << merged_v << endl;
                                    if (debug) cout << "1st phi, v, by CC: " << _phi << ", " << _v << endl;
                                }
 
                            } // one combo
                            if (Vec_layerid.size() != 6 && TEST_N == 0) return false;
                            else if (Vec_layerid.size() != 5 && TEST_N > 0) return false;
                            if (debug) cout << "before OrderCluster" << endl;
                            OrderClusters();
                            //                     0     1             2   3
                            // order is layer 1 , top, down, layer 0, top down
                            // ????? How to order layer 2 ?????
                            // OrderClusters() need modified !!!!!
                            //                     0     1             2   3              4   5
                            // order is layer 2 , top, down, layer 1, top down, layer 0, top down
 
                            StraightLine *l_outer_tmp;
 
                            if (debug) cout << "before IniPar" << endl;
                            //                                                 lay0top   virlay
                            //                                                 lay0down  virlay
                            if (TEST_N == 1 || TEST_N == 6) l_outer_tmp = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
                            else if (TEST_N == 2 || TEST_N == 3 || TEST_N == 4 || TEST_N == 5)
                                l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
                            else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
 
                            if (debug) cout << "before fit" << endl;
                            TMinuit *_fit = Fit(l_outer_tmp, fa, 1);
 
                            if (debug) cout << "chisq is : " << _fit->fAmin << endl;
                            if (debug) cout << "===========================================" << endl;
 
                            double _chi = _fit->fAmin;
                            if (debug) cout << "finish the fit and the chisq is " << _chi << endl;
                            if (debug) cout << "----------------------------------------" << endl;
                            if (_chi < Min_chi) {
                                Min_chi = _chi;
                                if (TEST_N == 0) C_00 = Vec_phi[5];
                                else C_00 = -99;   // Why set C_00 to -99 Guoaq--2020/11/01
                                C_01 = Vec_phi[4]; // It is because there is only 3 clusters
                                // selected if TEST_N!=0
                                C_10 = Vec_phi[3];
                                C_11 = Vec_phi[2];
                                C_20 = Vec_phi[1];
                                C_21 = Vec_phi[0];
                                com0 = Com[0];
                                com1 = Com[1];
                                com2 = Com[2];
                                com3 = Com[3];
                                com4 = Com[4];
                                com5 = Com[5];
                            }
                            delete _fit;
                            delete l_outer_tmp;
                        }
                    }
                }
            }
        }
    }
 
    if (debug) cout << "finish the second loop" << endl;
 
    iter = recCgemClusterCol->begin();
    ZL_00.clear();
    ZL_01.clear();
    ZL_10.clear();
    ZL_11.clear();
    ZL_20.clear();
    ZL_21.clear();
    L_00.clear();
    L_01.clear();
    L_10.clear();
    L_11.clear();
    L_20.clear();
    L_21.clear();
    L_00_s.clear();
    L_01_s.clear();
    L_10_s.clear();
    L_11_s.clear();
    L_20_s.clear();
    L_21_s.clear();
    Q_00.clear();
    Q_01.clear();
    Q_10.clear();
    Q_11.clear();
    Q_20.clear();
    Q_21.clear();
    _loop = 0;
    int ic_tot(0);
 
    if (debug) cout << "c00, c01, c11, c10 : " << C_00 << ", " << C_01 << ", " << C_11 << ", " << C_10 << endl;
    if (debug) cout << "ids is : " << com0 << ", " << com1 << ", " << com2 << ", " << com3 << ", " << com4 << ", " << com5 << endl;
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        if (flag != 2) continue;
        ic_tot++;
        int layerid = cluster->getlayerid();
        double _Q = cluster->getenergydeposit();
        // double t_phi=cluster->getrecphi_mTPC();
 
        double _phi = 99;
        double _Z = cluster->getRecZ();
        if (Switch_CCmTPC == 0) _phi = cluster->getrecphi();
        else if (Switch_CCmTPC == 1) _phi = cluster->getrecphi_mTPC();
        else if (Switch_CCmTPC == 2) {
            RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
            RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
            _phi = cluster_x->get_merge_phi();
        }
        if (debug)
            cout << "2D: phi_cc, phi_merge, index " << cluster->getrecphi() << ", " << cluster->get_merge_phi() << ", " << _loop
                << endl;
        int sheetid = cluster->getsheetid();
        if (_Q < MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
        
        // --- cuts Q>10fC, size>1
        //int clusterflage = cluster->getclusterflage();
        //RecCgemCluster* clusterV=*(recCgemClusterCol->begin()+clusterflage);
        //double QV=clusterV->getenergydeposit();
        //if(QV<10.) continue;
        //int size_V = 1+abs(clusterV->getclusterflage()-clusterV->getclusterflagb());
        //if(size_V<2) continue;
 
        if (fabs(C_00 - _phi) < 0.001 || fabs(C_01 - _phi) < 0.001 || fabs(C_11 - _phi) < 0.001 || fabs(C_10 - _phi) < 0.001 ||
                fabs(C_21 - _phi) < 0.001 || fabs(C_20 - _phi) < 0.001) { // is it dangerious to compare two float values ?
            if (layerid == 2 && sheetid == 1) {
                L_21.push_back(_loop);
                Q_21.push_back(_Q);
            }
            if (layerid == 2 && sheetid == 0) {
                L_20.push_back(_loop);
                Q_20.push_back(_Q);
            }
            if (layerid == 1 && sheetid == 1) {
                L_11.push_back(_loop);
                Q_11.push_back(_Q);
            }
            if (layerid == 1 && sheetid == 0) {
                L_10.push_back(_loop);
                Q_10.push_back(_Q);
            }
            if (layerid == 0 && sheetid == 0 && _phi > 0) {
                L_01.push_back(_loop);
                Q_01.push_back(_Q);
            }
            if (layerid == 0 && sheetid == 0 && _phi < 0) {
                L_00.push_back(_loop);
                Q_00.push_back(_Q);
            }
        }
    }
 
    if (debug) cout << "finish the third loop" << endl;
 
    Min_chi = 1E10;
 
    clst_21 = L_21.size();
    clst_11 = L_11.size();
    clst_01 = L_01.size();
    clst_20 = L_20.size();
    clst_10 = L_10.size();
    clst_00 = L_00.size();
 
    if (debug)
        cout << "clst_21, clst_11, clst_01, clst_20, clst_10, clst_00 : " << clst_21 << ", " << clst_11 << ", " << clst_01 << ", "
            << clst_20 << ", " << clst_10 << ", " << clst_00 << endl;
 
    L_00_s = GetNMaxQ(Q_00, L_00, Nmax);
    L_01_s = GetNMaxQ(Q_01, L_01, Nmax);
    L_10_s = GetNMaxQ(Q_10, L_10, Nmax);
    L_11_s = GetNMaxQ(Q_11, L_11, Nmax);
    L_20_s = GetNMaxQ(Q_20, L_20, Nmax);
    L_21_s = GetNMaxQ(Q_21, L_21, Nmax);
 
    if (TEST_N == 1) {
        L_21_s.clear();
        L_21_s.push_back(-1);
    } else if (TEST_N == 2) {
        L_11_s.clear();
        L_11_s.push_back(-1);
    } else if (TEST_N == 3) {
        L_01_s.clear();
        L_01_s.push_back(-1);
    } else if (TEST_N == 4) {
        L_00_s.clear();
        L_00_s.push_back(-1);
    } else if (TEST_N == 5) {
        L_10_s.clear();
        L_10_s.push_back(-1);
    } else if (TEST_N == 6) {
        L_20_s.clear();
        L_20_s.push_back(-1);
    }
 
    NVComb = L_00_s.size() * L_01_s.size() * L_10_s.size() * L_11_s.size() * L_20_s.size() * L_21_s.size();
 
    for (int i_21 = 0; i_21 < L_21_s.size(); i_21++) {
        for (int i_20 = 0; i_20 < L_20_s.size(); i_20++) {
            for (int i_11 = 0; i_11 < L_11_s.size(); i_11++) {
                for (int i_10 = 0; i_10 < L_10_s.size(); i_10++) {
                    for (int i_00 = 0; i_00 < L_00_s.size(); i_00++) {
                        for (int i_01 = 0; i_01 < L_01_s.size(); i_01++) {
 
                            Vec_layer.clear();
                            Vec_phi.clear();
                            Vec_Z.clear();
                            Vec_v.clear();
                            Vec_layerid.clear();
                            Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
                            Vec_flag.clear();
                            Vec_Q.clear();
                            Vec_sheetid.clear();
                            Vec_m_layer.clear();
                            Vec_m_phi.clear();
                            Vec_m_Z.clear();
                            Vec_m_v.clear();
                            Vec_m_layerid.clear();
                            Vec_m_flag.clear();
                            Vec_m_Q.clear();
                            Vec_m_sheetid.clear();
                            _loop = 0;
                            iter = recCgemClusterCol->begin();
                            for (; iter != recCgemClusterCol->end(); iter++) {
                                _loop++;
                                RecCgemCluster *cluster = (*iter);
                                int flag = cluster->getflag();
                                if (flag != 2) continue;
                                if (!(_loop == L_21_s[i_21] || _loop == L_20_s[i_20] || _loop == L_11_s[i_11] || _loop == L_10_s[i_10] ||
                                            _loop == L_00_s[i_00] || _loop == L_01_s[i_01]))
                                    continue;
 
                                int layerid = cluster->getlayerid();
                                int sheetid = cluster->getsheetid();
                                double _phi = cluster->getrecphi();
                                double _v = cluster->getrecv();
                                double _Z = cluster->getRecZ();
                                double _Q = cluster->getenergydeposit();
                                double t_phi = cluster->getrecphi_mTPC();
                                double t_v = cluster->getrecv_mTPC();
                                double t_Z = cluster->getRecZ_mTPC();
 
                                RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
                                RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
                                double merged_phi = cluster_x->get_merge_phi();
                                RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
                                double merged_v = cluster_v->get_merge_v();
 
                                int vir_layerid = GetVirLay(layerid, _phi);
 
                                Vec_flag.push_back(flag);
                                Vec_layerid.push_back(layerid);
                                Vec_sheetid.push_back(sheetid);
                                Vec_layer.push_back(R_layer[layerid]);
                                Vec_Virlayerid.push_back(vir_layerid);
 
                                Vec_m_phi.push_back(t_phi);
                                Vec_m_v.push_back(t_v);
                                Vec_m_Z.push_back(t_Z);
 
                                if (Switch_CCmTPC == 0) {
                                    Vec_phi.push_back(_phi);
                                    Vec_v.push_back(_v);
                                    Vec_Z.push_back(_Z);
                                } else if (Switch_CCmTPC == 1) {
                                    Vec_phi.push_back(t_phi);
                                    Vec_v.push_back(t_v);
                                    Vec_Z.push_back(t_Z);
                                } else if (Switch_CCmTPC == 2) {
 
                                    CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
                                    double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
                                    Vec_phi.push_back(merged_phi);
                                    Vec_v.push_back(merged_v);
                                    Vec_Z.push_back(merged_Z);
                                }
 
                            } // one combo
                            if (Vec_layerid.size() != 6 && TEST_N == 0) return false;
                            else if (Vec_layerid.size() != 5 && TEST_N > 0) return false;
                            OrderClusters();
 
                            // discard the combination if V strip is re-used !
                            if (Vec_v[0] == Vec_v[1] || Vec_v[0] == Vec_v[2] || Vec_v[0] == Vec_v[3] || Vec_v[0] == Vec_v[4] ||
                                    Vec_v[0] == Vec_v[5] || Vec_v[1] == Vec_v[2] || Vec_v[1] == Vec_v[3] || Vec_v[1] == Vec_v[4] ||
                                    Vec_v[1] == Vec_v[5] || Vec_v[2] == Vec_v[3] || Vec_v[2] == Vec_v[4] || Vec_v[2] == Vec_v[5] ||
                                    Vec_v[3] == Vec_v[4] || Vec_v[3] == Vec_v[5] || Vec_v[4] == Vec_v[5])
                                continue;
                            StraightLine *l_outer_tmp;
 
                            if (TEST_N == 1 || TEST_N == 6) l_outer_tmp = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
                            else if (TEST_N == 2 || TEST_N == 3 || TEST_N == 4 || TEST_N == 5)
                                l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
                            else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
 
                            TMinuit *_fit = Fit(l_outer_tmp, fa, 2);
                            double _chi = _fit->fAmin;
                            if (_chi < Min_chi) {
                                Min_chi = _chi;
                                iC_00 = L_00_s[i_00];
                                iC_01 = L_01_s[i_01];
                                iC_10 = L_10_s[i_10];
                                iC_11 = L_11_s[i_11];
                                iC_20 = L_20_s[i_20];
                                iC_21 = L_21_s[i_21];
                            }
                            delete _fit;
                            delete l_outer_tmp;
                        }
                    }
                }
            }
        }
    }
 
    Vec_layer.clear();
    Vec_phi.clear();
    Vec_Z.clear();
    Vec_v.clear();
    Vec_layerid.clear();
    Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
    Vec_flag.clear();
    Vec_Q.clear();
    Vec_sheetid.clear();
    Vec_xCluSize.clear();
    Vec_vCluSize.clear();
    Vec_XQ.clear();
    Vec_VQ.clear();
    Vec_XstripQ.clear();
    Vec_VstripQ.clear();
    Vec_XstripT.clear();
    Vec_VstripT.clear();
    Vec_XstripTf.clear();
    Vec_VstripTf.clear();
    Vec_XstripID.clear();
    Vec_VstripID.clear();
    Vec_m_layer.clear();
    Vec_m_phi.clear();
    Vec_m_Z.clear();
    Vec_m_v.clear();
    Vec_m_layerid.clear();
    Vec_m_flag.clear();
    Vec_m_Q.clear();
    Vec_m_sheetid.clear();
    if (debug) cout << "finish the forth loop" << endl;
    _loop = 0;
    ic = 0;
    iter = recCgemClusterCol->begin();
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        int layerid = cluster->getlayerid();
        int sheetid = cluster->getsheetid();
        double _phi = cluster->getrecphi();
        double _v = cluster->getrecv();
        double _Z = cluster->getRecZ();
        double _Q = cluster->getenergydeposit();
        double t_phi = cluster->getrecphi_mTPC();
        double t_v = cluster->getrecv_mTPC();
        double t_Z = cluster->getRecZ_mTPC();
 
        int clusterid = cluster->getclusterid();
        // int clusterflagb=cluster->getclusterflagb();
        // int clusterflage=cluster->getclusterflage();
        double x = R_layer[layerid] * cos(_phi);
        double y = R_layer[layerid] * sin(_phi);
        double z = cluster->getRecZ();
        int vir_layerid = GetVirLay(layerid, _phi);
        if (_loop > MAX_Clusters) break;
        if (flag != 2) continue;
 
        RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
        RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
        double merged_phi = cluster_x->get_merge_phi();
        RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
        double merged_v = cluster_v->get_merge_v();
 
        tr_id_cluster[ic] = clusterid;
        tr_x_cluster[ic] = x;
        tr_y_cluster[ic] = y;
        tr_z_cluster[ic] = z;
        tr_Q_cluster[ic] = _Q;
        tr_phi_cluster[ic] = _phi;
        tr_layer_cluster[ic] = layerid;
        tr_sheet_cluster[ic] = sheetid;
        tr_Is_selected_cluster[ic] = 0;
        // if ((_loop == iC_00 || _loop == iC_10 || _loop == iC_20 || _loop == iC_01 || _loop == iC_11 || _loop == iC_21))
        // {tr_Is_selected_cluster[ic] = 1;}
 
        // cout << "Test : ic = " << ic << "  tr_layer_cluster = " << tr_layer_cluster[ic] << "  tr_sheet_cluster = " <<
        // tr_sheet_cluster[ic] << "  tr_phi_cluster = " << tr_phi_cluster[ic] << "  tr_Is_selected_cluster = " <<
        // tr_Is_selected_cluster[ic] << endl;
        ic++;
        if (!(_loop == iC_00 || _loop == iC_10 || _loop == iC_20 || _loop == iC_01 || _loop == iC_11 || _loop == iC_21)) continue;
        tr_Is_selected_cluster[ic - 1] = 1;
        SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(), "/Event/Digi/CgemDigiCol");
        int clusterflagb = cluster->getclusterflagb();
        int clusterflage = cluster->getclusterflage();
        RecCgemClusterCol::iterator id_XCluster = recCgemClusterCol->begin() + clusterflagb;
        RecCgemCluster *xcluster = (*id_XCluster);
        int xclustersize = abs(xcluster->getclusterflage() - xcluster->getclusterflagb()) + 1;
        double X_Q = xcluster->getenergydeposit();
        RecCgemClusterCol::iterator id_VCluster = recCgemClusterCol->begin() + clusterflage;
        RecCgemCluster *vcluster = (*id_VCluster);
        int vclustersize = abs(vcluster->getclusterflage() - vcluster->getclusterflagb()) + 1;
        double V_Q = vcluster->getenergydeposit();
        int xclusterb = xcluster->getclusterflagb();
        int xclustere = xcluster->getclusterflage();
        int vclusterb = vcluster->getclusterflagb();
        int vclustere = vcluster->getclusterflage();
        int layerid_str, sheetid_str, stripid_str, time_chnnel;
        int nxstrips = 0;
        int nvstrips = 0;
        double energydeposit;
        double Q_fC, T_ns, Tf_ns;
        CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
        Vec_x_stripsQ.clear();
        Vec_v_stripsQ.clear();
        Vec_x_stripsT.clear();
        Vec_v_stripsT.clear();
        Vec_x_stripsTf.clear();
        Vec_v_stripsTf.clear();
        Vec_x_stripsid.clear();
        Vec_v_stripsid.clear();
        for (; iter_digi != cgemDigiCol->end(); iter_digi++) {
            layerid_str = CgemID::layer((*iter_digi)->identify());
            sheetid_str = CgemID::sheet((*iter_digi)->identify());
            stripid_str = CgemID::strip((*iter_digi)->identify());
            bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
            int flagxv;
            if (striptype == true) flagxv = 0;
            else flagxv = 1;
            if (layerid == layerid_str && sheetid == sheetid_str) {
                Q_fC = (*iter_digi)->getCharge_fc();
                T_ns = (*iter_digi)->getTime_ns();
                Tf_ns = get_Time(iter_digi);
                if (flagxv == 0) {
                    if (stripid_str >= xclusterb && stripid_str <= xclustere) {
                        // cout<<"stripid_str  =   "<<stripid_str<<"   xclusterb   =
                        // "<<xclusterb<<"  xclustere  =  "<<xclustere<<"  "<<flagxv<<endl;
                        if (xclustersize < 100) {
                            Vec_x_stripsid.push_back(stripid_str);
                            Vec_x_stripsQ.push_back(Q_fC);
                            Vec_x_stripsT.push_back(T_ns);
                            Vec_x_stripsTf.push_back(Tf_ns);
                            nxstrips++;
                        }
                    }
                }
                if (flagxv == 1) {
                    if (stripid_str >= vclusterb && stripid_str <= vclustere) {
                        if (vclustersize < 100) {
                            Vec_v_stripsid.push_back(stripid_str);
                            Vec_v_stripsQ.push_back(Q_fC);
                            Vec_v_stripsT.push_back(T_ns);
                            Vec_v_stripsTf.push_back(Tf_ns);
                            nvstrips++;
                        }
                    }
                }
            }
        }
 
        Vec_flag.push_back(flag);
        Vec_layerid.push_back(layerid);
        Vec_Virlayerid.push_back(vir_layerid);
        Vec_sheetid.push_back(sheetid);
        Vec_layer.push_back(R_layer[layerid]);
        Vec_xCluSize.push_back(xclustersize);
        Vec_vCluSize.push_back(vclustersize);
        Vec_XQ.push_back(X_Q);
        Vec_VQ.push_back(V_Q);
        Vec_XstripQ.push_back(Vec_x_stripsQ);
        Vec_VstripQ.push_back(Vec_v_stripsQ);
        Vec_XstripT.push_back(Vec_x_stripsT);
        Vec_VstripT.push_back(Vec_v_stripsT);
        Vec_XstripTf.push_back(Vec_x_stripsTf);
        Vec_VstripTf.push_back(Vec_v_stripsTf);
        Vec_XstripID.push_back(Vec_x_stripsid);
        Vec_VstripID.push_back(Vec_v_stripsid);
 
        Vec_m_phi.push_back(t_phi);
        Vec_m_v.push_back(t_v);
        Vec_m_Z.push_back(t_Z);
 
        if (Switch_CCmTPC == 0) {
            Vec_phi.push_back(_phi);
            Vec_v.push_back(_v);
            Vec_Z.push_back(_Z);
        } else if (Switch_CCmTPC == 1) {
            Vec_phi.push_back(t_phi);
            Vec_v.push_back(t_v);
            Vec_Z.push_back(t_Z);
        } else if (Switch_CCmTPC == 2) {
 
            CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
            double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
            Vec_phi.push_back(merged_phi);
            Vec_v.push_back(merged_v);
            Vec_Z.push_back(merged_Z);
        }
 
        cluster->setTrkId(0);
        vecclusters.push_back(cluster);
 
    } // one combo
    NC = vecclusters.size();
    tr_ncluster = ic;
    if (debug) cout << "finish the fifth loop" << endl;
 
    for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
        if (Vec_layerid[i_cl] == 2) _clst_2++;
        if (Vec_layerid[i_cl] == 1) _clst_1++;
        if (Vec_layerid[i_cl] == 0) _clst_0++;
    }
    if (debug) cout << "Vec_layerid.size is " << Vec_layerid.size() << endl;
 
    if (Vec_layerid.size() != 6 && TEST_N == 0) return false;
    else if (Vec_layerid.size() != 5 && TEST_N > 0) return false;
    // OrderClusters();
    OrderClusterSizeQ();
 
    if (TEST_N == 1 || TEST_N == 6) l_outer = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
    else if (TEST_N == 2 || TEST_N == 3 || TEST_N == 4 || TEST_N == 5)
        l_outer = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
    else if (TEST_N == 0) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
    if (debug) {
        cout << "Event " << event << ", Clu. ID are " << iC_00 << ", " << iC_01 << ", " << iC_10 << ", " << iC_11 << ", " << iC_20
            << ", " << iC_21 << ",  Loop_MaxQ Min_chi = " << Min_chi << endl;
    }
 
    if (Min_chi < Chisq_cut) return true;
    else return false;
}

Referenced by execute().

MC_truth()

vector< double > CgemLineFit::MC_truth

(

)

Definition at line 3136 of file CgemLineFit.cxx.

                                     {
    int nTruth[3] = {0, 0, 0};
    int trkID_Layer1[100], trkID_Layer2[100], trkID_Layer3[100];
    int pdg_Layer1[100], pdg_Layer2[100], pdg_Layer3[100];
    double x_pre_Layer1[100], x_pre_Layer2[100], x_pre_Layer3[100];
    double y_pre_Layer1[100], y_pre_Layer2[100], y_pre_Layer3[100];
    double z_pre_Layer1[100], z_pre_Layer2[100], z_pre_Layer3[100];
    double x_post_Layer1[100], x_post_Layer2[100], x_post_Layer3[100];
    double y_post_Layer1[100], y_post_Layer2[100], y_post_Layer3[100];
    double z_post_Layer1[100], z_post_Layer2[100], z_post_Layer3[100];
    double phi_Layer1[100], phi_Layer2[100], phi_Layer3[100];
    double z_Layer1[100], z_Layer2[100], z_Layer3[100];
    double V_Layer1[100], V_Layer2[100], V_Layer3[100];
    double px_pre_Layer1[100], px_pre_Layer2[100], px_pre_Layer3[100];
    double py_pre_Layer1[100], py_pre_Layer2[100], py_pre_Layer3[100];
    double pz_pre_Layer1[100], pz_pre_Layer2[100], pz_pre_Layer3[100];
    double en_Layer1[100], en_Layer2[100], en_Layer3[100];
    SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
    /*  if (!cgemMcHitCol)
        {
        std::cout << "Could not retrieve cgemMcHitCol" << std::endl;
        return;
        } */
    Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
    for (; iter_truth != cgemMcHitCol->end(); iter_truth++) {
        int iLayer = (*iter_truth)->GetLayerID();
        int trkID = (*iter_truth)->GetTrackID();
        int pdg = (*iter_truth)->GetPDGCode();
        double x_pre = (*iter_truth)->GetPositionXOfPrePoint(); // in mm
        double y_pre = (*iter_truth)->GetPositionYOfPrePoint();
        double z_pre = (*iter_truth)->GetPositionZOfPrePoint();
        double x_post = (*iter_truth)->GetPositionXOfPostPoint();
        double y_post = (*iter_truth)->GetPositionYOfPostPoint();
        double z_post = (*iter_truth)->GetPositionZOfPostPoint();
        double x_middle = 0.5 * (x_pre + x_post);
        double y_middle = 0.5 * (y_pre + y_post);
        double z_middle = 0.5 * (z_pre + z_post);
        double r_middle = sqrt(x_middle * x_middle + y_middle * y_middle); // cout<<"r_middle = "<<r_middle<<endl;
        double phi_middle = atan2(y_middle, x_middle);
        double px_pre = (*iter_truth)->GetMomentumXOfPrePoint();
        double py_pre = (*iter_truth)->GetMomentumYOfPrePoint();
        double pz_pre = (*iter_truth)->GetMomentumZOfPrePoint();
        double en = (*iter_truth)->GetTotalEnergyDeposit();
        CgemGeoReadoutPlane *readoutPlane = NULL;
        // cout<<"phi_middle, z_middle = "<<phi_middle<<", "<<z_middle<<endl;
        // cout<<"iLayer = "<<iLayer<<endl;
        for (int j = 0; j < myNSheets[iLayer]; j++) {
            // cout<<"j = "<<j <<endl;
            readoutPlane = m_cgemGeomSvc->getReadoutPlane(iLayer, j);
            // cout<<"OnThePlane: "<<readoutPlane->OnThePlane(phi_middle,z_middle)<<endl;
            if (readoutPlane->OnThePlane(phi_middle, z_middle)) break;
            readoutPlane = NULL;
        }
        double V_mid = 9999;
        if (readoutPlane == NULL) {
            cout << "CgemClusterCreate::fillMCTruth() readoutPlane not found with phi_middle = " << phi_middle << ", layer = " << iLayer
                << endl;
        } else {
            // readoutPlane->print();
            // cout<<"phi_middle, z_middle = "<<phi_middle<<", "<<z_middle<<endl;
            V_mid = readoutPlane->getVFromPhiZ(phi_middle, z_middle);
            // cout<<"V_mid = "<<V_mid<<endl;
            // double V_mid = readoutPlane->GetVFromLocalXZ(phi_middle,z_middle);
        }
        switch (iLayer) {
            case 0:
                if (nTruth[0] >= 100) break;
                trkID_Layer1[nTruth[0]] = trkID;
                pdg_Layer1[nTruth[0]] = pdg;
                x_pre_Layer1[nTruth[0]] = x_pre;
                y_pre_Layer1[nTruth[0]] = y_pre;
                z_pre_Layer1[nTruth[0]] = z_pre;
                x_post_Layer1[nTruth[0]] = x_post;
                y_post_Layer1[nTruth[0]] = y_post;
                z_post_Layer1[nTruth[0]] = z_post;
                phi_Layer1[nTruth[0]] = atan2(y_post + y_pre, x_post + x_pre);
                z_Layer1[nTruth[0]] = z_middle;
                V_Layer1[nTruth[0]] = V_mid;
                if (atan2(y_post + y_pre, x_post + x_pre) > 0) {
                    x_0_up_mc2 = atan2(y_post + y_pre, x_post + x_pre);
                    z_0_up_mc2 = z_middle;
                    v_0_up_mc2 = V_mid;
                } else {
                    x_0_down_mc2 = atan2(y_post + y_pre, x_post + x_pre);
                    z_0_down_mc2 = z_middle;
                    v_0_down_mc2 = V_mid;
                }
                px_pre_Layer1[nTruth[0]] = px_pre;
                py_pre_Layer1[nTruth[0]] = py_pre;
                pz_pre_Layer1[nTruth[0]] = pz_pre;
                en_Layer1[nTruth[0]] = en;
                break;
            case 1:
                if (nTruth[1] >= 100) break;
                trkID_Layer2[nTruth[1]] = trkID;
                pdg_Layer2[nTruth[1]] = pdg;
                x_pre_Layer2[nTruth[1]] = x_pre;
                y_pre_Layer2[nTruth[1]] = y_pre;
                z_pre_Layer2[nTruth[1]] = z_pre;
                x_post_Layer2[nTruth[1]] = x_post;
                y_post_Layer2[nTruth[1]] = y_post;
                z_post_Layer2[nTruth[1]] = z_post;
                phi_Layer2[nTruth[1]] = atan2(y_post + y_pre, x_post + x_pre);
                z_Layer2[nTruth[1]] = z_middle;
                V_Layer2[nTruth[1]] = V_mid;
                if (atan2(y_post + y_pre, x_post + x_pre) > 0) {
                    x_1_up_mc2 = atan2(y_post + y_pre, x_post + x_pre);
                    z_1_up_mc2 = z_middle;
                    v_1_up_mc2 = V_mid;
                } else {
                    x_1_down_mc2 = atan2(y_post + y_pre, x_post + x_pre);
                    z_1_down_mc2 = z_middle;
                    v_1_down_mc2 = V_mid;
                }
                px_pre_Layer2[nTruth[1]] = px_pre;
                py_pre_Layer2[nTruth[1]] = py_pre;
                pz_pre_Layer2[nTruth[1]] = pz_pre;
                en_Layer2[nTruth[1]] = en;
                break;
            case 2:
                if (nTruth[2] >= 100) break;
                trkID_Layer3[nTruth[2]] = trkID;
                pdg_Layer3[nTruth[2]] = pdg;
                x_pre_Layer3[nTruth[2]] = x_pre;
                y_pre_Layer3[nTruth[2]] = y_pre;
                z_pre_Layer3[nTruth[2]] = z_pre;
                x_post_Layer3[nTruth[2]] = x_post;
                y_post_Layer3[nTruth[2]] = y_post;
                z_post_Layer3[nTruth[2]] = z_post;
                phi_Layer3[nTruth[2]] = atan2(y_post + y_pre, x_post + x_pre);
                z_Layer3[nTruth[2]] = z_middle;
                V_Layer3[nTruth[2]] = V_mid;
                if (atan2(y_post + y_pre, x_post + x_pre) > 0) {
                    x_2_up_mc2 = atan2(y_post + y_pre, x_post + x_pre);
                    v_2_up_mc2 = V_mid;
                    z_2_up_mc2 = z_middle;
                } else {
                    x_2_down_mc2 = atan2(y_post + y_pre, x_post + x_pre);
                    v_2_down_mc2 = V_mid;
                    z_2_down_mc2 = z_middle;
                }
                px_pre_Layer3[nTruth[2]] = px_pre;
                py_pre_Layer3[nTruth[2]] = py_pre;
                pz_pre_Layer3[nTruth[2]] = pz_pre;
                en_Layer3[nTruth[2]] = en;
                break;
            default: cout << "wrong layer ID for CGEM: " << iLayer << endl;
        }
        nTruth[iLayer]++;
    } // loop hit truth
    for (int i = 0; i < 3; i++)
        if (nTruth[i] > 100) nTruth[i] = 100;
    SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
 
    Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
    HepLorentzVector p4_mu(0, 0, 0, 0);
    HepLorentzVector p4_pos(0, 0, 0, 0);
    for (; iter_mc != mcParticleCol->end(); iter_mc++) {
        if (!(*iter_mc)->decayFromGenerator()) continue;
        HepLorentzVector p4_mc(0, 0, 0, 0);
 
        if (13 == abs((*iter_mc)->particleProperty())) {
 
            p4_mu = (*iter_mc)->initialFourMomentum();
            p4_pos = (*iter_mc)->initialPosition();
        }
    }
    vector<double> _mc;
    Hep3Vector _BP(p4_pos.x() * 10, p4_pos.y() * 10, p4_pos.z() * 10);
    pos_x = _BP[0];
    pos_y = _BP[1];
    pos_z = _BP[2];
    _mc = Get4Par(p4_mu, _BP);
 
    _mc.push_back(p4_mu.px());
    _mc.push_back(p4_mu.py());
    _mc.push_back(p4_mu.pz());
    if (HitPos(p4_mu, _BP)) _mc.push_back(1);
    else _mc.push_back(0);
 
    return _mc;
}

Referenced by Get_MCInf().

Min_diff()

double CgemLineFit::Min_diff ( double arg1, double arg2 )

static

Definition at line 3955 of file CgemLineFit.cxx.

                                                     {
    to_0_2pi(arg1);
    to_0_2pi(arg2);
    double diff_raw = fabs(arg1 - arg2);
 
    if (diff_raw > acos(-1)) return (2 * acos(-1) - diff_raw);
    else return diff_raw;
}

Referenced by dV(), dx(), Min_dis(), Min_dis_down(), and Min_dis_up().

Min_dis()

double CgemLineFit::Min_dis	(	int	R_id,
		double	x,
		double	z,
		StraightLine *	l1 )

Definition at line 3964 of file CgemLineFit.cxx.

                                                                          {
    // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
    double phiV_up[2], phiV_down[2];
    HepPoint3D Up, Down;
 
    int vir_R_id = GetVirLay(R_id, x);
 
    if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
    else Mp->getPointAligned(vir_R_id, l1, Up, Down, phiV_up, phiV_down);
 
    double ds1 = pow((phiV_up[1] - v), 2) + pow(R_layer[R_id] * Min_diff(phiV_up[0], x), 2);
 
    double ds2 = pow((phiV_down[1] - v), 2) + pow(R_layer[R_id] * Min_diff(phiV_down[0], x), 2);
 
    return min(ds1, ds2);
}

Referenced by Filter().

Min_dis_down()

double CgemLineFit::Min_dis_down	(	int	R_id,
		double	x,
		double	z,
		StraightLine *	l1 )

Definition at line 3993 of file CgemLineFit.cxx.

                                                                               {
 
    // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
    double phiV_up[2], phiV_down[2];
    HepPoint3D Up, Down;
    int vir_R_id = GetVirLay(R_id, x);
    if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
    else Mp->getPointAligned(vir_R_id, l1, Up, Down, phiV_up, phiV_down);
 
    double ds2 = pow((phiV_down[1] - v), 2) + pow(R_layer[R_id] * Min_diff(phiV_down[0], x), 2);
 
    return ds2;
}

Referenced by Filter().

Min_dis_up()

double CgemLineFit::Min_dis_up	(	int	R_id,
		double	x,
		double	z,
		StraightLine *	l1 )

Definition at line 3981 of file CgemLineFit.cxx.

                                                                             {
    // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
    double phiV_up[2], phiV_down[2];
    HepPoint3D Up, Down;
    int vir_R_id = GetVirLay(R_id, x);
    if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
    else Mp->getPointAligned(vir_R_id, l1, Up, Down, phiV_up, phiV_down);
    double ds1 = pow((phiV_up[1] - z), 2) + pow(R_layer[R_id] * Min_diff(phiV_up[0], x), 2);
 
    return ds1;
}

Referenced by Filter().

OrderClusters()

void CgemLineFit::OrderClusters

(

)

Definition at line 2834 of file CgemLineFit.cxx.

                                {
    int noswap(0);
    while (noswap != Vec_flag.size() - 1) {
        noswap = 0;
        for (int i = 0; i < Vec_flag.size() - 1; i++) {
            //cout << "In OrderClusters()  i = " << i << "  Vec_flag.size() = " << Vec_flag.size() << endl;
            if (Vec_flag[i] < Vec_flag[i + 1]) {
                Swap(i, i + 1);
                //cout << "layerid = " << Vec_layerid[i] << "  sheetid = " << Vec_sheetid[i] << "  flag = " << Vec_flag[i] << "  phi = "
                //<< Vec_phi[i] << endl;
            } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] < Vec_layerid[i + 1]) {
                Swap(i, i + 1);
                //cout << "layerid = " << Vec_layerid[i] << "  sheetid = " << Vec_sheetid[i] << "  flag = " << Vec_flag[i] << "  phi = "
                //<< Vec_phi[i] << endl;
            } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] == Vec_layerid[i + 1] && (Vec_phi[i] < Vec_phi[i + 1])) {
                Swap(i, i + 1);
                //cout << "layerid = " << Vec_layerid[i] << "  sheetid = " << Vec_sheetid[i] << "  flag = " << Vec_flag[i] << "  phi = "
                //<< Vec_phi[i] << endl;
            } else {
                //cout << "xxxxxxx" << endl;
                //cout << "(i)layerid = " << Vec_layerid[i] << "  sheetid = " << Vec_sheetid[i] << "Virlayerid = " << Vec_Virlayerid[i]
                //<< "  flag = " << Vec_flag[i] << "  phi = " << Vec_phi[i] << endl; cout << "(i+1)layerid = " << Vec_layerid[i+1] << "sheetid = " << Vec_sheetid[i+1] << "Virlayerid = " << Vec_Virlayerid[i+1] << "  flag = " << Vec_flag[i+1] << "  phi = "
                //<< Vec_phi[i+1] << endl;
                noswap++;
            }
        }
    }
}

Referenced by Data_Closest(), Data_Max(), Loop_All(), Loop_MaxQ(), and ToyMC().

OrderClusterSizeQ()

void CgemLineFit::OrderClusterSizeQ

(

)

Definition at line 2880 of file CgemLineFit.cxx.

                                    {
    int noswap(0);
    while (noswap != Vec_flag.size() - 1) {
        noswap = 0;
        for (int i = 0; i < Vec_flag.size() - 1; i++) {
            if (Vec_flag[i] < Vec_flag[i + 1]) {
                Swap(i, i + 1);
                swap(Vec_xCluSize[i], Vec_xCluSize[i + 1]);
                swap(Vec_vCluSize[i], Vec_vCluSize[i + 1]);
                swap(Vec_XQ[i], Vec_XQ[i + 1]);
                swap(Vec_VQ[i], Vec_VQ[i + 1]);
                swap(Vec_XstripQ[i], Vec_XstripQ[i + 1]);
                swap(Vec_VstripQ[i], Vec_VstripQ[i + 1]);
                swap(Vec_XstripT[i], Vec_XstripT[i + 1]);
                swap(Vec_VstripT[i], Vec_VstripT[i + 1]);
                swap(Vec_XstripTf[i], Vec_XstripTf[i + 1]);
                swap(Vec_VstripTf[i], Vec_VstripTf[i + 1]);
                swap(Vec_XstripID[i], Vec_XstripID[i + 1]);
                swap(Vec_VstripID[i], Vec_VstripID[i + 1]);
            } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] < Vec_layerid[i + 1]) {
                Swap(i, i + 1);
                swap(Vec_xCluSize[i], Vec_xCluSize[i + 1]);
                swap(Vec_vCluSize[i], Vec_vCluSize[i + 1]);
                swap(Vec_XQ[i], Vec_XQ[i + 1]);
                swap(Vec_VQ[i], Vec_VQ[i + 1]);
                swap(Vec_XstripQ[i], Vec_XstripQ[i + 1]);
                swap(Vec_VstripQ[i], Vec_VstripQ[i + 1]);
                swap(Vec_XstripT[i], Vec_XstripT[i + 1]);
                swap(Vec_VstripT[i], Vec_VstripT[i + 1]);
                swap(Vec_XstripTf[i], Vec_XstripTf[i + 1]);
                swap(Vec_VstripTf[i], Vec_VstripTf[i + 1]);
                swap(Vec_XstripID[i], Vec_XstripID[i + 1]);
                swap(Vec_VstripID[i], Vec_VstripID[i + 1]);
            } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] == Vec_layerid[i + 1] && (Vec_phi[i] < Vec_phi[i + 1])) {
                Swap(i, i + 1);
                swap(Vec_xCluSize[i], Vec_xCluSize[i + 1]);
                swap(Vec_vCluSize[i], Vec_vCluSize[i + 1]);
                swap(Vec_XQ[i], Vec_XQ[i + 1]);
                swap(Vec_VQ[i], Vec_VQ[i + 1]);
                swap(Vec_XstripQ[i], Vec_XstripQ[i + 1]);
                swap(Vec_VstripQ[i], Vec_VstripQ[i + 1]);
                swap(Vec_XstripT[i], Vec_XstripT[i + 1]);
                swap(Vec_VstripT[i], Vec_VstripT[i + 1]);
                swap(Vec_XstripTf[i], Vec_XstripTf[i + 1]);
                swap(Vec_VstripTf[i], Vec_VstripTf[i + 1]);
                swap(Vec_XstripID[i], Vec_XstripID[i + 1]);
                swap(Vec_VstripID[i], Vec_VstripID[i + 1]);
            } else {
                noswap++;
            }
        }
    }
}

Referenced by Loop_MaxQ().

RealPhi()

double CgemLineFit::RealPhi

(

double

SimuPhi

)

Definition at line 4394 of file CgemLineFit.cxx.

                                          {
    double RPhi;
    if (SimuPhi <= TMath::Pi()) RPhi = SimuPhi;
    else RPhi = SimuPhi - 2 * TMath::Pi();
    return RPhi;
}

Referenced by ToyMC().

Rec_ClusterQ()

void CgemLineFit::Rec_ClusterQ

(

)

Definition at line 3528 of file CgemLineFit.cxx.

                               {
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5) {
            x_2_up_Q = Vec_XQ[i];
            v_2_up_Q = Vec_VQ[i];
        }
        if (Vec_Virlayerid[i] == 4) {
            x_2_down_Q = Vec_XQ[i];
            v_2_down_Q = Vec_VQ[i];
        }
        if (Vec_Virlayerid[i] == 3) {
            x_1_up_Q = Vec_XQ[i];
            v_1_up_Q = Vec_VQ[i];
        }
        if (Vec_Virlayerid[i] == 2) {
            x_1_down_Q = Vec_XQ[i];
            v_1_down_Q = Vec_VQ[i];
        }
        if (Vec_Virlayerid[i] == 1) {
            x_0_up_Q = Vec_XQ[i];
            v_0_up_Q = Vec_VQ[i];
        }
        if (Vec_Virlayerid[i] == 0) {
            x_0_down_Q = Vec_XQ[i];
            v_0_down_Q = Vec_VQ[i];
        }
    }
}

Referenced by execute().

Rec_Clusters()

void CgemLineFit::Rec_Clusters

(

)

Definition at line 3422 of file CgemLineFit.cxx.

                               {
 
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
 
        if (Vec_Virlayerid[i] == 5) {
            x_2_up = Vec_phi[i];
            v_2_up = Vec_v[i];
            z_2_up = Vec_Z[i];
        }
        if (Vec_Virlayerid[i] == 4) {
            x_2_down = Vec_phi[i];
            v_2_down = Vec_v[i];
            z_2_down = Vec_Z[i];
        }
        if (Vec_Virlayerid[i] == 3) {
            x_1_up = Vec_phi[i];
            v_1_up = Vec_v[i];
            z_1_up = Vec_Z[i];
        }
 
        if (Vec_Virlayerid[i] == 2) {
            x_1_down = Vec_phi[i];
            v_1_down = Vec_v[i];
            z_1_down = Vec_Z[i];
        }
 
        if (Vec_Virlayerid[i] == 1) {
            x_0_up = Vec_phi[i];
            v_0_up = Vec_v[i];
            z_0_up = Vec_Z[i];
        }
 
        if (Vec_Virlayerid[i] == 0) {
            x_0_down = Vec_phi[i];
            v_0_down = Vec_v[i];
            z_0_down = Vec_Z[i];
        }
    }
}

Referenced by execute().

Rec_Clusters_mTPC()

void CgemLineFit::Rec_Clusters_mTPC

(

)

Definition at line 3462 of file CgemLineFit.cxx.

                                    {
 
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5) {
            x_2_up_m = Vec_m_phi[i];
            v_2_up_m = Vec_m_v[i];
            z_2_up_m = Vec_m_Z[i];
        }
        if (Vec_Virlayerid[i] == 4) {
            x_2_down_m = Vec_m_phi[i];
            v_2_down_m = Vec_m_v[i];
            z_2_down_m = Vec_m_Z[i];
        }
        if (Vec_Virlayerid[i] == 3) {
            x_1_up_m = Vec_m_phi[i];
            v_1_up_m = Vec_m_v[i];
            z_1_up_m = Vec_m_Z[i];
        }
 
        if (Vec_Virlayerid[i] == 2) {
            x_1_down_m = Vec_m_phi[i];
            v_1_down_m = Vec_m_v[i];
            z_1_down_m = Vec_m_Z[i];
        }
 
        if (Vec_Virlayerid[i] == 1) {
            x_0_up_m = Vec_m_phi[i];
            v_0_up_m = Vec_m_v[i];
            z_0_up_m = Vec_m_Z[i];
        }
        if (Vec_Virlayerid[i] == 0) {
            x_0_down_m = Vec_m_phi[i];
            v_0_down_m = Vec_m_v[i];
            z_0_down_m = Vec_m_Z[i];
        }
    }
}

Referenced by execute().

Rec_ClusterSize()

void CgemLineFit::Rec_ClusterSize

(

)

Definition at line 3499 of file CgemLineFit.cxx.

                                  {
    for (int i = 0; i < Vec_Virlayerid.size(); i++) {
        if (Vec_Virlayerid[i] == 5) {
            x_2_up_size = Vec_xCluSize[i];
            v_2_up_size = Vec_vCluSize[i];
        }
        if (Vec_Virlayerid[i] == 4) {
            x_2_down_size = Vec_xCluSize[i];
            v_2_down_size = Vec_vCluSize[i];
        }
        if (Vec_Virlayerid[i] == 3) {
            x_1_up_size = Vec_xCluSize[i];
            v_1_up_size = Vec_vCluSize[i];
        }
        if (Vec_Virlayerid[i] == 2) {
            x_1_down_size = Vec_xCluSize[i];
            v_1_down_size = Vec_vCluSize[i];
        }
        if (Vec_Virlayerid[i] == 1) {
            x_0_up_size = Vec_xCluSize[i];
            v_0_up_size = Vec_vCluSize[i];
        }
        if (Vec_Virlayerid[i] == 0) {
            x_0_down_size = Vec_xCluSize[i];
            v_0_down_size = Vec_vCluSize[i];
        }
    }
}

Referenced by execute().

registerTrack()

StatusCode CgemLineFit::registerTrack	(	RecMdcTrackCol *&	recMdcTrackCol,
		RecMdcHitCol *&	recMdcHitCol )

Definition at line 2574 of file CgemLineFit.cxx.

                                                                                                   {
    MsgStream log(msgSvc(), name());
    StatusCode sc;
    IDataProviderSvc *eventSvc = NULL;
    Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
    if (eventSvc) {
        log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
    } else {
        log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
        return StatusCode::FAILURE;
        // return StatusCode::SUCCESS;
    }
    IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc *>(eventSvc);
    ;
 
    // --- register RecMdcTrack
    recMdcTrackCol = new RecMdcTrackCol;
    DataObject *aRecMdcTrackCol;
    eventSvc->findObject("/Event/Recon/RecMdcTrackCol", aRecMdcTrackCol);
    if (aRecMdcTrackCol != NULL) {
        dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
        eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
    }
 
    sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
    if (sc.isFailure()) {
        log << MSG::FATAL << " Could not register RecMdcTrack collection" << endreq;
        return StatusCode::FAILURE;
    }
    log << MSG::INFO << "RecMdcTrackCol registered successfully!" << endreq;
 
    recMdcHitCol = new RecMdcHitCol;
    DataObject *aRecMdcHitCol;
    eventSvc->findObject("/Event/Recon/RecMdcHitCol", aRecMdcHitCol);
    if (aRecMdcHitCol != NULL) {
        dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
        eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
    }
 
    sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
    if (sc.isFailure()) {
        log << MSG::FATAL << " Could not register RecMdcHit collection" << endreq;
        return StatusCode::FAILURE;
    }
    log << MSG::INFO << "RecMdcHitCol registered successfully!" << endreq;
 
    return StatusCode::SUCCESS;
} // end of stroeTracks

Referenced by execute().

Store_Trk()

void CgemLineFit::Store_Trk	(	TMinuit *	fit,
		int	trkid )

Definition at line 4088 of file CgemLineFit.cxx.

                                                   {
    double trkpar[4] = {-999, -999, -999, -999};
    double trkparErr[4] = {-999, -999, -999, -999};
 
    Int_t ierflg, npari, nparx, istat3;
    Double_t fmin, edm, errdef;
 
    fit->mnstat(fmin, edm, errdef, npari, nparx, istat3);
    for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
    double emat[16];
    fit->mnemat(emat, 4);
 
    RecMdcTrack *recMdcTrack = new RecMdcTrack();
 
    double helixPar[5];
    helixPar[0] = trkpar[0];
    helixPar[1] = trkpar[1];
    helixPar[2] = 0.0;
    helixPar[3] = trkpar[2];
    helixPar[4] = trkpar[3];
    double errorMat[15];
    int k = 0;
    errorMat[0] = emat[0];
    errorMat[1] = emat[1];
    errorMat[2] = 0;
    errorMat[3] = emat[2];
    errorMat[4] = emat[3];
    errorMat[5] = emat[5];
    errorMat[6] = 0;
    errorMat[7] = emat[6];
    errorMat[8] = emat[7];
    errorMat[9] = 0;
    errorMat[10] = 0;
    errorMat[11] = 0;
    errorMat[12] = emat[10];
    errorMat[13] = emat[11];
    errorMat[14] = emat[15];
 
    recMdcTrack->setChi2(fit->fAmin);
    recMdcTrack->setError(errorMat);
    recMdcTrack->setHelix(helixPar);
    recMdcTrack->setTrackId(trkid);
    /////
    //SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
    //RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();
    /////
    if (trkid == 0) {
        //cout << "before setNcluster : " <<  recMdcTrack->getNcluster() << endl;
        recMdcTrack->setNcluster(NC);
        //cout << "set NC = " << NC << endl;
        //cout << "after setNcluster : " <<  recMdcTrack->getNcluster() << endl;
        //cout << "before setVecClusters : " << recMdcTrack->getNcluster() << endl;
        recMdcTrack->setVecClusters(vecclusters);
        //cout << "after setVecClusters : " << recMdcTrack->getNcluster() << endl;
        m_trackList_tds->push_back(recMdcTrack);
        /* SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
           RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();
           cout << "after m_trackList_tds->push_back() : " << (*iter)->getNcluster() << endl; */
    }
}

Referenced by execute(), and GetIntersection().

Swap()

void CgemLineFit::Swap	(	int	i,
		int	j )

Definition at line 2863 of file CgemLineFit.cxx.

                                   {
 
    swap(Vec_layer[i], Vec_layer[j]);
    swap(Vec_sheetid[i], Vec_sheetid[j]);
    swap(Vec_phi[i], Vec_phi[j]);
    swap(Vec_Z[i], Vec_Z[j]);
    swap(Vec_v[i], Vec_v[j]);
    swap(Vec_layerid[i], Vec_layerid[j]);
    swap(Vec_Virlayerid[i], Vec_Virlayerid[j]);
    swap(Vec_flag[i], Vec_flag[j]);
    if (Vec_m_phi.size() == Vec_flag.size() && (Method == 1 || Method == 2)) {
        swap(Vec_m_phi[i], Vec_m_phi[j]);
        swap(Vec_m_Z[i], Vec_m_Z[j]);
        swap(Vec_m_v[i], Vec_m_v[j]);
    }
}

Referenced by OrderClusters(), and OrderClusterSizeQ().

to_0_2pi()

void CgemLineFit::to_0_2pi ( double & arg )

static

Definition at line 3950 of file CgemLineFit.cxx.

                                      {
    while (arg < -1 * acos(-1)) arg += acos(-1) * 2;
    while (arg > acos(-1)) arg -= acos(-1) * 2;
}

Referenced by Fit(), Get_Clusters(), GetIntersection(), and Min_diff().

ToyMC()

bool CgemLineFit::ToyMC

(

)

Definition at line 795 of file CgemLineFit.cxx.

                        {
 
    double cl_00(0), cl_01(0), cl_11(0), cl_10(0);
    int _loop(0);
    if (!Get_MCInf()) return false;
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
    RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
    for (; iter != recCgemClusterCol->end(); iter++) {
        _loop++;
        RecCgemCluster *cluster = (*iter);
        int flag = cluster->getflag();
        if (flag != 2) continue;
 
        int layerid = cluster->getlayerid();
        int sheetid = cluster->getsheetid();
        double _phi = cluster->getrecphi();
        if (MC) _phi = RealPhi(_phi);
        double _v = cluster->getrecv();
        double _Z = cluster->getRecZ();
        double _Q = cluster->getenergydeposit();
        int vir_layerid = GetVirLay(layerid, _phi);
        if (layerid == 2 || ((!Run10_Flag) && layerid == 0 && _phi > acos(-1))) continue;
 
        Vec_flag.push_back(flag);
        Vec_layerid.push_back(layerid);
        Vec_Virlayerid.push_back(vir_layerid);
        Vec_sheetid.push_back(sheetid);
        Vec_layer.push_back(R_layer[layerid]);
 
        Vec_phi.push_back(_phi);
        Vec_v.push_back(_v);
        Vec_Z.push_back(_Z);
 
        Vec_Q.push_back(_Q);
        cluster->setTrkId(0);
        vecclusters.push_back(cluster);
    }
    int NC = vecclusters.size();
    if (debug) cout << "NC is " << NC << endl;
 
    if (Vec_layer.size() < 3) { return false; }
 
    for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
        if (Vec_layerid[i_cl] == 2) _clst_2++;
        if (Vec_layerid[i_cl] == 1) _clst_1++;
        if (Vec_layerid[i_cl] == 0) _clst_0++;
    }
 
    if (debug) cout << "clst1, clst0 are " << _clst_1 << ", " << _clst_0 << endl;
    OrderClusters();
    if (_clst_1 == 2) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2); // should be this
    else return false;
 
    if (_clst_0 > 2) Filter(0, l_outer);
 
    return true;
}

Referenced by execute().

Trans()

vector< double > CgemLineFit::Trans ( double par0, double par1, double par2, double par3 )

static

Definition at line 4019 of file CgemLineFit.cxx.

                                                                                    {
    vector<double> rotat;
    rotat.clear();
 
    rotat.push_back(par0);
    rotat.push_back(par1);
    rotat.push_back(par2);
    rotat.push_back(par3);
    return rotat;
}

Referenced by execute(), Get_MCInf(), and Get_OtherIniPar().

Member Data Documentation

debug

bool CgemLineFit::debug

Definition at line 63 of file CgemLineFit.h.

Referenced by CgemLineFit(), Data_Max(), execute(), GetIntersection(), Loop_MaxQ(), and ToyMC().

MAX_COMB

int CgemLineFit::MAX_COMB

Definition at line 65 of file CgemLineFit.h.

Referenced by CgemLineFit(), and Loop_All().

MC

bool CgemLineFit::MC

Definition at line 64 of file CgemLineFit.h.

Referenced by CgemLineFit(), and ToyMC().

MinQ_Clus2D

double CgemLineFit::MinQ_Clus2D

Definition at line 120 of file CgemLineFit.h.

Referenced by CgemLineFit(), Data_Max(), Loop_All(), and Loop_MaxQ().

Nmax

int CgemLineFit::Nmax

Definition at line 67 of file CgemLineFit.h.

Referenced by CgemLineFit(), GetNMaxQ(), and Loop_MaxQ().

NVComb

int CgemLineFit::NVComb

Definition at line 119 of file CgemLineFit.h.

Referenced by Data_Max(), execute(), Loop_All(), and Loop_MaxQ().

NXComb

int CgemLineFit::NXComb

Definition at line 118 of file CgemLineFit.h.

Referenced by Data_Max(), execute(), Loop_All(), and Loop_MaxQ().

Switch_CCmTPC

int CgemLineFit::Switch_CCmTPC

Definition at line 68 of file CgemLineFit.h.

Referenced by CgemLineFit(), Data_Closest(), Data_Max(), execute(), Loop_All(), and Loop_MaxQ().

TEST_N

int CgemLineFit::TEST_N

Definition at line 66 of file CgemLineFit.h.

Referenced by CgemLineFit(), execute(), Loop_All(), and Loop_MaxQ().

---

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

• CgemLineFit.h
• CgemLineFit.cxx