CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Public Attributes | List of all members
CgemLineFit Class Reference

#include <CgemLineFit.h>

+ Inheritance diagram for CgemLineFit:

Public Member Functions

 CgemLineFit (const std::string &name, ISvcLocator *pSvcLocator)
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
bool Data_Max ()
 
bool Data_Closest ()
 
bool Loop_All ()
 
bool Loop_MaxQ ()
 
bool ToyMC ()
 
void OrderClusters ()
 
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 IncidentAngle (StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[])
 
void GetIntersection ()
 
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 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 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 47 of file CgemLineFit.h.

Constructor & Destructor Documentation

◆ CgemLineFit()

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

Definition at line 111 of file CgemLineFit.cxx.

111 :
112 Algorithm(name, pSvcLocator)
113{
114 declareProperty("sigma", sigma = 0.13);// unit:mm
115 declareProperty("Run10_Flag", Run10_Flag = false);// true:4 clusters, false: 3 clusters
116 declareProperty("Align_Flag", Align_Flag = false);// Alignment
117 declareProperty("Switch_CCmTPC", Switch_CCmTPC = 0);// 0:Center charge 1:mTPC
118 declareProperty("Method", Method = 0); // 0: ToyMC 1: max_charge 2:closest to intersection 3: Loop_all 4: Loop_MaxQ
119 declareProperty("MAX_COMB", MAX_COMB = 50); // The maximum number of combinations for method 3
120 declareProperty("MAX_Clusters", MAX_Clusters = 500); // The maximum number of combinations for method 3
121 declareProperty("MaxClu_Sheet", Nmax = 3); //Max number of clusters on each sheet for method 4
122 declareProperty("Chisq_cut", Chisq_cut = 2000); //Max chisq for the selected set of clusters
123 declareProperty("TEST_N", TEST_N = 0); //set the sheet you want to study the efficiency / resolution. 1 to 4 from top to bottom. 0 use all the layers
124 declareProperty("Debug", debug = 0); //Switch of debug
125 declareProperty("MinQ_Clus2D", MinQ_Clus2D = 0); //Add charge cut on 2D cluster to reduce the noise
126 declareProperty("MC", MC = 0); //remove the hits on 3rd layer for MC
127}
CgemLineFit::Switch_CCmTPC
int Switch_CCmTPC
Definition: CgemLineFit.h:65
CgemLineFit::MinQ_Clus2D
double MinQ_Clus2D
Definition: CgemLineFit.h:110

Member Function Documentation

◆ Data_Closest()

bool CgemLineFit::Data_Closest ( )

Definition at line 814 of file CgemLineFit.cxx.

814 {
815 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
816 int _loop(0);
817 double maxQ_11(0),maxQ_10(0),maxQ_00(0);
818 int max_11(-1),max_10(-1),max_00(-1);
819 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
820 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
821
822 for(;iter!=recCgemClusterCol->end();iter++)
823 {
824 _loop++;
825 RecCgemCluster *cluster = (*iter);
826 int flag=cluster->getflag();
827 int layerid=cluster->getlayerid();
828 double _Q=cluster->getenergydeposit();
829 int sheetid=cluster->getsheetid();
830 double _phi=cluster->getrecphi();
831 if(flag!=2)continue;
832 if(layerid==1&&sheetid==1&&maxQ_11<_Q){
833 maxQ_11=_Q;max_11=_loop; }
834 if(layerid==1&&sheetid==0&&maxQ_10<_Q){
835 maxQ_10=_Q;max_10=_loop; }
836 if(layerid==1&&sheetid==1)cl_11++;
837 if(layerid==0&&sheetid==0&&_phi>0)cl_01++;
838 if(layerid==0&&sheetid==0&&_phi<0)cl_00++;
839 if(layerid==1&&sheetid==0)cl_10++;
840 }
841 clst_11=cl_11;
842 clst_01=cl_01;
843 clst_10=cl_10;
844 clst_00=cl_00;
845 _loop=0;
846 iter =recCgemClusterCol->begin();
847 for(;iter!=recCgemClusterCol->end();iter++)
848 {
849 _loop++;
850 RecCgemCluster *cluster = (*iter);
851 int flag=cluster->getflag();
852 if(flag!=2)continue;
853
854 int layerid=cluster->getlayerid();
855 int sheetid=cluster->getsheetid();
856 double _phi=cluster->getrecphi();
857 double _v=cluster->getrecv();
858 double _Z=cluster->getRecZ();
859 double _Q=cluster->getenergydeposit();
860 double t_phi=cluster->getrecphi_mTPC();
861 double t_v=cluster->getrecv_mTPC();
862 double t_Z=cluster->getRecZ_mTPC();
863
864 if(_loop!=max_11&&_loop!=max_10&&layerid!=0)continue;
865 Vec_flag.push_back(flag);
866 Vec_layerid.push_back(layerid);
867 Vec_sheetid.push_back(sheetid);
868 Vec_layer.push_back(R_layer[layerid]);
869 if(Switch_CCmTPC==1)
870 {
871 Vec_phi.push_back(t_phi);
872 Vec_v.push_back(t_v);
873 Vec_Z.push_back(t_Z);
874 Vec_m_phi.push_back(_phi);
875 Vec_m_v.push_back(_v);
876 Vec_m_Z.push_back(_Z);
877 }
878 else if(Switch_CCmTPC==0)
879 {
880 Vec_phi.push_back(_phi);
881 Vec_v.push_back(_v );
882 Vec_Z.push_back(_Z);
883
884 Vec_m_phi.push_back(t_phi);
885 Vec_m_v.push_back(t_v);
886 Vec_m_Z.push_back(t_Z);
887
888 }
889
890 Vec_Q.push_back(_Q);
891 cluster->setTrkId(0);
892 vecclusters.push_back(cluster);
893
894 }
895 NC = vecclusters.size();
896
897 if(Vec_layer.size()<3)
898 {
899 return false;
900 }
901 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
902 {
903 if(Vec_layerid[i_cl]==2)_clst_2++;
904 if(Vec_layerid[i_cl]==1)_clst_1++;
905 if(Vec_layerid[i_cl]==0)_clst_0++;
906 }
907
908 OrderClusters();
909
910 if(_clst_1==2)
911 l_outer=IniPar(Vec_phi[0],Vec_Z[0],1,Vec_phi[1],Vec_Z[1],1);
912 else return false;
913
914 if(!Layer_cross(0,l_outer))return false;
915 if(_clst_0>1||(Run10_Flag&&_clst_0>2))Filter(0,l_outer);
916
917 return true;
918}
_clst_2
int _clst_2(0)
Vec_Q
vector< double > Vec_Q
Definition: CgemLineFit.cxx:85
Vec_flag
vector< double > Vec_flag
Definition: CgemLineFit.cxx:85
R_layer
double R_layer[3]
Definition: CgemLineFit.cxx:107
Vec_Z
vector< double > Vec_Z
Definition: CgemLineFit.cxx:85
Vec_m_phi
vector< double > Vec_m_phi
Definition: CgemLineFit.cxx:86
Vec_m_Z
vector< double > Vec_m_Z
Definition: CgemLineFit.cxx:86
Vec_sheetid
vector< double > Vec_sheetid
Definition: CgemLineFit.cxx:85
Vec_m_v
vector< double > Vec_m_v
Definition: CgemLineFit.cxx:86
l_outer
StraightLine * l_outer
Definition: CgemLineFit.cxx:93
_clst_0
int _clst_0(0)
NC
int NC
Definition: CgemLineFit.cxx:108
Vec_layer
vector< double > Vec_layer
Definition: CgemLineFit.cxx:85
Vec_v
vector< double > Vec_v
Definition: CgemLineFit.cxx:85
_clst_1
int _clst_1(0)
Vec_phi
vector< double > Vec_phi
Definition: CgemLineFit.cxx:85
Vec_layerid
vector< double > Vec_layerid
Definition: CgemLineFit.cxx:85
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
_Z
#define _Z
Definition: cfortran.h:871
CgemLineFit::OrderClusters
void OrderClusters()
Definition: CgemLineFit.cxx:2061
CgemLineFit::Filter
void Filter(int layerid, StraightLine *l1)
Definition: CgemLineFit.cxx:2189
CgemLineFit::IniPar
StraightLine * IniPar(double phi1, double z1, int i, double phi2, double z2, int j)
Definition: CgemLineFit.cxx:2171
CgemLineFit::Layer_cross
bool Layer_cross(int R_id, StraightLine *l1)
Definition: CgemLineFit.cxx:2888
RecCgemCluster::getrecphi_mTPC
double getrecphi_mTPC(void) const
Definition: RecCgemCluster.h:65
RecCgemCluster::getenergydeposit
double getenergydeposit(void) const
Definition: RecCgemCluster.h:58
RecCgemCluster::getRecZ
double getRecZ(void) const
Definition: RecCgemCluster.h:61
RecCgemCluster::getrecv_mTPC
double getrecv_mTPC(void) const
Definition: RecCgemCluster.h:66
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition: RecCgemCluster.h:55
RecCgemCluster::getflag
int getflag(void) const
Definition: RecCgemCluster.h:57
RecCgemCluster::getRecZ_mTPC
double getRecZ_mTPC(void) const
Definition: RecCgemCluster.h:67
RecCgemCluster::getrecphi
double getrecphi(void) const
Definition: RecCgemCluster.h:59
RecCgemCluster::getrecv
double getrecv(void) const
Definition: RecCgemCluster.h:60
RecCgemCluster::setTrkId
void setTrkId(int trkid)
Definition: RecCgemCluster.h:36
RecCgemCluster::getsheetid
int getsheetid(void) const
Definition: RecCgemCluster.h:56

Referenced by execute().

◆ Data_Max()

bool CgemLineFit::Data_Max ( )

Definition at line 656 of file CgemLineFit.cxx.

656 {
657 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
658 int _loop(0);
659 double maxQ_11(0),maxQ_10(0),maxQ_00(0),maxQ_01(0);
660 int max_11(-1),max_10(-1),max_00(-1),max_01(-1);
661 double maxv_00(0),maxv_01;
662
663 double phi_11(-1),phi_10(-1),phi_00(-1),phi_01(-1);
664 double z_11(0),z_10(0),z_00(0),z_01(0);
665 double cid_11(-1),cid_10(-1),cid_00(-1),cid_01(-1);
666 double Xid_11(-1),Xid_10(-1),Xid_00(-1),Xid_01(-1);
667 double Vid_11(-1),Vid_10(-1),Vid_00(-1),Vid_01(-1);
668
669 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
670 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
671
672 int ic(0);
673 int ic_tot(0);
674 for(;iter!=recCgemClusterCol->end();iter++)
675 {
676 _loop++;
677 RecCgemCluster *cluster = (*iter);
678 int flag=cluster->getflag();
679 int layerid=cluster->getlayerid();
680 double _Q=cluster->getenergydeposit();
681 double _phi=cluster->getrecphi();
682 double _v=cluster->getrecv();
683 int sheetid=cluster->getsheetid();
684
685 int clusterid=cluster->getclusterid();
686 int clusterflagb=cluster->getclusterflagb();
687 int clusterflage=cluster->getclusterflage();
688 double x=R_layer[layerid]*cos(_phi);
689 double y=R_layer[layerid]*sin(_phi);
690 double z=cluster->getRecZ();
691 if(_loop>MAX_Clusters) break;
692
693 if(flag!=2)continue;
694 if(_Q<MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
695
696 tr_id_cluster[ic]=clusterid;
697 tr_x_cluster[ic]=x;
698 tr_y_cluster[ic]=y;
699 tr_z_cluster[ic]=z;
700 tr_Q_cluster[ic]=_Q;
701 tr_phi_cluster[ic]=_phi;
702 tr_layer_cluster[ic]=layerid;
703 tr_sheet_cluster[ic]=sheetid;
704
705 ic_tot++;
706 ic++;
707 if(layerid==1&&sheetid==1&&maxQ_11<_Q){
708 maxQ_11=_Q;max_11=_loop;phi_11=_phi;z_11=z;cid_11=clusterid;Xid_11=clusterflagb;Vid_11=clusterflage;}
709 if(layerid==1&&sheetid==0&&maxQ_10<_Q){
710 maxQ_10=_Q;max_10=_loop;phi_10=_phi;z_10=z;cid_10=clusterid;Xid_10=clusterflagb;Vid_10=clusterflage;}
711 if(layerid==0&&sheetid==0&&_phi>0&&maxQ_01<_Q&&maxv_00!=_v){
712 maxQ_01=_Q;max_01=_loop;phi_01=_phi;z_01=z;cid_01=clusterid;Xid_01=clusterflagb;Vid_01=clusterflage;}
713 if(layerid==0&&sheetid==0&&_phi<0&&maxQ_00<_Q&&maxv_00!=_v){
714 maxQ_00=_Q;max_00=_loop;phi_00=_phi;z_00=z;cid_00=clusterid;Xid_00=clusterflagb;Vid_00=clusterflage;}
715
716 if(layerid==1&&sheetid==1)cl_11++;
717 if(layerid==0&&sheetid==0&&_phi>0)cl_01++;
718 if(layerid==0&&sheetid==0&&_phi<0)cl_00++;
719 if(layerid==1&&sheetid==0)cl_10++;
720 }
721
722
723 tr_ncluster = ic;
724 for(int i=0; i<ic; i++)
725 {
726 tr_Is_selected_cluster[i]=0;
727 if(i==(max_00-1)&&maxQ_00>0) tr_Is_selected_cluster[i] =1;
728 if(i==(max_01-1)&&maxQ_01>0) tr_Is_selected_cluster[i] =1;
729 if(i==(max_10-1)&&maxQ_10>0) tr_Is_selected_cluster[i] =1;
730 if(i==(max_11-1)&&maxQ_11>0) tr_Is_selected_cluster[i] =1;
731 }
732
733
734 clst_11=cl_11;
735 clst_01=cl_01;
736 clst_10=cl_10;
737 clst_00=cl_00;
738 _loop=0;
739
740 iter =recCgemClusterCol->begin();
741 for(;iter!=recCgemClusterCol->end();iter++)
742 {
743 _loop++;
744 RecCgemCluster *cluster = (*iter);
745 int flag=cluster->getflag();
746 if(flag!=2)continue;
747
748 int layerid=cluster->getlayerid();
749 int sheetid=cluster->getsheetid();
750 double _phi=cluster->getrecphi();
751 double _v=cluster->getrecv();
752 double _Z=cluster->getRecZ();
753 double _Q=cluster->getenergydeposit();
754 double t_phi=cluster->getrecphi_mTPC();
755 double t_v=cluster->getrecv_mTPC();
756 double t_Z=cluster->getRecZ_mTPC();
757
758 if(!(_loop==max_00||_loop==max_10||_loop==max_01||_loop==max_11))continue;
759
760 Vec_flag.push_back(flag);
761 Vec_layerid.push_back(layerid);
762 Vec_sheetid.push_back(sheetid);
763 Vec_layer.push_back(R_layer[layerid]);
764
765 if(Switch_CCmTPC==1)
766 {
767 Vec_m_phi.push_back(t_phi);
768 Vec_m_v.push_back(t_v);
769 Vec_m_Z.push_back(t_Z);
770 }
771 else if(Switch_CCmTPC==0)
772 {
773 Vec_phi.push_back(_phi);
774 Vec_v.push_back(_v );
775 Vec_Z.push_back(_Z);
776
777 Vec_m_phi.push_back(t_phi);
778 Vec_m_v.push_back(t_v);
779 Vec_m_Z.push_back(t_Z);
780
781 }
782
783 Vec_Q.push_back(_Q);
784 cluster->setTrkId(0);
785 vecclusters.push_back(cluster);
786
787 }
788 NC = vecclusters.size();
789
790 if(Vec_layer.size()!=4) return false;
791 NXComb = 1;
792 NVComb = 1;
793
794 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
795 {
796 if(Vec_layerid[i_cl]==2)_clst_2++;
797 if(Vec_layerid[i_cl]==1)_clst_1++;
798 if(Vec_layerid[i_cl]==0)_clst_0++;
799 }
800
801 OrderClusters();
802
803 if(_clst_1==2)
804 l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
805 else return false;
806
807 return true;
808}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition: RecCgemCluster.h:53
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition: RecCgemCluster.h:71
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition: RecCgemCluster.h:72

Referenced by execute().

◆ Discard()

void CgemLineFit::Discard ( int  layer)

Definition at line 2163 of file CgemLineFit.cxx.

2164{
2165 for(int i=0;i<Vec_layer.size();i++)
2166 {
2167 if(Vec_layerid[i]==layer) erase(i);
2168 }
2169}
CgemLineFit::erase
void erase(int i)
Definition: CgemLineFit.cxx:2130

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 1951 of file CgemLineFit.cxx.

1952{
1953
1954 HepPoint3D Up,Down;
1955 StraightLine l1(dr,phi0,z0,tglam);
1956 double phiV_up[2],phiV_down[2];
1957 bool findinter = false;
1958 if(!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1959 else findinter = Mp->getPointAligned(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1960
1961 if(!findinter) return 9999999999999;
1962
1963 double dx1=Min_diff(x,phiV_down[0]);
1964 double dx2=Min_diff(x,phiV_up[0]);
1965
1966 if(dx1<dx2)return pow(fabs(phiV_down[1]-V),2.0);
1967 else return pow(fabs(phiV_up[1]-V),2.0);
1968
1969}
Mp
CgemMidDriftPlane * Mp
Definition: CgemLineFit.cxx:88
Align_FLAG
bool Align_FLAG(false)
CgemLineFit::Min_diff
static double Min_diff(double arg1, double arg2)
Definition: CgemLineFit.cxx:2821
CgemMidDriftPlane::getPointIdealGeom
bool getPointIdealGeom(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
Definition: CgemMidDriftPlane.cxx:24
CgemMidDriftPlane::getPointAligned
bool getPointAligned(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
Definition: CgemMidDriftPlane.cxx:164

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 1924 of file CgemLineFit.cxx.

1925{
1926 if(dr>R) return 9999999999999;
1927 StraightLine l1(dr,phi0,dz,tanl);
1928 double phiV_up[2],phiV_down[2];
1929 HepPoint3D Up,Down;
1930 double phiV_up1_old[2],phiV_down1_old[2];
1931 HepPoint3D Up1_old,Down1_old;
1932 bool findinter = false;
1933 if(!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1934 else findinter = Mp->getPointAligned(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1935
1936 if(!findinter) return 9999999999999;
1937
1938 phiV_up[0] = fmod((phiV_up[0]),(2.0*TMath::Pi()));
1939 phiV_down[0] = fmod((phiV_down[0]),(2.0*TMath::Pi()));
1940
1941 if(phiV_up[0]<0) phiV_up[0]+=2.0*TMath::Pi();
1942 if(phiV_down[0]<0) phiV_down[0]+=2.0*TMath::Pi();
1943
1944 double dx1=Min_diff(x,phiV_down[0]);
1945 double dx2=Min_diff(x,phiV_up[0]);
1946 double mdv=min(pow(dx1,2.0),pow(dx2,2.0));
1947 return mdv*R*R;
1948}
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27

Referenced by fcn(), and fcn3().

◆ Erase()

void CgemLineFit::Erase ( int  _layerid)

◆ erase()

void CgemLineFit::erase ( int  i)

Definition at line 2130 of file CgemLineFit.cxx.

2131{
2132 vector<double>::iterator _iter_v;
2133 _iter_v=Vec_v.begin();
2134 Vec_v.erase( _iter_v+i );
2135
2136 vector<double>::iterator _iter_phi;
2137 _iter_phi=Vec_phi.begin();
2138 Vec_phi.erase( _iter_phi+i );
2139
2140 vector<double>::iterator _iter_Z;
2141 _iter_Z=Vec_Z.begin();
2142 Vec_Z.erase( _iter_Z+i );
2143
2144 vector<double>::iterator _iter_layer;
2145 _iter_layer=Vec_layer.begin();
2146 Vec_layer.erase( _iter_layer+i );
2147
2148 vector<double >::iterator _iter_layerid;
2149 _iter_layerid=Vec_layerid.begin();
2150 Vec_layerid.erase( _iter_layerid+i );
2151
2152 vector<double >::iterator _iter_Virlayerid;
2153 _iter_Virlayerid=Vec_Virlayerid.begin();
2154 Vec_Virlayerid.erase( _iter_Virlayerid+i );
2155
2156 vector<double>::iterator _iter_flag;
2157 _iter_flag=Vec_flag.begin();
2158 Vec_flag.erase( _iter_flag+i );
2159
2160}
Vec_Virlayerid
vector< double > Vec_Virlayerid
Definition: CgemLineFit.cxx:85

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

◆ Erase_outer()

bool CgemLineFit::Erase_outer ( )

Definition at line 2105 of file CgemLineFit.cxx.

2106{
2107
2108 vector<double>_x,_v;
2109 vector<int>_iter;
2110 int wrong(-1);
2111 for(int i=0;i<3;i++)
2112 {
2113 _x.push_back(Vec_phi[i]);
2114 _v.push_back(Vec_v[i]);
2115 _iter.push_back(i);
2116 }
2117 if(_x[0]==_x[1]&&_v[0]==_v[2])wrong=0;
2118 else if(_x[0]==_x[1]&&_v[1]==_v[2])wrong=1;
2119 else if(_x[0]==_x[2]&&_v[0]==_v[1])wrong=0;
2120 else if(_x[0]==_x[2]&&_v[2]==_v[1])wrong=2;
2121 else if(_x[1]==_x[2]&&_v[1]==_v[0])wrong=1;
2122 else if(_x[1]==_x[2]&&_v[2]==_v[0])wrong=2;
2123 else {
2124 return false;
2125 }
2126 erase(_iter[wrong]);
2127 return true;
2128}

◆ execute()

StatusCode CgemLineFit::execute ( )

Definition at line 399 of file CgemLineFit.cxx.

399 {
400 MsgStream log(msgSvc(), name());
401 log << MSG::INFO << "in execute()" << endreq;
402
403 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
404 if (!eventHeader) {
405 log << MSG::FATAL << "Could not find Event Header" << endreq;
406 return StatusCode::FAILURE;
407 }
408
409 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
410 if (!recCgemClusterCol){
411 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
412 return StatusCode::FAILURE;
413 }
414
415 m_event=eventHeader->eventNumber();
416 m_run=eventHeader->runNumber();
417 event=m_event;
418 run=m_run;
419 tr_run=run;
420 tr_event=event;
421 if((m_event % 100)==0)cout<<" begin evt : "<<event <<endl;
422 //////////////////Collecting Clusters///////////
423
424 // release the memory
425 Vec_layer.clear();
426 Vec_phi.clear();
427 Vec_Z.clear();
428 Vec_v.clear();
429 Vec_layerid.clear();
430 Vec_Virlayerid.clear();
431 Vec_flag.clear();
432 Vec_Q.clear();
433 Vec_sheetid.clear();
434 Vec_m_layer.clear();
435 Vec_m_phi.clear();
436 Vec_m_Z.clear();
437 Vec_m_v.clear();
438 Vec_m_layerid.clear();
439 Vec_m_flag.clear();
440 Vec_m_Q.clear();
441 Vec_m_sheetid.clear();
442 vecclusters.clear();
443 SmartRefVector<RecCgemCluster>().swap(vecclusters);
444
445 _clst_0=0;_clst_1=0;_clst_2=0;
446
447 NXComb = -1;
448 NVComb = -1;
449
450 TStopwatch gtimer;
451 gtimer.Start();
452
453 if(debug) cout<<"start Straight line fit"<<endl;
454 if(Method==0){
455 if(!ToyMC()){
456 tr_Is_fitted = 0;
457 tr_chisq = -1;
458 m_tuple_track->write();
459 return StatusCode::SUCCESS;}}
460 else if(Method==1){
461 if(!Data_Max()){
462 tr_Is_fitted = 0;
463 tr_chisq = -1;
464 m_tuple_track->write();
465 return StatusCode::SUCCESS;}}
466 else if(Method==2){
467 if(!Data_Closest()){
468 tr_Is_fitted = 0;
469 tr_chisq = -1;
470 m_tuple_track->write();
471 return StatusCode::SUCCESS;}}
472 else if(Method==3){
473 if(!Loop_All()){
474 tr_Is_fitted = 0;
475 tr_chisq = -1;
476 m_tuple_track->write();
477 return StatusCode::SUCCESS;}}
478 else if(Method==4){
479 if(!Loop_MaxQ()){
480 tr_Is_fitted = 0;
481 tr_chisq = -1;
482 m_tuple_track->write();
483 return StatusCode::SUCCESS;}}
484
485 if(Vec_layerid.size()<4&&Run10_Flag)
486 return StatusCode::SUCCESS;
487
488 /////Set initial parameter according to the clusters on the outtermost layer and perform the final fit:
489 /////1, fit phi0, drho; 2, fit tanl z0; 3, fit 4 par. at the same time/////
490
491 TMinuit* fit=Fit(l_outer,fa,0);
492 double trkpar[4]={-999,-999,-999,-999};
493 double trkparErr[4]={-999,-999,-999,-999};
494 Int_t ierflg ,npari,nparx,istat3;
495 Double_t fmin,edm,errdef;
496 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat3);
497 for(int i=0; i<4; i++){
498 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
499 double emat[16];
500 fit->mnemat(emat,4);
501 chisq_3=fit->fAmin;
502 registerTrack(m_trackList_tds,m_hitList_tds);
503 Store_Trk(fit,0);
504 delete fit;
505 double chisq_last = chisq_3;
506
507 gtimer.Stop();
508 vector<double> _trkpar=Trans(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
509
510 par0=_trkpar[0];
511 par1=_trkpar[1];
512 par2=_trkpar[2];
513 par3=_trkpar[3];
514 Err_par0=trkparErr[0];
515 Err_par1=trkparErr[1];
516 Err_par2=trkparErr[2];
517 Err_par3=trkparErr[3];
518
519 status3=istat3;
520
521 clst_0=_clst_0;
522 clst_1=_clst_1;
523 clst_2=_clst_2;
524
525 R0=R_layer[0];
526 R1=R_layer[1];
527 R2=R_layer[2];
528
529 // the initial track parameters based on the clusters on the outtermost layer
530 ini_0=l_outer->dr();
531 ini_1=l_outer->phi0();
532 ini_2=l_outer->dz();
533 ini_3=l_outer->tanl();
534
535 if(debug) cout<<"Get_otherIniPar"<<endl;
536 // Get the initial track parameters based on the clusters on middle and innermost layers
537 if(TEST_N==0) Get_OtherIniPar(_clst_0,_clst_1,_clst_2);
538
539 // 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
540 // _clst_0=2;
541 // _clst_1=2;
542 // _clst_2=0;
543
544 if(debug) cout<<"Rec_Clusters"<<endl;
545 Rec_Clusters();
546 // Store the x,v,z information of clusters by mTPC method on the avialable layers
547 if(Switch_CCmTPC==1) Rec_Clusters_mTPC();
548
549 if(debug) cout<<"Fit_Clusters"<<endl;
550 // Store the intersections between fitted line and Cgem
551 Fit_Clusters(trkpar);
552
553 if(debug) cout<<"GetIntersection"<<endl;
554 // Store the intersections according to the extroplation of the clusters on other layers
555 GetIntersection();
556
557 if(debug) cout<<"GetRes"<<endl;
558 // Get the resolution
559 StraightLine l_fit(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
560 GetRes(&l_fit,TEST_N);
561 m_tuple->write();
562
563 tr_drho =trkpar[0];
564 tr_phi0 =trkpar[1];
565 tr_dz0 =trkpar[2];
566 tr_tanLam=trkpar[3];
567 tr_Is_fitted=1;
568 tr_chisq=chisq_last;
569 m_tuple_track->write();
570 if(debug) cout<<"Finish all"<<endl;
571
572 delete l_outer;
573 return StatusCode::SUCCESS;
574
575}
Vec_m_layer
vector< double > Vec_m_layer
Definition: CgemLineFit.cxx:86
Vec_m_layerid
vector< double > Vec_m_layerid
Definition: CgemLineFit.cxx:86
fa
int fa(63)
Vec_m_Q
vector< double > Vec_m_Q
Definition: CgemLineFit.cxx:86
Vec_m_flag
vector< double > Vec_m_flag
Definition: CgemLineFit.cxx:86
Vec_m_sheetid
vector< double > Vec_m_sheetid
Definition: CgemLineFit.cxx:86
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
CgemLineFit::Rec_Clusters_mTPC
void Rec_Clusters_mTPC()
Definition: CgemLineFit.cxx:2559
CgemLineFit::Loop_All
bool Loop_All()
Definition: CgemLineFit.cxx:925
CgemLineFit::Data_Max
bool Data_Max()
Definition: CgemLineFit.cxx:656
CgemLineFit::Fit_Clusters
void Fit_Clusters(double par[])
Definition: CgemLineFit.cxx:2610
CgemLineFit::Trans
static vector< double > Trans(double par0, double par1, double par2, double par3)
Definition: CgemLineFit.cxx:2904
CgemLineFit::Loop_MaxQ
bool Loop_MaxQ()
Definition: CgemLineFit.cxx:1351
CgemLineFit::Data_Closest
bool Data_Closest()
Definition: CgemLineFit.cxx:814
CgemLineFit::Fit
TMinuit * Fit(StraightLine *l_outer, int sheetflag, int typ)
Definition: CgemLineFit.cxx:1971
CgemLineFit::Get_OtherIniPar
void Get_OtherIniPar(int clst_0, int clst_1, int clst_2)
Definition: CgemLineFit.cxx:2453
CgemLineFit::GetIntersection
void GetIntersection()
Definition: CgemLineFit.cxx:2637
CgemLineFit::Store_Trk
void Store_Trk(TMinuit *fit, int trkid)
Definition: CgemLineFit.cxx:2987
CgemLineFit::Rec_Clusters
void Rec_Clusters()
Definition: CgemLineFit.cxx:2507
CgemLineFit::registerTrack
StatusCode registerTrack(RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
Definition: CgemLineFit.cxx:1806
StraightLine::dz
double dz(void) const
Definition: StraightLine.h:48
StraightLine::dr
double dr(void) const
returns an element of parameters.
Definition: StraightLine.h:45
StraightLine::tanl
double tanl(void) const
Definition: StraightLine.h:49
StraightLine::phi0
double phi0(void) const
Definition: StraightLine.h:46

◆ fcn()

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

Definition at line 1859 of file CgemLineFit.cxx.

1860{
1861 double chisq = 0;
1862 for (int i=0;i<Vec_Virlayerid.size();i++)
1863 {
1864
1865 if(Vec_Virlayerid[i]==5&&bool(sheet_flag&f21))continue;
1866 if(Vec_Virlayerid[i]==4&&bool(sheet_flag&f20))continue;
1867 if(Vec_Virlayerid[i]==3&&bool(sheet_flag&f11))continue;
1868 if(Vec_Virlayerid[i]==2&&bool(sheet_flag&f10))continue;
1869 if(Vec_Virlayerid[i]==1&&bool(sheet_flag&f01))continue;
1870 if(Vec_Virlayerid[i]==0&&bool(sheet_flag&f00))continue;
1871 double dx1=dx(Vec_Virlayerid[i],Vec_layer[i],par[0],par[1],dz_set,tanl_set,Vec_phi[i]);
1872 chisq+=dx1/sigma2;
1873 }
1874 f = chisq;
1875}
f21
int f21(1)
f11
int f11(2)
sigma2
double sigma2(0)
f01
int f01(8)
f20
int f20(32)
f00
int f00(4)
f10
int f10(16)
tanl_set
double tanl_set
Definition: CgemLineFit.cxx:106
dz_set
double dz_set
Definition: CgemLineFit.cxx:106
sheet_flag
int sheet_flag(0)
CgemLineFit::dx
static double dx(double layerid, double R, double dr, double phi0, double z0, double tanl, double x)
Definition: CgemLineFit.cxx:1924

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 1877 of file CgemLineFit.cxx.

1878{
1879 double chisq = 0;
1880
1881 for (int i=0;i<Vec_Virlayerid.size();i++)
1882 {
1883 if(Vec_Virlayerid[i]==5&&bool(sheet_flag&f21))continue;
1884 if(Vec_Virlayerid[i]==4&&bool(sheet_flag&f20))continue;
1885 if(Vec_Virlayerid[i]==3&&bool(sheet_flag&f11))continue;
1886 if(Vec_Virlayerid[i]==2&&bool(sheet_flag&f10))continue;
1887 if(Vec_Virlayerid[i]==1&&bool(sheet_flag&f01))continue;
1888 if(Vec_Virlayerid[i]==0&&bool(sheet_flag&f00))continue;
1889
1890 if(Vec_flag[i]!=2)continue;
1891 double dv1=dV(Vec_Virlayerid[i],Vec_layer[i],dr_set,phi0_set,par[0],par[1],Vec_phi[i],Vec_v[i]);
1892 chisq+=dv1/sigma2;
1893 }
1894 f=chisq;
1895}
phi0_set
double phi0_set
Definition: CgemLineFit.cxx:106
dr_set
double dr_set
Definition: CgemLineFit.cxx:106
CgemLineFit::dV
static double dV(int layer, double R, double dr, double phi0, double z0, double tglam, double x, double V)
Definition: CgemLineFit.cxx:1951

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 1898 of file CgemLineFit.cxx.

1899{
1900
1901 double chisq = 0;
1902 double dv1(0),dx1(0) ;
1903 for (int i=0;i<Vec_Virlayerid.size();i++)
1904 {
1905 if(Vec_flag[i]==2){
1906
1907 if(Vec_Virlayerid[i]==5&&bool(sheet_flag&f21))continue;
1908 if(Vec_Virlayerid[i]==4&&bool(sheet_flag&f20))continue;
1909 if(Vec_Virlayerid[i]==3&&bool(sheet_flag&f11))continue;
1910 if(Vec_Virlayerid[i]==2&&bool(sheet_flag&f10))continue;
1911 if(Vec_Virlayerid[i]==1&&bool(sheet_flag&f01))continue;
1912 if(Vec_Virlayerid[i]==0&&bool(sheet_flag&f00))continue;
1913
1914 dv1=dV(Vec_Virlayerid[i],Vec_layer[i],par[0],par[1],par[2],par[3],Vec_phi[i],Vec_v[i]);
1915 chisq+=(dv1/sigma2);
1916 dx1=dx(Vec_Virlayerid[i],Vec_layer[i],par[0],par[1],par[2],par[3],Vec_phi[i]);
1917 chisq+=dx1/sigma2;
1918 }
1919 }
1920 f=chisq;
1921}

Referenced by Fit().

◆ Filter()

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

Definition at line 2189 of file CgemLineFit.cxx.

2190{
2191 if(!Layer_cross(layerid,l1))
2192 {
2193 Discard(layerid);
2194 }
2195
2196 vector <int>number;
2197 vector <double>dst;
2198 vector <double>dst_up;
2199 vector <double>dst_down;
2200 vector <double>v_x;
2201 vector <double>v_v;
2202
2203 if(Run10_Flag)
2204 {
2205 for (int i =0; i<Vec_layer.size();i++)
2206 {
2207 if(Vec_layerid[i]!=layerid)continue;
2208 number.push_back(i);
2209 // obtain the distance between the up/down side of crosssections to the given position (phi, V)
2210 dst_up.push_back(Min_dis_up(layerid,Vec_phi[i],Vec_v[i],l1));
2211 dst_down.push_back(Min_dis_down(layerid,Vec_phi[i],Vec_v[i],l1));
2212
2213 }
2214
2215 for (int j =0; j<number.size();j++)
2216 {
2217 int i=number[j];
2218 }
2219
2220 // sort the clusters according to the up side distance
2221 int Loop;
2222 Loop=0;
2223 while(Loop!=number.size()-1){
2224 Loop=0;
2225 for (int j=0 ;j<number.size()-1;j++)
2226 {
2227 if(dst_up[j]>dst_up[j+1])
2228 {
2229 swap(dst_up[j],dst_up[j+1]);
2230 swap(dst_down[j],dst_down[j+1]);
2231 swap(number[j],number[j+1]);
2232 }
2233 else Loop++;
2234 }
2235 }
2236
2237 // set the index of the nearest cluster
2238 for (int j =0; j<number.size();j++)
2239 {
2240 int i=number[j];
2241 }
2242
2243 // remove the first member
2244 number.erase(number.begin());
2245 dst_up.erase(dst_up.begin());
2246 dst_down.erase(dst_down.begin());
2247
2248 // what's the difference to previous operation?
2249 for (int j =0; j<number.size();j++)
2250 {
2251 int i=number[j];
2252 }
2253
2254 // sort the clusters according to the down side distance
2255 Loop=0;
2256 while(Loop!=number.size()-1){
2257 Loop=0;
2258 for (int j=0 ;j<number.size()-1;j++)
2259 {
2260 if(dst_down[j]>dst_down[j+1])
2261 {
2262 swap(dst_down[j],dst_down[j+1]);
2263 swap(number[j],number[j+1]);
2264 }
2265 else Loop++;
2266 }
2267 }
2268
2269 for (int j =0; j<number.size();j++)
2270 {
2271 int i=number[j];
2272 }
2273
2274 number.erase(number.begin());
2275
2276 }
2277 else // if run10_flag is false
2278 {
2279 for (int i =0; i<Vec_layer.size();i++)
2280 {
2281 if(Vec_layerid[i]!=layerid)continue;
2282 number.push_back(i);
2283
2284 dst.push_back(Min_dis(layerid,Vec_phi[i],Vec_v[i],l1));
2285
2286 }
2287
2288 int Loop(0);
2289
2290 while(Loop!=number.size()-1){
2291 Loop=0;
2292 for (int j=0 ;j<number.size()-1;j++)
2293 {
2294 if(dst[j]>dst[j+1])
2295 {
2296 swap(dst[j],dst[j+1]);
2297 swap(number[j],number[j+1]);
2298 }
2299 else Loop++;
2300 }
2301 }
2302
2303
2304 number.erase(number.begin());
2305
2306 } // end of if run10_flag
2307
2308 for (int j =0; j<number.size();j++)
2309 {
2310 int i=number[j];
2311 }
2312
2313
2314 sort(number.begin(),number.end(),less<int>());
2315 for(int i=number.size()-1;i!=-1;i--)
2316 {
2317 erase(number[i]);
2318 }
2319
2320 return;
2321}
CgemLineFit::Min_dis_up
double Min_dis_up(int R_id, double x, double z, StraightLine *l1)
Definition: CgemLineFit.cxx:2853
CgemLineFit::Min_dis_down
double Min_dis_down(int R_id, double x, double z, StraightLine *l1)
Definition: CgemLineFit.cxx:2870
CgemLineFit::Min_dis
double Min_dis(int R_id, double x, double z, StraightLine *l1)
Definition: CgemLineFit.cxx:2830
CgemLineFit::Discard
void Discard(int layer)
Definition: CgemLineFit.cxx:2163

Referenced by Data_Closest(), and ToyMC().

◆ FilterClusters()

void CgemLineFit::FilterClusters ( )

Definition at line 2964 of file CgemLineFit.cxx.

2965{
2966 int max_1_0(0),max_0_0(0);
2967 int _size=Vec_flag.size();
2968 max_1_0=_size;max_0_0=_size;
2969
2970 for(int i=0;i<Vec_flag.size();i++)
2971 {
2972 if(1==Vec_layerid[i]&&0==Vec_sheetid[i])
2973 {max_1_0==i;break;}
2974 }
2975
2976 for(int j=0;j<Vec_flag.size();j++)
2977 {
2978 if(0==Vec_layerid[j]&&0==Vec_sheetid[j])
2979 {max_0_0==j;break;}
2980 }
2981 for(int k=Vec_flag.size()-1;k>max_0_0;k--)
2982 {erase(k);}
2983 for(int m=Vec_flag.size()-2;m>max_1_0;m--)
2984 {erase(m);}
2985}

◆ finalize()

StatusCode CgemLineFit::finalize ( )

Definition at line 577 of file CgemLineFit.cxx.

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

◆ Fit()

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

Definition at line 1971 of file CgemLineFit.cxx.

1972{
1973 // typ == 0 : 2D fit on XY plane and SZ plane indipendently, then use the fit results as initial values to perform 3D fit!
1974 // typ == 1 : 2D fit on XY plane
1975 // typ == 2 : 3D fit on XY and SZ plane
1976
1977 double trkpar[4]={-999,-999,-999,-999};
1978 double trkparErr[4]={-999,-999,-999,-999};
1979
1980 sheet_flag=~_sheet_flag;
1981 dz_set=l_outer->dz();
1982 tanl_set=l_outer->tanl();
1983 Int_t ierflg ,npari,nparx,istat1,istat2,istat3;
1984 Double_t fmin,edm,errdef;
1985 Double_t arglist[10];
1986
1987 arglist[0] = 50;
1988 arglist[1] = 0.01;
1989
1990 if(typ==2){
1991 TMinuit *Min3 = new TMinuit(4);
1992 Min3 -> SetFCN(fcn3);
1993 Min3 -> SetPrintLevel(-1);
1994 Min3 -> mnparm(0, "dr" , l_outer->dr(), 0.1 , -1* R_layer[2], R_layer[2], ierflg);
1995 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
1996 Min3 -> mnparm(2, "dz" , l_outer->dz(), 0.1 , -0.5*length, 0.5*length, ierflg);
1997 Min3 -> mnparm(3, "tanL" , l_outer->tanl(), 0.01, -9999, 9999, ierflg);
1998 Min3 -> mnexcm("MIGRAD", arglist, 2, ierflg);
1999
2000 return Min3;
2001 }
2002
2003 if(typ==0){
2004
2005 arglist[0] = 2000000;
2006 arglist[1] = 0.0001;
2007 }
2008 TMinuit *Min = new TMinuit(2);
2009 Min -> SetPrintLevel(-1);
2010 Min -> SetFCN(fcn);
2011 Min -> SetErrorDef(1.0); // For Chisq
2012 Min -> mnparm(0, "dr" ,l_outer->dr(), 0.001, -1*R_layer[2], R_layer[2], ierflg);
2013 Min -> mnparm(1, "phi0" ,l_outer->phi0(), 0.001, -100*acos(-1), 100*acos(-1), ierflg);
2014 arglist[0] = 2000000;
2015 arglist[1] = 0.001;
2016 arglist[0] = 20;
2017 arglist[1] = 0.1;
2018 Min -> mnexcm("MIGRAD", arglist, 2, ierflg);
2019 Min -> mnstat(fmin, edm, errdef, npari, nparx, istat1);
2020 Min -> GetParameter(0, trkpar[0], trkparErr[0]);
2021 Min -> GetParameter(1, trkpar[1], trkparErr[1]);
2022
2023 to_0_2pi(trkpar[1]);
2024 dr_set=trkpar[0];
2025 phi0_set=trkpar[1];
2026 chisq_1=Min->fAmin;
2027 if(typ==1)return Min;
2028 else delete Min;
2029
2030 /////////////////////////Fit to z0 & tanlam////////////
2031 TMinuit *Min2 = new TMinuit(2);
2032 Min2 -> SetFCN(fcn2);
2033 Min2 -> SetPrintLevel(-1);
2034 Min2 -> mnparm(0, "dz" , l_outer->dz(), 0.001 , -0.5*length, 0.5*length, ierflg);
2035 Min2 -> mnparm(1, "tanL" , l_outer->tanl(), 0.001 , -9999 , 9999, ierflg);
2036 Min2 -> mnexcm("MIGRAD", arglist, 2, ierflg);
2037 Min2 -> mnstat(fmin, edm, errdef, npari, nparx, istat2);
2038 for(int i=0; i<2; i++){
2039 Min2 -> GetParameter(i, trkpar[i+2], trkparErr[i+2]);
2040 }
2041 chisq_2=Min2->fAmin;
2042
2043 delete Min2;
2044 //////////////////re-Fit 4 parameters to get correct EDM /////////////
2045
2046
2047 TMinuit *Min3 = new TMinuit(4);
2048 Min3 -> SetFCN(fcn3);
2049 Min3 -> SetPrintLevel(-1);
2050 Min3 -> mnparm(0, "dr" , trkpar[0], 0.001 , -1* R_layer[2], R_layer[2], ierflg);
2051 Min3 -> mnparm(1, "phi0" , trkpar[1], 0.001 , -100*acos(-1), 100*acos(-1), ierflg);
2052 Min3 -> mnparm(2, "dz" , trkpar[2], 0.001 , -0.5*length, 0.5*length, ierflg);
2053 Min3 -> mnparm(3, "tanL" , trkpar[3], 0.001 , -9999, 9999, ierflg);
2054 Min3 -> mnexcm("MIGRAD", arglist, 2, ierflg);
2055
2056 if(typ==0) return Min3;
2057 else delete Min3;
2058}
CgemLineFit::fcn2
static void fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: CgemLineFit.cxx:1877
CgemLineFit::fcn3
static void fcn3(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: CgemLineFit.cxx:1898
CgemLineFit::fcn
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: CgemLineFit.cxx:1859
CgemLineFit::to_0_2pi
static void to_0_2pi(double &arg)
Definition: CgemLineFit.cxx:2814

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

◆ Fit_Clusters()

void CgemLineFit::Fit_Clusters ( double  par[])

Definition at line 2610 of file CgemLineFit.cxx.

2611{
2612 vector<double>Vec_truth;
2613 Vec_truth.push_back(par[0]);
2614 Vec_truth.push_back(par[1]);
2615 Vec_truth.push_back(par[2]);
2616 Vec_truth.push_back(par[3]);
2617
2618 vector<double>clusters=Get_Clusters(Vec_truth);
2619
2620 x_0_up_f= clusters[0];
2621 v_0_up_f= clusters[1];
2622 x_0_down_f= clusters[2];
2623 v_0_down_f= clusters[3];
2624
2625 x_1_up_f= clusters[4];
2626 v_1_up_f= clusters[5];
2627 x_1_down_f= clusters[6];
2628 v_1_down_f= clusters[7];
2629
2630 x_2_up_f= clusters[8];
2631 v_2_up_f= clusters[9];
2632 x_2_down_f= clusters[10];
2633 v_2_down_f= clusters[11];
2634
2635}
CgemLineFit::Get_Clusters
vector< double > Get_Clusters(vector< double >Vec_truth)
Definition: CgemLineFit.cxx:2759

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 2400 of file CgemLineFit.cxx.

2401{
2402 vector<double>_V_par;
2403 _V_par.clear();
2404 double xb=p4.px();
2405 double yb=p4.py();
2406 double zb=p4.pz();
2407 double xa=bp[0];
2408 double ya=bp[1];
2409 double za=bp[2];
2410 double k=-1*(xa*xb+ya*yb)/(xb*xb+yb*yb);
2411 double xc=k*xb+xa;
2412 double yc=k*yb+ya;
2413 double zc=k*zb+za;
2414 HepPoint3D a(xa,ya,za);
2415 HepPoint3D c(xc,yc,zc);
2416 StraightLine l(a,c);
2417 _V_par.push_back(l.dr());
2418 _V_par.push_back(l.phi0());
2419 _V_par.push_back(l.dz());
2420 _V_par.push_back(l.tanl());
2421 return _V_par;
2422}

Referenced by MC_truth().

◆ Get_Clusters()

vector< double > CgemLineFit::Get_Clusters ( vector< double >  Vec_truth)

Definition at line 2759 of file CgemLineFit.cxx.

2760{
2761 vector<double>clusters;
2762 StraightLine l(Vec_truth[0],Vec_truth[1],Vec_truth[2],Vec_truth[3]);
2763 double phiV_up[2],phiV_down[2];
2764 HepPoint3D Up,Down;
2765
2766 double phiV_up_Tmp[2],phiV_down_Tmp[2];
2767 HepPoint3D Up_Tmp,Down_Tmp;
2768
2769 if(!Align_Flag) Mp->getPointIdealGeom(0*2,l,Up,Down,phiV_up,phiV_down);
2770 else
2771 {
2772 Mp->getPointAligned(0,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2773 Mp->getPointAligned(1,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2774 }
2775 to_0_2pi(phiV_up[0]);
2776 to_0_2pi(phiV_down[0]);
2777
2778 clusters.push_back(phiV_up[0]);
2779 clusters.push_back(phiV_up[1]);
2780 clusters.push_back(phiV_down[0]);
2781 clusters.push_back(phiV_down[1]);
2782
2783 if(!Align_Flag) Mp->getPointIdealGeom(1*2,l,Up,Down,phiV_up,phiV_down);
2784 else
2785 {
2786 Mp->getPointAligned(2,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2787 Mp->getPointAligned(3,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2788 }
2789 to_0_2pi(phiV_up[0]);
2790 to_0_2pi(phiV_down[0]);
2791
2792 clusters.push_back(phiV_up[0]);
2793 clusters.push_back(phiV_up[1]);
2794 clusters.push_back(phiV_down[0]);
2795 clusters.push_back(phiV_down[1]);
2796
2797 if(!Align_Flag) Mp->getPointIdealGeom(2*2,l,Up,Down,phiV_up,phiV_down);
2798 else
2799 {
2800 Mp->getPointAligned(4,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2801 Mp->getPointAligned(5,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2802 }
2803 to_0_2pi(phiV_up[0]);
2804 to_0_2pi(phiV_down[0]);
2805
2806 clusters.push_back(phiV_up[0]);
2807 clusters.push_back(phiV_up[1]);
2808 clusters.push_back(phiV_down[0]);
2809 clusters.push_back(phiV_down[1]);
2810
2811 return clusters;
2812}

Referenced by Fit_Clusters(), and Get_MCInf().

◆ Get_MCInf()

bool CgemLineFit::Get_MCInf ( )

Definition at line 2326 of file CgemLineFit.cxx.

2326 {
2327 vector<double>Vec_truth;
2328 Vec_truth.clear();
2329 Vec_truth=MC_truth();
2330
2331 vector<double> _Vec_truth=Trans(Vec_truth[0],Vec_truth[1],Vec_truth[2],Vec_truth[3]);
2332
2333 MC_par0=_Vec_truth[0];
2334 MC_par1=_Vec_truth[1];
2335 MC_par2=_Vec_truth[2];
2336 MC_par3=_Vec_truth[3];
2337 MC_px=Vec_truth[4];
2338 MC_py=Vec_truth[5];
2339 MC_pz=Vec_truth[6];
2340
2341 vector<double>V_MC_clusters=Get_Clusters(Vec_truth);
2342 x_0_up_mc=V_MC_clusters[0];
2343 v_0_up_mc=V_MC_clusters[1];
2344 x_0_down_mc=V_MC_clusters[2];
2345 v_0_down_mc=V_MC_clusters[3];
2346
2347 x_1_up_mc=V_MC_clusters[4];
2348 v_1_up_mc=V_MC_clusters[5];
2349 x_1_down_mc=V_MC_clusters[6];
2350 v_1_down_mc=V_MC_clusters[7];
2351
2352 x_2_up_mc=V_MC_clusters[8];
2353 v_2_up_mc=V_MC_clusters[9];
2354 x_2_down_mc=V_MC_clusters[10];
2355 v_2_down_mc=V_MC_clusters[11];
2356 return true;
2357}
CgemLineFit::MC_truth
vector< double > MC_truth()
Definition: CgemLineFit.cxx:2360

Referenced by ToyMC().

◆ Get_OtherIniPar()

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

Definition at line 2453 of file CgemLineFit.cxx.

2454{
2455
2456 if(clst_2==0)
2457 {
2458
2459 if(clst_0>=2&&clst_1>=2){
2460 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
2461
2462 vector<double> inii_par=Trans(l_ini_i->dr(),l_ini_i->phi0(),l_ini_i->dz(),l_ini_i->tanl());
2463 inii_0=inii_par[0];
2464 inii_1=inii_par[1];
2465 inii_2=inii_par[2];
2466 inii_3=inii_par[3];
2467 delete l_ini_i;
2468 }
2469 if(clst_1>=2){
2470 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
2471 vector<double> inim_par=Trans(l_ini_m->dr(),l_ini_m->phi0(),l_ini_m->dz(),l_ini_m->tanl());
2472
2473 inim_0=inim_par[0];
2474 inim_1=inim_par[1];
2475 inim_2=inim_par[2];
2476 inim_3=inim_par[3];
2477 delete l_ini_m;
2478 }
2479
2480 }
2481 else
2482 {
2483 if(clst_0>=2&&clst_1>=2){
2484 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
2485
2486 vector<double> inii_par=Trans(l_ini_i->dr(),l_ini_i->phi0(),l_ini_i->dz(),l_ini_i->tanl());
2487 inii_0=inii_par[0];
2488 inii_1=inii_par[1];
2489 inii_2=inii_par[2];
2490 inii_3=inii_par[3];
2491 delete l_ini_i;
2492 }
2493 if(clst_1>=2){
2494 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
2495 vector<double> inim_par=Trans(l_ini_m->dr(),l_ini_m->phi0(),l_ini_m->dz(),l_ini_m->tanl());
2496
2497 inim_0=inim_par[0];
2498 inim_1=inim_par[1];
2499 inim_2=inim_par[2];
2500 inim_3=inim_par[3];
2501 delete l_ini_m;
2502 }
2503 }
2504}

Referenced by execute().

◆ GetIntersection()

void CgemLineFit::GetIntersection ( )

Definition at line 2637 of file CgemLineFit.cxx.

2638{
2639
2640 if(debug) cout<<"Start Getintersection"<<endl;
2641 Double_t fmin,edm,errdef;
2642 Int_t ierflg ,npari,nparx,istat;
2643 double trkpar[4]={-999,-999,-999,-999};
2644 double trkparErr[4]={-999,-999,-999,-999};
2645
2646 double phiV_up[2],phiV_down[2];
2647 HepPoint3D Up,Down;
2648 double phiV_up_Tmp[2],phiV_down_Tmp[2];
2649 HepPoint3D Up_Tmp,Down_Tmp;
2650
2651 double Ang[3];
2652 StraightLine* l1a=IniPar(Vec_phi[2],Vec_Z[2],1,Vec_phi[3],Vec_Z[3],0);
2653 TMinuit* fit=Fit(l1a,fa-f11,0);
2654 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat);
2655 inter_1_flag=istat;
2656 for(int i=0; i<4; i++){
2657 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2658 StraightLine* l1b=new StraightLine(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
2659 Store_Trk(fit,1000);
2660 delete fit;
2661
2662 if(Align_FLAG)
2663 {
2664 Mp->getPointAligned(2,l1b,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2665 Mp->getPointAligned(3,l1b,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2666 }
2667 else Mp->getPointIdealGeom(1*2,l1b,Up,Down,phiV_up,phiV_down);
2668
2669 to_0_2pi(phiV_up[0]);
2670 to_0_2pi(phiV_down[0]);
2671 inter_1_x=phiV_up[0];
2672 inter_1_V=phiV_up[1];
2673 inter_1_z=Up[2];
2674
2675
2676 inter_4_x=phiV_down[0];
2677 inter_4_V=phiV_down[1];
2678 inter_4_z=Down[2];
2679
2680
2681 HepPoint3D l1_up(cos(Vec_phi[0])*R_layer[1],sin(Vec_phi[0])*R_layer[1],Vec_Z[0]);
2682 HepPoint3D l1_down(cos(Vec_phi[1])*R_layer[1],sin(Vec_phi[1])*R_layer[1],Vec_Z[1]);
2683
2684 IncidentAngle(l1b,l1_up,pl_10->getStereoAngle(),Ang);
2685 ang_1=Ang[0];
2686 ang_1_x=Ang[1];
2687 ang_1_v=Ang[2];
2688
2689 IncidentAngle(l1b,l1_down,pl_11->getStereoAngle(),Ang);
2690 ang_4=Ang[0];
2691 ang_4_x=Ang[1];
2692 ang_4_v=Ang[2];
2693
2694 delete l1a;
2695 delete l1b;
2696
2697 StraightLine * l2a=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
2698 fit=Fit(l2a,fa-f01,0);
2699 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat);
2700 inter_2_flag=istat;
2701 for(int i=0; i<4; i++){
2702 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2703 Store_Trk(fit,2000);
2704 delete fit;
2705 StraightLine* l2b=new StraightLine(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
2706 if(Align_FLAG)
2707 {
2708 Mp->getPointAligned(0,l2b,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2709 Mp->getPointAligned(1,l2b,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2710 }
2711 else Mp->getPointIdealGeom(0*2,l2b,Up,Down,phiV_up,phiV_down);
2712 to_0_2pi(phiV_up[0]);
2713 inter_2_x=phiV_up[0];
2714 inter_2_V=phiV_up[1];
2715 inter_2_z=Up[2];
2716 HepPoint3D l0_up(cos(Vec_phi[2])*R_layer[0],sin(Vec_phi[2])*R_layer[0],Vec_Z[2]);
2717 IncidentAngle(l2b,l0_up,pl_00->getStereoAngle(),Ang);
2718 ang_2=Ang[0];
2719 ang_2_x=Ang[1];
2720 ang_2_v=Ang[2];
2721 delete l2a;
2722 delete l2b;
2723
2724 if(Run10_Flag)
2725 {
2726
2727 StraightLine * l3a=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
2728 fit=Fit(l3a,fa-f00,0);
2729 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat);
2730 inter_3_flag=istat;
2731 for(int i=0; i<4; i++){
2732 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2733 Store_Trk(fit,3000);
2734 delete fit;
2735 StraightLine* l3b=new StraightLine(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
2736
2737 if(Align_FLAG)
2738 {
2739 Mp->getPointAligned(0,l3b,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2740 Mp->getPointAligned(1,l3b,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2741 }
2742 else Mp->getPointIdealGeom(0*2,l3b,Up,Down,phiV_up,phiV_down);
2743 to_0_2pi(phiV_down[0]);
2744 inter_3_x=phiV_down[0];
2745 inter_3_V=phiV_down[1];
2746 inter_3_z=Down[2];
2747
2748 HepPoint3D l0_down(cos(Vec_phi[1])*R_layer[1],sin(Vec_phi[1])*R_layer[1],Vec_Z[1]);
2749 IncidentAngle(l3b,l0_down,pl_00->getStereoAngle(),Ang);
2750 ang_3=Ang[0];
2751 ang_3_x=Ang[1];
2752 ang_3_v=Ang[2];
2753 delete l3a;
2754 delete l3b;
2755 }
2756
2757}
pl_00
CgemGeoReadoutPlane * pl_00
Definition: CgemLineFit.cxx:95
pl_11
CgemGeoReadoutPlane * pl_11
Definition: CgemLineFit.cxx:97
pl_10
CgemGeoReadoutPlane * pl_10
Definition: CgemLineFit.cxx:96
CgemLineFit::IncidentAngle
void IncidentAngle(StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[])
Definition: CgemLineFit.cxx:3039

Referenced by execute().

◆ GetNMaxQ()

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

Definition at line 2924 of file CgemLineFit.cxx.

2925{
2926 // chose the N X clusters with the largest charge and save corresponding index
2927 vector<int> L_11_s;
2928 vector<double> Q_11_s;
2929 double t, tt;
2930 if (L_11.size()>Nmax)
2931 {
2932 for(int i=0;i<Nmax;i++)
2933 {
2934 for(int j=i+1;j<L_11.size();j++)
2935 {
2936 if(Q_11[i]<Q_11[j])
2937 {
2938 t=Q_11[i];
2939 Q_11[i]=Q_11[j];
2940 Q_11[j]=t;
2941
2942 tt=L_11[i];
2943 L_11[i]=L_11[j];
2944 L_11[j]=tt;
2945 }
2946 }
2947 Q_11_s.push_back(Q_11[i]);
2948 L_11_s.push_back(L_11[i]);
2949 }
2950
2951 }
2952 else
2953 {
2954 L_11_s.assign(L_11.begin(), L_11.end());
2955 Q_11_s.assign(Q_11.begin(), Q_11.end());
2956 }
2957 return L_11_s;
2958}
t
int t()
Definition: t.c:1

Referenced by Loop_MaxQ().

◆ GetVirLay()

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

Definition at line 2917 of file CgemLineFit.cxx.

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

Referenced by 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 2425 of file CgemLineFit.cxx.

2426{
2427 double L(500),W(160),H(600);
2428
2429 vector<double>_V_par;
2430 _V_par.clear();
2431 double xb=p4.px();
2432 double yb=p4.py();
2433 double zb=p4.pz();
2434 double xa=bp[0];
2435 double ya=bp[1];
2436 double za=bp[2];
2437 double k=H/yb*-1;
2438 double xc=k*xb+xa;
2439 double yc=k*yb+ya;
2440 double zc=k*zb+za;
2441 hit_x=xc;
2442 hit_y=yc;
2443 hit_z=zc;
2444
2445 if(fabs(xc)>W/2)return false;
2446 if(fabs(zc)>L/2)return false;
2447
2448 else return true;
2449}
H
IMPLICIT REAL *A H
Definition: myXsection.h:1

Referenced by MC_truth().

◆ InAngle()

void CgemLineFit::InAngle ( StraightLine  sl)

◆ IncidentAngle()

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

Definition at line 3039 of file CgemLineFit.cxx.

3040{
3041
3042 const Hep3Vector Normals(cluster[0],cluster[1],0);
3043
3044 HepPoint3D x1=cosmic_ray->x(1);
3045 HepPoint3D x2=cosmic_ray->x(-1);
3046 const Hep3Vector CosmicRay(x2[0]-x1[0],x2[1]-x1[1],x2[2]-x1[2]);
3047 const Hep3Vector z(0,0,1);
3048 const Hep3Vector phi=z.cross(Normals);
3049 Hep3Vector _phi=phi;
3050
3051 const Hep3Vector V=_phi.rotate(-1*theta,Normals);
3052
3053 const Hep3Vector N_V=Normals.cross(V);
3054 const Hep3Vector N_phi=Normals.cross(phi);
3055
3056 Hep3Vector C_phi=CosmicRay-CosmicRay.project(N_phi);
3057 ang[0]=Normals.angle(C_phi);
3058
3059 Hep3Vector C_V=CosmicRay-CosmicRay.project(N_V);
3060 ang[1]=Normals.angle(C_V);
3061 ang[2]=Normals.angle(CosmicRay);
3062
3063
3064}
StraightLine::x
HepPoint3D x(double s=0.) const
returns position after moving s in downwoards
Definition: StraightLine.cxx:81

Referenced by GetIntersection().

◆ IniPar()

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

Definition at line 2171 of file CgemLineFit.cxx.

2172{
2173 int Gi = int(i/2);
2174 int Gj = int(j/2);
2175
2176 HepPoint3D up(R_layer[Gi]*cos(phi1),R_layer[Gi]*sin(phi1),z1);
2177 HepPoint3D down(R_layer[Gj]*cos(phi2),R_layer[Gj]*sin(phi2),z2);
2178
2179 if(Align_Flag)up=Al->point_invTransform(i ,up);// added by Guoaq-2020/03/13
2180 if(Align_Flag)down=Al->point_invTransform(j ,down);
2181
2182 StraightLine* l=new StraightLine(up,down);
2183 return l;
2184
2185}
Al
CgemGeoAlign * Al
Definition: CgemLineFit.cxx:89
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
phi1
Double_t phi1
Definition: DataBase/tau_mode.c:7
CgemGeoAlign::point_invTransform
HepPoint3D point_invTransform(int layer_vir, HepPoint3D pos)
Definition: CgemGeoAlign.cxx:322

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

◆ initialize()

StatusCode CgemLineFit::initialize ( )

Definition at line 130 of file CgemLineFit.cxx.

130 {
131 MsgStream log(msgSvc(), name());
132 sigma2=sigma*sigma;
133 Align_FLAG=Align_Flag;
134 log << MSG::INFO << "in initialize()" << endreq;
135 StatusCode sc;
136 NTuplePtr cosmic(ntupleSvc(), "LineFit/cosmic");
137 if ( cosmic ) m_tuple = cosmic;
138 else {
139 m_tuple = ntupleSvc()->book ("LineFit/cosmic", CLID_ColumnWiseTuple, "N-Tuple example");
140 if ( m_tuple ) {
141
142 sc = m_tuple->addItem ("run",run);
143 sc = m_tuple->addItem ("event",event);
144 sc = m_tuple->addItem ("R0",R0);
145 sc = m_tuple->addItem ("R1",R1);
146 sc = m_tuple->addItem ("R2",R2);
147 sc = m_tuple->addItem ("clst_0",clst_0);
148 sc = m_tuple->addItem ("clst_1",clst_1);
149 sc = m_tuple->addItem ("clst_2",clst_2);
150 sc = m_tuple->addItem ("clst_01",clst_01);
151 sc = m_tuple->addItem ("clst_00",clst_00);
152 sc = m_tuple->addItem ("clst_11",clst_11);
153 sc = m_tuple->addItem ("clst_10",clst_10);
154 sc = m_tuple->addItem ("status1",status1);
155 sc = m_tuple->addItem ("status2",status2);
156 sc = m_tuple->addItem ("status3",status3);
157 sc = m_tuple->addItem ("ini_0",ini_0);
158 sc = m_tuple->addItem ("ini_1",ini_1);
159 sc = m_tuple->addItem ("ini_2",ini_2);
160 sc = m_tuple->addItem ("ini_3",ini_3);
161 sc = m_tuple->addItem ("inim_0",inim_0);
162 sc = m_tuple->addItem ("inim_1",inim_1);
163 sc = m_tuple->addItem ("inim_2",inim_2);
164 sc = m_tuple->addItem ("inim_3",inim_3);
165 sc = m_tuple->addItem ("inii_0",inii_0);
166 sc = m_tuple->addItem ("inii_1",inii_1);
167 sc = m_tuple->addItem ("inii_2",inii_2);
168 sc = m_tuple->addItem ("inii_3",inii_3);
169 sc = m_tuple->addItem ("par0",par0);
170 sc = m_tuple->addItem ("par1",par1);
171 sc = m_tuple->addItem ("par2",par2);
172 sc = m_tuple->addItem ("par3",par3);
173 sc = m_tuple->addItem ("MC_par0",MC_par0);
174 sc = m_tuple->addItem ("MC_par1",MC_par1);
175 sc = m_tuple->addItem ("MC_par2",MC_par2);
176 sc = m_tuple->addItem ("MC_par3",MC_par3);
177 sc = m_tuple->addItem ("MC_px",MC_px);
178 sc = m_tuple->addItem ("MC_py",MC_py);
179 sc = m_tuple->addItem ("MC_pz",MC_pz);
180 sc = m_tuple->addItem ("Err_par0",Err_par0);
181 sc = m_tuple->addItem ("Err_par1",Err_par1);
182 sc = m_tuple->addItem ("Err_par2",Err_par2);
183 sc = m_tuple->addItem ("Err_par3",Err_par3);
184
185 sc = m_tuple->addItem ("x_2_up",x_2_up);
186 sc = m_tuple->addItem ("v_2_up",v_2_up);
187 sc = m_tuple->addItem ("x_1_up",x_1_up);
188 sc = m_tuple->addItem ("v_1_up",v_1_up);
189 sc = m_tuple->addItem ("x_0_up",x_0_up);
190 sc = m_tuple->addItem ("v_0_up",v_0_up);
191 sc = m_tuple->addItem ("x_2_down",x_2_down);
192 sc = m_tuple->addItem ("v_2_down",v_2_down);
193 sc = m_tuple->addItem ("z_1_up",z_1_up);
194 sc = m_tuple->addItem ("z_2_up",z_2_up);
195 sc = m_tuple->addItem ("z_0_up",z_0_up);
196 sc = m_tuple->addItem ("z_1_down",z_1_down);
197 sc = m_tuple->addItem ("z_2_down",z_2_down);
198 sc = m_tuple->addItem ("z_0_down",z_0_down);
199 sc = m_tuple->addItem ("x_1_down",x_1_down);
200 sc = m_tuple->addItem ("v_1_down",v_1_down);
201 sc = m_tuple->addItem ("x_0_down",x_0_down);
202 sc = m_tuple->addItem ("v_0_down",v_0_down);
203
204
205 sc = m_tuple->addItem ("x_2_up_m",x_2_up_m);
206 sc = m_tuple->addItem ("x_1_up_m",x_1_up_m);
207 sc = m_tuple->addItem ("x_0_up_m",x_0_up_m);
208 sc = m_tuple->addItem ("x_0_down_m",x_0_down_m);
209 sc = m_tuple->addItem ("x_1_down_m",x_1_down_m);
210 sc = m_tuple->addItem ("x_2_down_m",x_2_down_m);
211 sc = m_tuple->addItem ("v_2_up_m",v_2_up_m);
212 sc = m_tuple->addItem ("v_1_up_m",v_1_up_m);
213 sc = m_tuple->addItem ("v_0_up_m",v_0_up_m);
214 sc = m_tuple->addItem ("v_0_down_m",v_0_down_m);
215 sc = m_tuple->addItem ("v_1_down_m",v_1_down_m);
216 sc = m_tuple->addItem ("v_2_down_m",v_2_down_m);
217 sc = m_tuple->addItem ("z_2_up_m",z_2_up_m);
218 sc = m_tuple->addItem ("z_1_up_m",z_1_up_m);
219 sc = m_tuple->addItem ("z_0_up_m",z_0_up_m);
220 sc = m_tuple->addItem ("z_0_down_m",z_0_down_m);
221 sc = m_tuple->addItem ("z_1_down_m",z_1_down_m);
222 sc = m_tuple->addItem ("z_2_down_m",z_2_down_m);
223
224 sc = m_tuple->addItem ("x_2_up_f",x_2_up_f);
225 sc = m_tuple->addItem ("x_1_up_f",x_1_up_f);
226 sc = m_tuple->addItem ("x_0_up_f",x_0_up_f);
227 sc = m_tuple->addItem ("x_0_down_f",x_0_down_f);
228 sc = m_tuple->addItem ("x_1_down_f",x_1_down_f);
229 sc = m_tuple->addItem ("x_2_down_f",x_2_down_f);
230 sc = m_tuple->addItem ("v_2_up_f",v_2_up_f);
231 sc = m_tuple->addItem ("v_1_up_f",v_1_up_f);
232 sc = m_tuple->addItem ("v_0_up_f",v_0_up_f);
233 sc = m_tuple->addItem ("v_0_down_f",v_0_down_f);
234 sc = m_tuple->addItem ("v_1_down_f",v_1_down_f);
235 sc = m_tuple->addItem ("v_2_down_f",v_2_down_f);
236
237 sc = m_tuple->addItem ("x_2_up_mc",x_2_up_mc);
238 sc = m_tuple->addItem ("x_1_up_mc",x_1_up_mc);
239 sc = m_tuple->addItem ("x_0_up_mc",x_0_up_mc);
240 sc = m_tuple->addItem ("x_0_down_mc",x_0_down_mc);
241 sc = m_tuple->addItem ("x_1_down_mc",x_1_down_mc);
242 sc = m_tuple->addItem ("x_2_down_mc",x_2_down_mc);
243 sc = m_tuple->addItem ("v_2_up_mc",v_2_up_mc);
244 sc = m_tuple->addItem ("v_1_up_mc",v_1_up_mc);
245 sc = m_tuple->addItem ("v_0_up_mc",v_0_up_mc);
246 sc = m_tuple->addItem ("v_0_down_mc",v_0_down_mc);
247 sc = m_tuple->addItem ("v_1_down_mc",v_1_down_mc);
248 sc = m_tuple->addItem ("v_2_down_mc",v_2_down_mc);
249
250
251 sc = m_tuple->addItem ("inter_1_x",inter_1_x);
252 sc = m_tuple->addItem ("inter_2_x",inter_2_x);
253 sc = m_tuple->addItem ("inter_3_x",inter_3_x);
254 sc = m_tuple->addItem ("inter_4_x",inter_4_x);
255
256 sc = m_tuple->addItem ("inter_1_V",inter_1_V);
257 sc = m_tuple->addItem ("inter_2_V",inter_2_V);
258 sc = m_tuple->addItem ("inter_3_V",inter_3_V);
259 sc = m_tuple->addItem ("inter_4_V",inter_4_V);
260
261 sc = m_tuple->addItem ("inter_1_flag",inter_1_flag);
262 sc = m_tuple->addItem ("inter_2_flag",inter_2_flag);
263 sc = m_tuple->addItem ("inter_3_flag",inter_3_flag);
264 sc = m_tuple->addItem ("inter_4_flag",inter_4_flag);
265
266 sc = m_tuple->addItem ("inter_1_z",inter_1_z);
267 sc = m_tuple->addItem ("inter_2_z",inter_2_z);
268 sc = m_tuple->addItem ("inter_3_z",inter_3_z);
269 sc = m_tuple->addItem ("inter_4_z",inter_4_z);
270
271
272 sc = m_tuple->addItem ("chisq_1",chisq_1);
273 sc = m_tuple->addItem ("chisq_2",chisq_2);
274 sc = m_tuple->addItem ("chisq_3",chisq_3);
275 sc = m_tuple->addItem ("pos_x",pos_x);
276 sc = m_tuple->addItem ("pos_y",pos_y);
277 sc = m_tuple->addItem ("pos_z",pos_z);
278
279 sc = m_tuple->addItem ("hit_x",hit_x);
280 sc = m_tuple->addItem ("hit_y",hit_y);
281 sc = m_tuple->addItem ("hit_z",hit_z);
282
283 sc = m_tuple->addItem ("ang_1",ang_1);
284 sc = m_tuple->addItem ("ang_1_x",ang_1_x);
285 sc = m_tuple->addItem ("ang_1_v",ang_1_v);
286
287 sc = m_tuple->addItem ("ang_2",ang_2);
288 sc = m_tuple->addItem ("ang_2_x",ang_2_x);
289 sc = m_tuple->addItem ("ang_2_v",ang_2_v);
290
291 sc = m_tuple->addItem ("ang_3",ang_3);
292 sc = m_tuple->addItem ("ang_3_x",ang_3_x);
293 sc = m_tuple->addItem ("ang_3_v",ang_3_v);
294
295 sc = m_tuple->addItem ("ang_4",ang_4);
296 sc = m_tuple->addItem ("ang_4_x",ang_4_x);
297 sc = m_tuple->addItem ("ang_4_v",ang_4_v);
298 sc = m_tuple->addItem ("Res_phi",Res_phi);
299 sc = m_tuple->addItem ("Res_V",Res_V);
300 sc = m_tuple->addItem ("Test_phi",Test_phi);
301 sc = m_tuple->addItem ("Test_V",Test_V);
302 }
303 else {
304 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple) << endmsg;
305 return StatusCode::FAILURE;
306 }
307 }
308
309
310
311 NTuplePtr Track(ntupleSvc(), "LineFit/Track");
312 if ( Track ) m_tuple_track = Track;
313 else {
314 m_tuple_track = ntupleSvc()->book ("LineFit/Track", CLID_ColumnWiseTuple, "N-Tuple example");
315 if ( m_tuple_track ) {
316 sc = m_tuple_track->addItem ("run", tr_run);
317 sc = m_tuple_track->addItem ("event", tr_event);
318 sc = m_tuple_track->addItem ("drho", tr_drho);
319 sc = m_tuple_track->addItem ("phi0", tr_phi0);
320 sc = m_tuple_track->addItem ("dz0", tr_dz0);
321 sc = m_tuple_track->addItem ("chisq", tr_chisq);
322 sc = m_tuple_track->addItem ("Is_fitted", tr_Is_fitted);
323 sc = m_tuple_track->addItem ("tanLam", tr_tanLam);
324 sc = m_tuple_track->addItem ("ncluster", tr_ncluster, 0, 5000);
325 sc = m_tuple_track->addItem ("id_cluster", tr_ncluster, tr_id_cluster);
326 sc = m_tuple_track->addItem ("x_cluster", tr_ncluster, tr_x_cluster);
327 sc = m_tuple_track->addItem ("y_cluster", tr_ncluster, tr_y_cluster);
328 sc = m_tuple_track->addItem ("z_cluster", tr_ncluster, tr_z_cluster);
329 sc = m_tuple_track->addItem ("Q_cluster", tr_ncluster, tr_Q_cluster);
330 sc = m_tuple_track->addItem ("phi_cluster", tr_ncluster, tr_phi_cluster);
331 sc = m_tuple_track->addItem ("layer_cluster", tr_ncluster, tr_layer_cluster);
332 sc = m_tuple_track->addItem ("sheet_cluster", tr_ncluster, tr_sheet_cluster);
333 sc = m_tuple_track->addItem ("Is_selected_cluster", tr_ncluster, tr_Is_selected_cluster);
334
335 }
336 else {
337 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple_track) << endmsg;
338 return StatusCode::FAILURE;
339 }
340 }
341
342 IRawDataProviderSvc* irawDataProviderSvc;
343 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
344 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
345 if ( sc.isFailure() ){
346 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
347 return StatusCode::FAILURE;
348 }
349
350 ICgemGeomSvc* icgemGeomSvc;
351 sc = service ("CgemGeomSvc", icgemGeomSvc);
352 m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
353
354 if ( sc.isFailure() ){
355 log << MSG::FATAL << "Could not load CgemGeomSvc!" << endreq;
356 return StatusCode::FAILURE;
357 }
358
359 R_layer[0]=(m_cgemGeomSvc->getCgemLayer(0)->getMiddleROfGapD());
360 R_layer[1]=(m_cgemGeomSvc->getCgemLayer(1)->getMiddleROfGapD());
361 R_layer[2]=(m_cgemGeomSvc->getCgemLayer(2)->getMiddleROfGapD());
362 length= m_cgemGeomSvc->getLengthOfCgem();
363 Mp=m_cgemGeomSvc->getMidDriftPtr() ;
364 Al=m_cgemGeomSvc->getAlignPtr() ;
365
366 GeoLayer0=m_cgemGeomSvc->getCgemLayer(0);
367 GeoLayer1=m_cgemGeomSvc->getCgemLayer(1);
368 GeoLayer2=m_cgemGeomSvc->getCgemLayer(2);
369
370 pl_00 =m_cgemGeomSvc->getReadoutPlane( 0, 0) ;
371 pl_10=m_cgemGeomSvc->getReadoutPlane( 1, 0) ;
372 pl_11 =m_cgemGeomSvc->getReadoutPlane( 1, 1) ;
373 pl_20=m_cgemGeomSvc->getReadoutPlane( 2, 0) ;
374 pl_21 =m_cgemGeomSvc->getReadoutPlane( 2, 1) ;
375
376
377 cout << "CgemLineFit alignment: is it on? " << Align_Flag << endl;
378 if(Align_Flag) {
379 for(int ilay=0; ilay<6; ilay++) {
380 cout << "LAYER " << ilay+1 << endl;
381 cout << "shift dx " << Al->getDx(ilay) << " dy " << Al->getDy(ilay) << " dz " << Al->getDz(ilay) << endl;
382 cout << "rotation Rx " << Al->getRx(ilay) << " Ry " << Al->getRy(ilay) << " Rz " << Al->getRz(ilay) << endl;
383 cout << "midplane radius " << Mp->getR(int(ilay/2)) << endl;
384 }
385 }
386
387 ICgemCalibFunSvc* icgemCalibSvc;
388 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
389 m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
390 if ( sc.isFailure() ){
391 log << MSG::FATAL << "Could not load CgemCalibFunSvc!" << endreq;
392 return StatusCode::FAILURE;
393 }
394
395 return StatusCode::SUCCESS;
396}
GeoLayer2
CgemGeoLayer * GeoLayer2
Definition: CgemLineFit.cxx:92
GeoLayer0
CgemGeoLayer * GeoLayer0
Definition: CgemLineFit.cxx:90
GeoLayer1
CgemGeoLayer * GeoLayer1
Definition: CgemLineFit.cxx:91
pl_20
CgemGeoReadoutPlane * pl_20
Definition: CgemLineFit.cxx:98
pl_21
CgemGeoReadoutPlane * pl_21
Definition: CgemLineFit.cxx:99
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
Track
#define Track
Definition: TTrackBase.h:30
CgemGeoAlign::getDx
double getDx(int layer_vir)
Definition: CgemGeoAlign.h:69
CgemGeoAlign::getRx
double getRx(int layer_vir)
Definition: CgemGeoAlign.h:72
CgemGeoAlign::getRy
double getRy(int layer_vir)
Definition: CgemGeoAlign.h:73
CgemGeoAlign::getRz
double getRz(int layer_vir)
Definition: CgemGeoAlign.h:74
CgemGeoAlign::getDz
double getDz(int layer_vir)
Definition: CgemGeoAlign.h:71
CgemGeoAlign::getDy
double getDy(int layer_vir)
Definition: CgemGeoAlign.h:70
CgemGeoLayer::getMiddleROfGapD
double getMiddleROfGapD() const
Definition: CgemGeoLayer.h:105
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition: CgemGeomSvc.h:48
CgemGeomSvc::getLengthOfCgem
double getLengthOfCgem() const
Definition: CgemGeomSvc.h:42
CgemGeomSvc::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition: CgemGeomSvc.h:140
CgemGeomSvc::getMidDriftPtr
CgemMidDriftPlane * getMidDriftPtr() const
Definition: CgemGeomSvc.h:72
CgemGeomSvc::getAlignPtr
CgemGeoAlign * getAlignPtr() const
Definition: CgemGeomSvc.h:69
CgemMidDriftPlane::getR
double getR(int layer)
Definition: CgemMidDriftPlane.h:19

◆ Layer_cross()

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

Definition at line 2888 of file CgemLineFit.cxx.

2889{
2890
2891 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2892 double phiV_up[2],phiV_down[2];
2893
2894 HepPoint3D Up,Down;
2895 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2896 else Mp->getPointAligned(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2897 if(Up[0]>9999||Up[1]>9999||Down[0]>9999||Down[1]>9999)
2898 return false;
2899 else return true;
2900
2901}

Referenced by Data_Closest(), and Filter().

◆ Loop_All()

bool CgemLineFit::Loop_All ( )

Definition at line 925 of file CgemLineFit.cxx.

925 {
926 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
927 int _loop(0);
928 double maxQ_11(0),maxQ_10(0),maxQ_00(0);
929 int max_11(-1),max_10(-1),max_00(-1);
930 double C_00(0),C_10(0),C_11(0),C_01(0);
931 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;//
932 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;//
933 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;//
934 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;//
935
936 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
937 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
938 vector<int>L_00,L_01,L_10,L_11;
939 int ic(0);
940 for(;iter!=recCgemClusterCol->end();iter++)
941 {
942 _loop++;
943 RecCgemCluster *cluster = (*iter);
944
945 int flag=cluster->getflag();
946 int layerid=cluster->getlayerid();
947 double _Q=cluster->getenergydeposit();
948 double _phi=cluster->getrecphi();
949 double _v=cluster->getrecv();
950 int sheetid=cluster->getsheetid();
951 int clusterid=cluster->getclusterid();
952 int clusterflagb=cluster->getclusterflagb();
953 int clusterflage=cluster->getclusterflage();
954 double x=R_layer[layerid]*cos(_phi);
955 double y=R_layer[layerid]*sin(_phi);
956 double z=cluster->getRecZ();
957 if(_loop>MAX_Clusters) break;
958 ic++;
959 if(flag!=0)continue;
960
961 if(layerid==1&&sheetid==1&&TEST_N!=1){
962 L_11.push_back(_loop); }
963 if(layerid==1&&sheetid==0&&TEST_N!=4){
964 L_10.push_back(_loop); }
965 if(layerid==0&&sheetid==0&&_phi>0&&TEST_N!=2){
966 L_01.push_back(_loop); }
967 if(layerid==0&&sheetid==0&&_phi<0&&TEST_N!=3){
968 L_00.push_back(_loop); }
969 }
970
971 if(TEST_N==1) L_11.push_back(-1);
972 if(TEST_N==2) L_01.push_back(-1);
973 if(TEST_N==3) L_00.push_back(-1);
974 if(TEST_N==4) L_10.push_back(-1);
975
976 double Min_chi(1E10);
977
978 if((L_00.size()*L_01.size()*L_10.size()*L_11.size())>MAX_COMB)
979 {
980 return false;
981 }
982
983 NXComb = L_00.size()*L_01.size()*L_10.size()*L_11.size();
984
985 for(int i_11=0;i_11<L_11.size();i_11++){
986 for(int i_10=0;i_10<L_10.size();i_10++){
987 for(int i_00=0;i_00<L_00.size();i_00++){
988 for(int i_01=0;i_01<L_01.size();i_01++){
989
990 Vec_layer.clear();
991 Vec_phi.clear();
992 Vec_Z.clear();
993 Vec_v.clear();
994 Vec_layerid.clear();
995 Vec_Virlayerid.clear();
996 Vec_flag.clear();
997 Vec_Q.clear();
998 Vec_sheetid.clear();
999 Vec_m_layer.clear();
1000 Vec_m_phi.clear();
1001 Vec_m_Z.clear();
1002 Vec_m_v.clear();
1003 Vec_m_layerid.clear();
1004 Vec_m_flag.clear();
1005 Vec_m_Q.clear();
1006 Vec_m_sheetid.clear();
1007 _loop=0;
1008 iter =recCgemClusterCol->begin();
1009 for(;iter!=recCgemClusterCol->end();iter++)
1010 {
1011 _loop++;
1012 RecCgemCluster *cluster = (*iter);
1013 int flag=cluster->getflag();
1014 if(!(_loop==L_11[i_11]||_loop==L_10[i_10]||_loop==L_00[i_00]||_loop==L_01[i_01]))continue;
1015 int layerid=cluster->getlayerid();
1016 int sheetid=cluster->getsheetid();
1017 double _phi=cluster->getrecphi();
1018 double _v=cluster->getrecv();
1019 double _Z=cluster->getRecZ();
1020 double _Q=cluster->getenergydeposit();
1021 double t_phi=cluster->getrecphi_mTPC();
1022 double t_v=cluster->getrecv_mTPC();
1023 double t_Z=cluster->getRecZ_mTPC();
1024
1025 int vir_layerid = GetVirLay(layerid, _phi);
1026
1027 Vec_flag.push_back(flag);
1028 Vec_layerid.push_back(layerid);
1029 Vec_Virlayerid.push_back(vir_layerid);
1030 Vec_sheetid.push_back(sheetid);
1031 Vec_layer.push_back(R_layer[layerid]);
1032 if(Switch_CCmTPC==1)
1033 {
1034 Vec_phi.push_back(t_phi);
1035 Vec_v.push_back(t_v);
1036 Vec_Z.push_back(t_Z);
1037 Vec_m_phi.push_back(_phi);
1038 Vec_m_v.push_back(_v);
1039 Vec_m_Z.push_back(_Z);
1040 }
1041 else if(Switch_CCmTPC==0)
1042 {
1043 Vec_phi.push_back(_phi);
1044 Vec_v.push_back(_v );
1045 Vec_Z.push_back(_Z);
1046
1047 Vec_m_phi.push_back(t_phi);
1048 Vec_m_v.push_back(t_v);
1049 Vec_m_Z.push_back(t_Z);
1050
1051 }
1052 }// one combo
1053 if(Vec_layerid.size()!=4&&TEST_N==0)continue;
1054 else if(Vec_layerid.size()!=3&&TEST_N>0)continue;
1055
1056 OrderClusters();
1057
1058 StraightLine* l_outer_tmp;
1059
1060 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);
1061 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);
1062 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1063
1064 TMinuit* _fit=Fit(l_outer_tmp,fa,1);
1065 double _chi=_fit->fAmin;
1066 if(_chi<Min_chi)
1067 {
1068 Min_chi=_chi;
1069 if(TEST_N==0) C_00=Vec_phi[3];
1070 else C_00=-99;
1071 C_01=Vec_phi[2];
1072 C_10=Vec_phi[1];
1073 C_11=Vec_phi[0];
1074
1075 }
1076 delete _fit;
1077 delete l_outer_tmp;
1078 }
1079 }
1080 }
1081 }
1082
1083
1084 iter =recCgemClusterCol->begin();
1085 L_00.clear();L_01.clear();L_10.clear();L_11.clear();
1086 _loop=0;
1087 int ic_tot(0);
1088 for(;iter!=recCgemClusterCol->end();iter++)
1089 {
1090 _loop++;
1091 RecCgemCluster *cluster = (*iter);
1092 int flag=cluster->getflag();
1093 if(flag!=2)continue;
1094 ic_tot++;
1095 int layerid=cluster->getlayerid();
1096 double _Q=cluster->getenergydeposit();
1097 double _phi=cluster->getrecphi();
1098 int sheetid=cluster->getsheetid();
1099 if(_Q<MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
1100 if(C_00==_phi||C_01==_phi||C_11==_phi||C_10==_phi){
1101 if(layerid==1&&sheetid==1){
1102 L_11.push_back(_loop); }
1103 if(layerid==1&&sheetid==0){
1104 L_10.push_back(_loop); }
1105 if(layerid==0&&sheetid==0&&_phi>0){
1106 L_01.push_back(_loop); }
1107 if(layerid==0&&sheetid==0&&_phi<0){
1108 L_00.push_back(_loop); }
1109 }
1110
1111 }
1112 Min_chi=1E10;
1113
1114 if(TEST_N==1)L_11.push_back(-1);
1115 if(TEST_N==2)L_01.push_back(-1);
1116 if(TEST_N==3)L_00.push_back(-1);
1117 if(TEST_N==4)L_10.push_back(-1);
1118
1119 clst_11=L_11.size();
1120 clst_01=L_01.size();
1121 clst_10=L_10.size();
1122 clst_00=L_00.size();
1123 if((L_00.size()*L_01.size()*L_10.size()*L_11.size())>MAX_COMB)
1124 {cout<<"xv size much "<<endl;
1125 return false;
1126 }
1127
1128 NVComb = L_00.size()*L_01.size()*L_10.size()*L_11.size();
1129
1130 for(int i_11=0;i_11<L_11.size();i_11++){
1131 for(int i_10=0;i_10<L_10.size();i_10++){
1132 for(int i_00=0;i_00<L_00.size();i_00++){
1133 for(int i_01=0;i_01<L_01.size();i_01++){
1134
1135 Vec_layer.clear();
1136 Vec_phi.clear();
1137 Vec_Z.clear();
1138 Vec_v.clear();
1139 Vec_layerid.clear();
1140 Vec_Virlayerid.clear();
1141 Vec_flag.clear();
1142 Vec_Q.clear();
1143 Vec_sheetid.clear();
1144 Vec_m_layer.clear();
1145 Vec_m_phi.clear();
1146 Vec_m_Z.clear();
1147 Vec_m_v.clear();
1148 Vec_m_layerid.clear();
1149 Vec_m_flag.clear();
1150 Vec_m_Q.clear();
1151 Vec_m_sheetid.clear();
1152 _loop=0;
1153 iter =recCgemClusterCol->begin();
1154 for(;iter!=recCgemClusterCol->end();iter++)
1155 {
1156 _loop++;
1157 RecCgemCluster *cluster = (*iter);
1158 int flag=cluster->getflag();
1159 if(flag!=2)continue;
1160
1161 int layerid=cluster->getlayerid();
1162 int sheetid=cluster->getsheetid();
1163 double _phi=cluster->getrecphi();
1164 double _v=cluster->getrecv();
1165 double _Z=cluster->getRecZ();
1166 double _Q=cluster->getenergydeposit();
1167 double t_phi=cluster->getrecphi_mTPC();
1168 double t_v=cluster->getrecv_mTPC();
1169 double t_Z=cluster->getRecZ_mTPC();
1170 int vir_layerid = GetVirLay(layerid, _phi);
1171
1172 if(!(_loop==L_11[i_11]||_loop==L_10[i_10]||_loop==L_00[i_00]||_loop==L_01[i_01]))continue;
1173
1174 Vec_flag.push_back(flag);
1175 Vec_layerid.push_back(layerid);
1176 Vec_Virlayerid.push_back(vir_layerid);
1177 Vec_sheetid.push_back(sheetid);
1178 Vec_layer.push_back(R_layer[layerid]);
1179
1180 if(Switch_CCmTPC==1)
1181 {
1182 Vec_phi.push_back(t_phi);
1183 Vec_v.push_back(t_v);
1184 Vec_Z.push_back(t_Z);
1185
1186 Vec_m_phi.push_back(_phi);
1187 Vec_m_v.push_back(_v);
1188 Vec_m_Z.push_back(_Z);
1189 }
1190 else if(Switch_CCmTPC==0)
1191 {
1192 Vec_phi.push_back(_phi);
1193 Vec_v.push_back(_v );
1194 Vec_Z.push_back(_Z);
1195
1196 Vec_m_phi.push_back(t_phi);
1197 Vec_m_v.push_back(t_v);
1198 Vec_m_Z.push_back(t_Z);
1199
1200 }
1201
1202 }// one combo
1203
1204 if(Vec_layerid.size()!=4&&TEST_N==0)continue;
1205 else if(Vec_layerid.size()!=3&&TEST_N>0)continue;
1206
1207 OrderClusters();
1208
1209 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;
1210 StraightLine* l_outer_tmp;
1211
1212 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);
1213 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);
1214 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1215
1216 TMinuit* _fit=Fit(l_outer_tmp,fa,2);
1217 double _chi=_fit->fAmin;
1218 if(_chi<Min_chi)
1219 {
1220 Min_chi=_chi;
1221 C_00=L_00[i_00];
1222 C_01=L_01[i_01];
1223 C_10=L_10[i_10];
1224 C_11=L_11[i_11];
1225 }
1226 delete _fit;
1227 delete l_outer_tmp;
1228 }
1229 }
1230 }
1231 }
1232
1233
1234 Vec_layer.clear();
1235 Vec_phi.clear();
1236 Vec_Z.clear();
1237 Vec_v.clear();
1238 Vec_layerid.clear();
1239 Vec_Virlayerid.clear();
1240 Vec_flag.clear();
1241 Vec_Q.clear();
1242 Vec_sheetid.clear();
1243 Vec_m_layer.clear();
1244 Vec_m_phi.clear();
1245 Vec_m_Z.clear();
1246 Vec_m_v.clear();
1247 Vec_m_layerid.clear();
1248 Vec_m_flag.clear();
1249 Vec_m_Q.clear();
1250 Vec_m_sheetid.clear();
1251 _loop=0;
1252 ic=0;
1253 iter =recCgemClusterCol->begin();
1254 for(;iter!=recCgemClusterCol->end();iter++)
1255 {
1256 _loop++;
1257 RecCgemCluster *cluster = (*iter);
1258 int flag=cluster->getflag();
1259 int layerid=cluster->getlayerid();
1260 int sheetid=cluster->getsheetid();
1261 double _phi=cluster->getrecphi();
1262 double _v=cluster->getrecv();
1263 double _Z=cluster->getRecZ();
1264 double _Q=cluster->getenergydeposit();
1265 double t_phi=cluster->getrecphi_mTPC();
1266 double t_v=cluster->getrecv_mTPC();
1267 double t_Z=cluster->getRecZ_mTPC();
1268 int clusterid=cluster->getclusterid();
1269 int clusterflagb=cluster->getclusterflagb();
1270 int clusterflage=cluster->getclusterflage();
1271 double x=R_layer[layerid]*cos(_phi);
1272 double y=R_layer[layerid]*sin(_phi);
1273 double z=cluster->getRecZ();
1274 int vir_layerid = GetVirLay(layerid, _phi);
1275
1276 if(_loop>MAX_Clusters) break;
1277 if(flag!=2)continue;
1278
1279 tr_id_cluster[ic]=clusterid;
1280 tr_x_cluster[ic]=x;
1281 tr_y_cluster[ic]=y;
1282 tr_z_cluster[ic]=z;
1283 tr_Q_cluster[ic]=_Q;
1284 tr_phi_cluster[ic]=_phi;
1285 tr_layer_cluster[ic]=layerid;
1286 tr_sheet_cluster[ic]=sheetid;
1287 tr_Is_selected_cluster[ic] =0;
1288 ic++;
1289 if(!(_loop==C_00||_loop==C_10||_loop==C_01||_loop==C_11))continue;
1290
1291 tr_Is_selected_cluster[ic-1] =1;
1292
1293 Vec_flag.push_back(flag);
1294 Vec_layerid.push_back(layerid);
1295 Vec_Virlayerid.push_back(vir_layerid);
1296 Vec_sheetid.push_back(sheetid);
1297 Vec_layer.push_back(R_layer[layerid]);
1298 if(Switch_CCmTPC==1)
1299 {
1300 Vec_phi.push_back(t_phi);
1301 Vec_v.push_back(t_v);
1302 Vec_Z.push_back(t_Z);
1303 Vec_m_phi.push_back(_phi);
1304 Vec_m_v.push_back(_v);
1305 Vec_m_Z.push_back(_Z);
1306 }
1307 else if(Switch_CCmTPC==0)
1308 {
1309 Vec_phi.push_back(_phi);
1310 Vec_v.push_back(_v );
1311 Vec_Z.push_back(_Z);
1312
1313 Vec_m_phi.push_back(t_phi);
1314 Vec_m_v.push_back(t_v);
1315 Vec_m_Z.push_back(t_Z);
1316
1317 }
1318
1319 cluster->setTrkId(0);
1320 vecclusters.push_back(cluster);
1321
1322 }// one combo
1323 NC = vecclusters.size();
1324 tr_ncluster = ic;
1325
1326 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
1327 {
1328 if(Vec_layerid[i_cl]==2) _clst_2++;
1329 if(Vec_layerid[i_cl]==1) _clst_1++;
1330 if(Vec_layerid[i_cl]==0) _clst_0++;
1331 }
1332
1333 if(Vec_layerid.size()!=4&&TEST_N==0) return false;
1334 else if(Vec_layerid.size()!=3&&TEST_N>0) return false;
1335
1336 OrderClusters();
1337
1338 if(TEST_N==1||TEST_N==4) l_outer=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1339 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);
1340 else if(TEST_N==0) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1341
1342 if(Min_chi<Chisq_cut)
1343 return true;
1344 else return false;
1345}
CgemLineFit::GetVirLay
int GetVirLay(int geolay, double phi)
Definition: CgemLineFit.cxx:2917

Referenced by execute().

◆ Loop_MaxQ()

bool CgemLineFit::Loop_MaxQ ( )

Definition at line 1351 of file CgemLineFit.cxx.

1351 {
1352 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
1353 int _loop(0);
1354 double maxQ_11(0),maxQ_10(0),maxQ_00(0);
1355 int max_11(-1),max_10(-1),max_00(-1);
1356 double C_00(0),C_10(0),C_11(0),C_01(0);
1357 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;//
1358 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;//
1359 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;//
1360 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;//
1361
1362 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
1363 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
1364 vector<int>L_00,L_01,L_10,L_11;
1365 vector<double>Q_00,Q_01,Q_10,Q_11;
1366 int ic(0);
1367 if(debug) cout<<"start Loop_MaxQ method"<<endl;
1368
1369 for(;iter!=recCgemClusterCol->end();iter++)
1370 {
1371 _loop++;
1372 RecCgemCluster *cluster = (*iter);
1373
1374 int flag=cluster->getflag();
1375 int layerid=cluster->getlayerid();
1376 double _Q=cluster->getenergydeposit();
1377 double _phi=cluster->getrecphi();
1378 double _v=cluster->getrecv();
1379 int sheetid=cluster->getsheetid();
1380 int clusterid=cluster->getclusterid();
1381 int clusterflagb=cluster->getclusterflagb();
1382 int clusterflage=cluster->getclusterflage();
1383 double x=R_layer[layerid]*cos(_phi);
1384 double y=R_layer[layerid]*sin(_phi);
1385 double z=cluster->getRecZ();
1386 if(debug) cout<<"layerid is "<<layerid<<endl;
1387 if(_loop>MAX_Clusters) break;
1388 ic++;
1389 if(flag!=0)continue;
1390
1391 if(layerid==1&&sheetid==1){
1392 L_11.push_back(_loop); Q_11.push_back(_Q); }
1393 if(layerid==1&&sheetid==0){
1394 L_10.push_back(_loop); Q_10.push_back(_Q); }
1395 if(layerid==0&&sheetid==0&&_phi>0){
1396 L_01.push_back(_loop); Q_01.push_back(_Q); }
1397 if(layerid==0&&sheetid==0&&_phi<0){
1398 L_00.push_back(_loop); Q_00.push_back(_Q); }
1399
1400 }
1401
1402 // chose the N X clusters with the largest charge and save corresponding index
1403 vector<int>L_00_s,L_01_s,L_10_s,L_11_s;
1404
1405 L_00_s = GetNMaxQ(Q_00,L_00,Nmax);
1406 L_01_s = GetNMaxQ(Q_01,L_01,Nmax);
1407 L_10_s = GetNMaxQ(Q_10,L_10,Nmax);
1408 L_11_s = GetNMaxQ(Q_11,L_11,Nmax);
1409
1410 if(TEST_N==1){L_11_s.clear();L_11_s.push_back(-1);}
1411 else if(TEST_N==2){L_01_s.clear();L_01_s.push_back(-1);}
1412 else if(TEST_N==3){L_00_s.clear();L_00_s.push_back(-1);}
1413 else if(TEST_N==4){L_10_s.clear();L_10_s.push_back(-1);}
1414
1415 if(debug) cout<<"finish the first loop"<<endl;
1416
1417 double Min_chi(1E10);
1418 NXComb = L_00_s.size()*L_01_s.size()*L_10_s.size()*L_11_s.size();
1419
1420 for(int i_11=0;i_11<L_11_s.size();i_11++){
1421 for(int i_10=0;i_10<L_10_s.size();i_10++){
1422 for(int i_00=0;i_00<L_00_s.size();i_00++){
1423 for(int i_01=0;i_01<L_01_s.size();i_01++){
1424
1425 Vec_layer.clear();
1426 Vec_phi.clear();
1427 Vec_Z.clear();
1428 Vec_v.clear();
1429 Vec_layerid.clear();
1430 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1431 Vec_flag.clear();
1432 Vec_Q.clear();
1433 Vec_sheetid.clear();
1434 Vec_m_layer.clear();
1435 Vec_m_phi.clear();
1436 Vec_m_Z.clear();
1437 Vec_m_v.clear();
1438 Vec_m_layerid.clear();
1439 Vec_m_flag.clear();
1440 Vec_m_Q.clear();
1441 Vec_m_sheetid.clear();
1442 _loop=0;
1443 iter =recCgemClusterCol->begin();
1444 for(;iter!=recCgemClusterCol->end();iter++)
1445 {
1446 _loop++;
1447 RecCgemCluster *cluster = (*iter);
1448 int flag=cluster->getflag();
1449 if(!(_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;
1450 int layerid=cluster->getlayerid();
1451 int sheetid=cluster->getsheetid();
1452 double _phi=cluster->getrecphi();
1453 double _v=cluster->getrecv();
1454 double _Z=cluster->getRecZ();
1455 double _Q=cluster->getenergydeposit();
1456 double t_phi=cluster->getrecphi_mTPC();
1457 double t_v=cluster->getrecv_mTPC();
1458 double t_Z=cluster->getRecZ_mTPC();
1459
1460 int vir_layerid = GetVirLay(layerid, _phi);
1461
1462 Vec_flag.push_back(flag);
1463 Vec_layerid.push_back(layerid);
1464 Vec_Virlayerid.push_back(vir_layerid);
1465 Vec_sheetid.push_back(sheetid);
1466 Vec_layer.push_back(R_layer[layerid]);
1467 if(Switch_CCmTPC==1)
1468 {
1469 Vec_phi.push_back(t_phi);
1470 Vec_v.push_back(t_v);
1471 Vec_Z.push_back(t_Z);
1472
1473 Vec_m_phi.push_back(_phi);
1474 Vec_m_v.push_back(_v);
1475 Vec_m_Z.push_back(_Z);
1476 }
1477 else if(Switch_CCmTPC==0)
1478 {
1479 Vec_phi.push_back(_phi);
1480 Vec_v.push_back(_v );
1481 Vec_Z.push_back(_Z);
1482
1483 Vec_m_phi.push_back(t_phi);
1484 Vec_m_v.push_back(t_v);
1485 Vec_m_Z.push_back(t_Z);
1486
1487 }
1488 }// one combo
1489 if(Vec_layerid.size()!=4&&TEST_N==0)return false;
1490 else if(Vec_layerid.size()!=3&&TEST_N>0)return false;
1491 if(debug) cout<<"before OrderCluster"<<endl;
1492 OrderClusters();
1493 // 0 1 2 3
1494 // order is layer 1 , top, down, layer 0, top down
1495
1496 StraightLine* l_outer_tmp;
1497
1498 if(debug) cout<<"before IniPar"<<endl;
1499 // lay0top virlay lay0down virlay
1500 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);
1501 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);
1502 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1503
1504 if(debug) cout<<"before fit"<<endl;
1505 TMinuit* _fit=Fit(l_outer_tmp,fa,1);
1506 double _chi=_fit->fAmin;
1507 if(debug) cout<<"finish the fit and the chisq is "<<_chi<<endl;
1508 if(debug) cout<<"----------------------------------------"<<endl;
1509 if(_chi<Min_chi)
1510 {
1511 Min_chi=_chi;
1512 if(TEST_N==0) C_00=Vec_phi[3];
1513 else C_00=-99; //Why set C_00 to -99 Guoaq--2020/11/01
1514 C_01=Vec_phi[2]; //It is because there is only 3 clusters selected if TEST_N!=0
1515 C_10=Vec_phi[1];
1516 C_11=Vec_phi[0];
1517
1518 }
1519 delete _fit;
1520 delete l_outer_tmp;
1521 }
1522 }
1523 }
1524 }
1525
1526 if(debug) cout<<"finish the second loop"<<endl;
1527
1528 iter =recCgemClusterCol->begin();
1529 L_00.clear();L_01.clear();L_10.clear();L_11.clear();
1530 L_00_s.clear();L_01_s.clear();L_10_s.clear();L_11_s.clear();
1531 Q_00.clear();Q_01.clear();Q_10.clear();Q_11.clear();
1532 _loop=0;
1533 int ic_tot(0);
1534 for(;iter!=recCgemClusterCol->end();iter++)
1535 {
1536 _loop++;
1537 RecCgemCluster *cluster = (*iter);
1538 int flag=cluster->getflag();
1539 if(flag!=2)continue;
1540 ic_tot++;
1541 int layerid=cluster->getlayerid();
1542 double _Q=cluster->getenergydeposit();
1543 double _phi=cluster->getrecphi();
1544 int sheetid=cluster->getsheetid();
1545 if(_Q<MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
1546 if(C_00==_phi||C_01==_phi||C_11==_phi||C_10==_phi){
1547 if(layerid==1&&sheetid==1){
1548 L_11.push_back(_loop); Q_11.push_back(_Q); }
1549 if(layerid==1&&sheetid==0){
1550 L_10.push_back(_loop); Q_10.push_back(_Q); }
1551 if(layerid==0&&sheetid==0&&_phi>0){
1552 L_01.push_back(_loop); Q_01.push_back(_Q); }
1553 if(layerid==0&&sheetid==0&&_phi<0){
1554 L_00.push_back(_loop); Q_00.push_back(_Q); }
1555 }
1556 }
1557
1558 if(debug) cout<<"finish the third loop"<<endl;
1559
1560 Min_chi=1E10;
1561
1562 clst_11=L_11.size();
1563 clst_01=L_01.size();
1564 clst_10=L_10.size();
1565 clst_00=L_00.size();
1566
1567 L_00_s = GetNMaxQ(Q_00,L_00,Nmax);
1568 L_01_s = GetNMaxQ(Q_01,L_01,Nmax);
1569 L_10_s = GetNMaxQ(Q_10,L_10,Nmax);
1570 L_11_s = GetNMaxQ(Q_11,L_11,Nmax);
1571
1572 if(TEST_N==1) {L_11_s.clear();L_11_s.push_back(-1);}
1573 else if(TEST_N==2){L_01_s.clear();L_01_s.push_back(-1);}
1574 else if(TEST_N==3){L_00_s.clear();L_00_s.push_back(-1);}
1575 else if(TEST_N==4){L_10_s.clear();L_10_s.push_back(-1);}
1576
1577 NVComb = L_00_s.size()*L_01_s.size()*L_10_s.size()*L_11_s.size();
1578
1579 for(int i_11=0;i_11<L_11_s.size();i_11++){
1580 for(int i_10=0;i_10<L_10_s.size();i_10++){
1581 for(int i_00=0;i_00<L_00_s.size();i_00++){
1582 for(int i_01=0;i_01<L_01_s.size();i_01++){
1583
1584 Vec_layer.clear();
1585 Vec_phi.clear();
1586 Vec_Z.clear();
1587 Vec_v.clear();
1588 Vec_layerid.clear();
1589 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1590 Vec_flag.clear();
1591 Vec_Q.clear();
1592 Vec_sheetid.clear();
1593 Vec_m_layer.clear();
1594 Vec_m_phi.clear();
1595 Vec_m_Z.clear();
1596 Vec_m_v.clear();
1597 Vec_m_layerid.clear();
1598 Vec_m_flag.clear();
1599 Vec_m_Q.clear();
1600 Vec_m_sheetid.clear();
1601 _loop=0;
1602 iter =recCgemClusterCol->begin();
1603 for(;iter!=recCgemClusterCol->end();iter++)
1604 {
1605 _loop++;
1606 RecCgemCluster *cluster = (*iter);
1607 int flag=cluster->getflag();
1608 if(flag!=2)continue;
1609 if(!(_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;
1610
1611 int layerid=cluster->getlayerid();
1612 int sheetid=cluster->getsheetid();
1613 double _phi=cluster->getrecphi();
1614 double _v=cluster->getrecv();
1615 double _Z=cluster->getRecZ();
1616 double _Q=cluster->getenergydeposit();
1617 double t_phi=cluster->getrecphi_mTPC();
1618 double t_v=cluster->getrecv_mTPC();
1619 double t_Z=cluster->getRecZ_mTPC();
1620 int vir_layerid = GetVirLay(layerid, _phi);
1621
1622 Vec_flag.push_back(flag);
1623 Vec_layerid.push_back(layerid);
1624 Vec_sheetid.push_back(sheetid);
1625 Vec_layer.push_back(R_layer[layerid]);
1626 Vec_Virlayerid.push_back(vir_layerid);
1627
1628
1629 if(Switch_CCmTPC==1)
1630 {
1631 Vec_phi.push_back(t_phi);
1632 Vec_v.push_back(t_v);
1633 Vec_Z.push_back(t_Z);
1634
1635 Vec_m_phi.push_back(_phi);
1636 Vec_m_v.push_back(_v);
1637 Vec_m_Z.push_back(_Z);
1638 }
1639 else if(Switch_CCmTPC==0)
1640 {
1641 Vec_phi.push_back(_phi);
1642 Vec_v.push_back(_v );
1643 Vec_Z.push_back(_Z);
1644
1645 Vec_m_phi.push_back(t_phi);
1646 Vec_m_v.push_back(t_v);
1647 Vec_m_Z.push_back(t_Z);
1648
1649 }
1650
1651 }// one combo
1652 if(Vec_layerid.size()!=4&&TEST_N==0)return false;
1653 else if(Vec_layerid.size()!=3&&TEST_N>0)return false;
1654 OrderClusters();
1655
1656 // discard the combination if V strip is re-used !
1657 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;
1658 StraightLine* l_outer_tmp;
1659
1660 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);
1661 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);
1662 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1663
1664 TMinuit* _fit=Fit(l_outer_tmp,fa,2);
1665 double _chi=_fit->fAmin;
1666 if(_chi<Min_chi)
1667 {
1668 Min_chi=_chi;
1669 C_00=L_00_s[i_00];
1670 C_01=L_01_s[i_01];
1671 C_10=L_10_s[i_10];
1672 C_11=L_11_s[i_11];
1673
1674 }
1675 delete _fit;
1676 delete l_outer_tmp;
1677 }
1678 }
1679 }
1680 }
1681
1682 Vec_layer.clear();
1683 Vec_phi.clear();
1684 Vec_Z.clear();
1685 Vec_v.clear();
1686 Vec_layerid.clear();
1687 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1688 Vec_flag.clear();
1689 Vec_Q.clear();
1690 Vec_sheetid.clear();
1691 Vec_m_layer.clear();
1692 Vec_m_phi.clear();
1693 Vec_m_Z.clear();
1694 Vec_m_v.clear();
1695 Vec_m_layerid.clear();
1696 Vec_m_flag.clear();
1697 Vec_m_Q.clear();
1698 Vec_m_sheetid.clear();
1699 if(debug) cout<<"finish the forth loop"<<endl;
1700 _loop=0;
1701 ic=0;
1702 iter =recCgemClusterCol->begin();
1703 for(;iter!=recCgemClusterCol->end();iter++)
1704 {
1705 _loop++;
1706 RecCgemCluster *cluster = (*iter);
1707 int flag=cluster->getflag();
1708 int layerid=cluster->getlayerid();
1709 int sheetid=cluster->getsheetid();
1710 double _phi=cluster->getrecphi();
1711 double _v=cluster->getrecv();
1712 double _Z=cluster->getRecZ();
1713 double _Q=cluster->getenergydeposit();
1714 double t_phi=cluster->getrecphi_mTPC();
1715 double t_v=cluster->getrecv_mTPC();
1716 double t_Z=cluster->getRecZ_mTPC();
1717 int clusterid=cluster->getclusterid();
1718 int clusterflagb=cluster->getclusterflagb();
1719 int clusterflage=cluster->getclusterflage();
1720 double x=R_layer[layerid]*cos(_phi);
1721 double y=R_layer[layerid]*sin(_phi);
1722 double z=cluster->getRecZ();
1723 int vir_layerid = GetVirLay(layerid, _phi);
1724 if(_loop>MAX_Clusters) break;
1725 if(flag!=2)continue;
1726
1727 tr_id_cluster[ic]=clusterid;
1728 tr_x_cluster[ic]=x;
1729 tr_y_cluster[ic]=y;
1730 tr_z_cluster[ic]=z;
1731 tr_Q_cluster[ic]=_Q;
1732 tr_phi_cluster[ic]=_phi;
1733 tr_layer_cluster[ic]=layerid;
1734 tr_sheet_cluster[ic]=sheetid;
1735 tr_Is_selected_cluster[ic] =0;
1736
1737 ic++;
1738
1739 if(!(_loop==C_00||_loop==C_10||_loop==C_01||_loop==C_11)) continue;
1740 tr_Is_selected_cluster[ic-1] =1;
1741
1742 Vec_flag.push_back(flag);
1743 Vec_layerid.push_back(layerid);
1744 Vec_Virlayerid.push_back(vir_layerid);
1745 Vec_sheetid.push_back(sheetid);
1746 Vec_layer.push_back(R_layer[layerid]);
1747 if(Switch_CCmTPC==1)
1748 {
1749 Vec_phi.push_back(t_phi);
1750 Vec_v.push_back(t_v);
1751 Vec_Z.push_back(t_Z);
1752
1753 Vec_m_phi.push_back(_phi);
1754 Vec_m_v.push_back(_v);
1755 Vec_m_Z.push_back(_Z);
1756 }
1757 else if(Switch_CCmTPC==0)
1758 {
1759 Vec_phi.push_back(_phi);
1760 Vec_v.push_back(_v );
1761 Vec_Z.push_back(_Z);
1762
1763 Vec_m_phi.push_back(t_phi);
1764 Vec_m_v.push_back(t_v);
1765 Vec_m_Z.push_back(t_Z);
1766
1767 }
1768
1769 cluster->setTrkId(0);
1770 vecclusters.push_back(cluster);
1771
1772 }// one combo
1773 NC = vecclusters.size();
1774 tr_ncluster = ic;
1775 if(debug) cout<<"finish the fifth loop"<<endl;
1776
1777 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
1778 {
1779 if(Vec_layerid[i_cl]==2) _clst_2++;
1780 if(Vec_layerid[i_cl]==1) _clst_1++;
1781 if(Vec_layerid[i_cl]==0) _clst_0++;
1782 }
1783 if(debug) cout<<"Vec_layerid.size is "<<Vec_layerid.size()<<endl;
1784
1785 if(Vec_layerid.size()!=4&&TEST_N==0)return false;
1786 else if(Vec_layerid.size()!=3&&TEST_N>0)return false;
1787 OrderClusters();
1788
1789 if(TEST_N==1||TEST_N==4) l_outer=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1790 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);
1791 else if(TEST_N==0) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1792
1793 if(debug) cout<<"Event "<<event<<", Clu. ID are "<<C_00<<", "<<C_01<<", "<<C_10<<", "<<C_11<<", Loop_MaxQ Min_chi = "<<Min_chi<<endl;
1794
1795 if(Min_chi<Chisq_cut)
1796 return true;
1797 else return false;
1798}
CgemLineFit::GetNMaxQ
static vector< int > GetNMaxQ(vector< double > Q_11, vector< int > L_11, int Nmax)
Definition: CgemLineFit.cxx:2924

Referenced by execute().

◆ MC_truth()

vector< double > CgemLineFit::MC_truth ( )

Definition at line 2360 of file CgemLineFit.cxx.

2361{
2362 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
2363
2364 Event::McParticleCol::iterator iter_mc=mcParticleCol->begin();
2365 HepLorentzVector p4_mu(0,0,0,0);
2366 HepLorentzVector p4_pos(0,0,0,0);
2367 for(;iter_mc!=mcParticleCol->end();iter_mc++)
2368 {
2369 if (!(*iter_mc)->decayFromGenerator()) continue;
2370 HepLorentzVector p4_mc(0,0,0,0);
2371
2372 if (13==(*iter_mc)->particleProperty())
2373 {
2374
2375 p4_mu = (*iter_mc)->initialFourMomentum();
2376 p4_pos = (*iter_mc)->initialPosition();
2377 }
2378 }
2379 vector<double>_mc;
2380 Hep3Vector _BP(p4_pos.x()*10,p4_pos.y()*10,p4_pos.z()*10);
2381 pos_x=_BP[0];
2382 pos_y=_BP[1];
2383 pos_z=_BP[2];
2384 _mc=Get4Par(p4_mu,_BP);
2385
2386 _mc.push_back(p4_mu.px());
2387 _mc.push_back(p4_mu.py());
2388 _mc.push_back(p4_mu.pz());
2389 if(HitPos(p4_mu,_BP))
2390 _mc.push_back(1);
2391 else
2392 _mc.push_back(0);
2393
2394 return _mc;
2395}
CgemLineFit::Get4Par
vector< double > Get4Par(HepLorentzVector p4, Hep3Vector bp)
Definition: CgemLineFit.cxx:2400
CgemLineFit::HitPos
bool HitPos(HepLorentzVector p4, Hep3Vector bp)
Definition: CgemLineFit.cxx:2425

Referenced by Get_MCInf().

◆ Min_diff()

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

Definition at line 2821 of file CgemLineFit.cxx.

2822{
2823 to_0_2pi(arg1);to_0_2pi(arg2);
2824 double diff_raw=fabs(arg1-arg2);
2825
2826 if(diff_raw>acos(-1)) return (2*acos(-1)-diff_raw);
2827 else return diff_raw;
2828}

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 2830 of file CgemLineFit.cxx.

2831{
2832 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2833 double phiV_up[2],phiV_down[2];
2834 HepPoint3D Up,Down;
2835
2836 int vir_R_id = GetVirLay(R_id, x);
2837
2838 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2839 else Mp->getPointAligned(vir_R_id,l1,Up,Down,phiV_up,phiV_down);
2840
2841 double ds1=
2842 pow((phiV_up[1]-v),2)+
2843 pow(R_layer[R_id]*Min_diff(phiV_up[0],x),2);
2844
2845 double ds2=
2846 pow((phiV_down[1]-v),2)+
2847 pow(R_layer[R_id]*Min_diff(phiV_down[0],x),2);
2848
2849 return min(ds1,ds2);
2850}
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35

Referenced by Filter().

◆ Min_dis_down()

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

Definition at line 2870 of file CgemLineFit.cxx.

2871{
2872
2873 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2874 double phiV_up[2],phiV_down[2];
2875 HepPoint3D Up,Down;
2876 int vir_R_id = GetVirLay(R_id, x);
2877 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2878 else Mp->getPointAligned(vir_R_id,l1,Up,Down,phiV_up,phiV_down);
2879
2880 double ds2=
2881 pow((phiV_down[1]-v),2)+
2882 pow(R_layer[R_id]*Min_diff(phiV_down[0],x),2);
2883
2884 return ds2;
2885}

Referenced by Filter().

◆ Min_dis_up()

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

Definition at line 2853 of file CgemLineFit.cxx.

2854{
2855 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2856 double phiV_up[2],phiV_down[2];
2857 HepPoint3D Up,Down;
2858 int vir_R_id = GetVirLay(R_id, x);
2859 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2860 else Mp->getPointAligned(vir_R_id,l1,Up,Down,phiV_up,phiV_down);
2861 double ds1=
2862 pow((phiV_up[1]-z),2)+
2863 pow(R_layer[R_id]*Min_diff(phiV_up[0],x),2);
2864
2865
2866 return ds1;
2867}

Referenced by Filter().

◆ OrderClusters()

void CgemLineFit::OrderClusters ( )

Definition at line 2061 of file CgemLineFit.cxx.

2062{
2063 int noswap(0);
2064 while(noswap!=Vec_flag.size()-1)
2065 {
2066 noswap=0;
2067 for(int i=0;i<Vec_flag.size()-1;i++)
2068 {
2069 if(Vec_flag[i]<Vec_flag[i+1])
2070 { Swap(i,i+1);
2071 }
2072 else if(Vec_flag[i]==Vec_flag[i+1]&&Vec_layerid[i]<Vec_layerid[i+1])
2073 { Swap(i,i+1);
2074 }
2075 else if(Vec_flag[i]==Vec_flag[i+1]&&Vec_layerid[i]==Vec_layerid[i+1]&&(Vec_phi[i]<Vec_phi[i+1]))
2076 {
2077 Swap(i,i+1);}
2078 else {noswap++;
2079 }
2080 }
2081 }
2082}
CgemLineFit::Swap
void Swap(int i, int j)
Definition: CgemLineFit.cxx:2084

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

◆ RealPhi()

double CgemLineFit::RealPhi ( double  SimuPhi)

Definition at line 3138 of file CgemLineFit.cxx.

3139{
3140 double RPhi;
3141 if(SimuPhi <= TMath::Pi()) RPhi = SimuPhi;
3142 else RPhi = SimuPhi-2*TMath::Pi();
3143 return RPhi;
3144
3145}

Referenced by ToyMC().

◆ Rec_Clusters()

void CgemLineFit::Rec_Clusters ( )

Definition at line 2507 of file CgemLineFit.cxx.

2508{
2509
2510 for (int i=0;i<Vec_Virlayerid.size();i++)
2511 {
2512
2513 if(Vec_Virlayerid[i]==5)
2514 {
2515 x_2_up = Vec_phi[i];
2516 v_2_up = Vec_v[i];
2517 z_2_up = Vec_Z[i];
2518 }
2519 if(Vec_Virlayerid[i]==4)
2520 {
2521 x_2_down = Vec_phi[i];
2522 v_2_down = Vec_v[i];
2523 z_2_down = Vec_Z[i];
2524 }
2525 if(Vec_Virlayerid[i]==3)
2526 {
2527 x_1_up = Vec_phi[i];
2528 v_1_up = Vec_v[i];
2529 z_1_up = Vec_Z[i];
2530 }
2531
2532 if(Vec_Virlayerid[i]==2)
2533 {
2534 x_1_down = Vec_phi[i];
2535 v_1_down = Vec_v[i];
2536 z_1_down = Vec_Z[i];
2537 }
2538
2539 if(Vec_Virlayerid[i]==1)
2540 {
2541 x_0_up = Vec_phi[i];
2542 v_0_up = Vec_v[i];
2543 z_0_up = Vec_Z[i];
2544 }
2545
2546 if(Vec_Virlayerid[i]==0)
2547 {
2548 x_0_up = Vec_phi[i];
2549 v_0_up = Vec_v[i];
2550 z_0_up = Vec_Z[i];
2551 }
2552
2553 }
2554
2555}

Referenced by execute().

◆ Rec_Clusters_mTPC()

void CgemLineFit::Rec_Clusters_mTPC ( )

Definition at line 2559 of file CgemLineFit.cxx.

2560{
2561
2562 for (int i=0;i<Vec_Virlayerid.size();i++)
2563 {
2564
2565 if(Vec_Virlayerid[i]==5)
2566 {
2567 x_2_up_m = Vec_m_phi[i];
2568 v_2_up_m = Vec_m_v[i];
2569 z_2_up_m = Vec_m_Z[i];
2570 }
2571 if(Vec_Virlayerid[i]==4)
2572 {
2573 x_2_down_m = Vec_m_phi[i];
2574 v_2_down_m = Vec_m_v[i];
2575 z_2_down_m = Vec_m_Z[i];
2576 }
2577 if(Vec_Virlayerid[i]==3)
2578 {
2579 x_1_up_m = Vec_m_phi[i];
2580 v_1_up_m = Vec_m_v[i];
2581 z_1_up_m = Vec_m_Z[i];
2582 }
2583
2584 if(Vec_Virlayerid[i]==2)
2585 {
2586 x_1_down_m = Vec_m_phi[i];
2587 v_1_down_m = Vec_m_v[i];
2588 z_1_down_m = Vec_m_Z[i];
2589 }
2590
2591 if(Vec_Virlayerid[i]==1)
2592 {
2593 x_0_up_m = Vec_m_phi[i];
2594 v_0_up_m = Vec_m_v[i];
2595 z_0_up_m = Vec_m_Z[i];
2596 }
2597
2598 if(Vec_Virlayerid[i]==0)
2599 {
2600 x_0_up_m = Vec_m_phi[i];
2601 v_0_up_m = Vec_m_v[i];
2602 z_0_up_m = Vec_m_Z[i];
2603 }
2604
2605 }
2606
2607}

Referenced by execute().

◆ registerTrack()

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

Definition at line 1806 of file CgemLineFit.cxx.

1807{
1808 MsgStream log(msgSvc(), name());
1809 StatusCode sc;
1810 IDataProviderSvc* eventSvc = NULL;
1811 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1812 if (eventSvc) {
1813 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1814 } else {
1815 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1816 return StatusCode::FAILURE;
1817 //return StatusCode::SUCCESS;
1818 }
1819 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1820
1821
1822 // --- register RecMdcTrack
1823 recMdcTrackCol = new RecMdcTrackCol;
1824 DataObject *aRecMdcTrackCol;
1825 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1826 if(aRecMdcTrackCol != NULL) {
1827 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1828 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1829 }
1830
1831 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
1832 if(sc.isFailure()) {
1833 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1834 return StatusCode::FAILURE;
1835 }
1836 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1837
1838
1839 recMdcHitCol = new RecMdcHitCol;
1840 DataObject *aRecMdcHitCol;
1841 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1842 if(aRecMdcHitCol != NULL) {
1843 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1844 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1845 }
1846
1847 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
1848 if(sc.isFailure()) {
1849 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1850 return StatusCode::FAILURE;
1851 }
1852 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1853
1854 return StatusCode::SUCCESS;
1855}//end of stroeTracks
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: RecMdcHit.h:99
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: RecMdcTrack.h:102

Referenced by execute().

◆ Store_Trk()

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

Definition at line 2987 of file CgemLineFit.cxx.

2988{
2989 double trkpar[4]={-999,-999,-999,-999};
2990 double trkparErr[4]={-999,-999,-999,-999};
2991
2992 Int_t ierflg ,npari,nparx,istat3;
2993 Double_t fmin,edm,errdef;
2994
2995 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat3);
2996 for(int i=0; i<4; i++){
2997 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2998 double emat[16];
2999 fit->mnemat(emat,4);
3000
3001 RecMdcTrack* recMdcTrack = new RecMdcTrack();
3002
3003 double helixPar[5];
3004 helixPar[0]=trkpar[0];
3005 helixPar[1]=trkpar[1];
3006 helixPar[2]=0.0;
3007 helixPar[3]=trkpar[2];
3008 helixPar[4]=trkpar[3];
3009 double errorMat[15];
3010 int k = 0;
3011 errorMat[0]=emat[0];
3012 errorMat[1]=emat[1];
3013 errorMat[2]=0;
3014 errorMat[3]=emat[2];
3015 errorMat[4]=emat[3];
3016 errorMat[5]=emat[5];
3017 errorMat[6]=0;
3018 errorMat[7]=emat[6];
3019 errorMat[8]=emat[7];
3020 errorMat[9]=0;
3021 errorMat[10]=0;
3022 errorMat[11]=0;
3023 errorMat[12]=emat[10];
3024 errorMat[13]=emat[11];
3025 errorMat[14]=emat[15];
3026
3027 recMdcTrack->setChi2(fit->fAmin);
3028 recMdcTrack->setError(errorMat);
3029 recMdcTrack->setHelix(helixPar);
3030 recMdcTrack->setTrackId(trkid);
3031 if(trkid==0){
3032 recMdcTrack->setNcluster(NC);
3033 recMdcTrack->setVecClusters(vecclusters);
3034 m_trackList_tds->push_back(recMdcTrack);
3035 }
3036
3037}
DstMdcTrack::setTrackId
void setTrackId(const int trackId)
Definition: DstMdcTrack.h:84
DstMdcTrack::setError
void setError(double err[15])
Definition: DstMdcTrack.cxx:118
DstMdcTrack::setHelix
void setHelix(double helix[5])
Definition: DstMdcTrack.cxx:104
DstMdcTrack::setChi2
void setChi2(const double chi)
Definition: DstMdcTrack.h:98
RecMdcTrack::setVecClusters
void setVecClusters(ClusterRefVec vecclusters)
Definition: RecMdcTrack.cxx:86
RecMdcTrack::setNcluster
void setNcluster(int ncluster)
Definition: RecMdcTrack.h:83

Referenced by execute(), and GetIntersection().

◆ Swap()

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

Definition at line 2084 of file CgemLineFit.cxx.

2085{
2086
2087 swap(Vec_layer[i],Vec_layer[j]);
2088 swap(Vec_sheetid[i],Vec_sheetid[j]);
2089 swap(Vec_phi[i],Vec_phi[j]);
2090 swap(Vec_Z[i],Vec_Z[j]);
2091 swap(Vec_v[i],Vec_v[j]);
2092 swap(Vec_layerid[i],Vec_layerid[j]);
2093 swap(Vec_Virlayerid[i],Vec_Virlayerid[j]);
2094 swap(Vec_flag[i],Vec_flag[j]);
2095 if(Vec_m_phi.size()==Vec_flag.size()&&(Method==1||Method==2)){
2096 swap(Vec_m_phi[i],Vec_m_phi[j]);
2097 swap(Vec_m_Z[i],Vec_m_Z[j]);
2098 swap(Vec_m_v[i],Vec_m_v[j]);
2099 }
2100}

Referenced by OrderClusters().

◆ to_0_2pi()

void CgemLineFit::to_0_2pi ( double &  arg)
static

Definition at line 2814 of file CgemLineFit.cxx.

2815{
2816
2817 arg = fmod((arg),(2.0*TMath::Pi()));
2818 if(arg<0) arg += 2.0*TMath::Pi();
2819}
arg
double arg(const EvtComplex &c)
Definition: EvtComplex.hh:227

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

◆ ToyMC()

bool CgemLineFit::ToyMC ( )

Definition at line 587 of file CgemLineFit.cxx.

587 {
588
589 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
590 int _loop(0);
591 if(!Get_MCInf()) return false;
592 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
593 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
594 for(;iter!=recCgemClusterCol->end();iter++)
595 {
596 _loop++;
597 RecCgemCluster *cluster = (*iter);
598 int flag=cluster->getflag();
599 if(flag!=2) continue;
600
601 int layerid=cluster->getlayerid();
602 int sheetid=cluster->getsheetid();
603 double _phi=cluster->getrecphi();
604 if(MC) _phi = RealPhi(_phi);
605 double _v=cluster->getrecv();
606 double _Z=cluster->getRecZ();
607 double _Q=cluster->getenergydeposit();
608 int vir_layerid = GetVirLay(layerid, _phi);
609 if(layerid==2||((!Run10_Flag)&&layerid==0&&_phi>acos(-1)))continue;
610
611 Vec_flag.push_back(flag);
612 Vec_layerid.push_back(layerid);
613 Vec_Virlayerid.push_back(vir_layerid);
614 Vec_sheetid.push_back(sheetid);
615 Vec_layer.push_back(R_layer[layerid]);
616
617 Vec_phi.push_back(_phi);
618 Vec_v.push_back(_v );
619 Vec_Z.push_back(_Z);
620
621 Vec_Q.push_back(_Q);
622 cluster->setTrkId(0);
623 vecclusters.push_back(cluster);
624
625 }
626 int NC = vecclusters.size();
627 if(debug) cout<<"NC is "<<NC<<endl;
628
629 if(Vec_layer.size()<3)
630 {
631 return false;
632 }
633
634 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
635 {
636 if(Vec_layerid[i_cl]==2) _clst_2++;
637 if(Vec_layerid[i_cl]==1) _clst_1++;
638 if(Vec_layerid[i_cl]==0) _clst_0++;
639 }
640
641 if(debug) cout<<"clst1, clst0 are "<<_clst_1<<", "<<_clst_0<<endl;
642 OrderClusters();
643 if(_clst_1==2)
644 l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2); // should be this
645 else return false;
646
647 if(_clst_0>2)Filter(0,l_outer);
648
649 return true;
650}
CgemLineFit::Get_MCInf
bool Get_MCInf()
Definition: CgemLineFit.cxx:2326
CgemLineFit::RealPhi
double RealPhi(double SimuPhi)
Definition: CgemLineFit.cxx:3138

Referenced by execute().

◆ Trans()

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

Definition at line 2904 of file CgemLineFit.cxx.

2905{
2906 vector<double> rotat;
2907 rotat.clear();
2908
2909 rotat.push_back(par0);
2910 rotat.push_back(par1);
2911 rotat.push_back(par2);
2912 rotat.push_back(par3);
2913 return rotat;
2914
2915}

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

Member Data Documentation

◆ debug

bool CgemLineFit::debug

Definition at line 60 of file CgemLineFit.h.

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

◆ MAX_COMB

int CgemLineFit::MAX_COMB

Definition at line 62 of file CgemLineFit.h.

Referenced by CgemLineFit(), and Loop_All().

◆ MC

bool CgemLineFit::MC

Definition at line 61 of file CgemLineFit.h.

Referenced by CgemLineFit(), and ToyMC().

◆ MinQ_Clus2D

double CgemLineFit::MinQ_Clus2D

Definition at line 110 of file CgemLineFit.h.

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

◆ Nmax

int CgemLineFit::Nmax

Definition at line 64 of file CgemLineFit.h.

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

◆ NVComb

int CgemLineFit::NVComb

Definition at line 109 of file CgemLineFit.h.

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

◆ NXComb

int CgemLineFit::NXComb

Definition at line 108 of file CgemLineFit.h.

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

◆ Switch_CCmTPC

int CgemLineFit::Switch_CCmTPC

Definition at line 65 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 63 of file CgemLineFit.h.

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

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