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

#include <CgemMilleAlign.h>

+ Inheritance diagram for CgemMilleAlign:

Public Member Functions

 CgemMilleAlign ()
 
 ~CgemMilleAlign ()
 
void initialize (TObjArray *hlist, ICgemGeomSvc *cgemGeomSvc, ICgemCalibFunSvc *cgemFunSvc)
 
void setParam (CgemAliParams &param)
 
bool fillHist (CgemAliEvent *event, CgemAlignPar *alignPar)
 
bool updateAlignPar (CgemAlignPar *alignPar)
 
void clear ()
 
- Public Member Functions inherited from CgemAlignBase
 CgemAlignBase ()
 
virtual ~CgemAlignBase ()
 
virtual void initialize (TObjArray *hlist, ICgemGeomSvc *cgemGeomSvc, ICgemCalibFunSvc *cgemFunSvc)=0
 
virtual void setParam (CgemAliParams &param)=0
 
virtual bool fillHist (CgemAliEvent *event, CgemAlignPar *alignPar)=0
 
virtual bool updateAlignPar (CgemAlignPar *alignPar)=0
 
virtual void clear ()=0
 

Detailed Description

Definition at line 29 of file CgemMilleAlign.h.

Constructor & Destructor Documentation

◆ CgemMilleAlign()

CgemMilleAlign::CgemMilleAlign ( )

Definition at line 33 of file CgemMilleAlign.cxx.

33 {
34 cout << "CgemMilleAlign::CgemMilleAlign()" << endl;
35 for(int lay=0; lay<LAYERNMAX; lay++){
36 m_resiCut_phi[lay] = 7.0; //how to optimize this number guoaq-19-03-21
37 m_resiCut_v[lay] = 7.0; //how to optimize this number guoaq-19-03-21
38 m_docaCut[lay] = 0.5;
39 m_docaCut[lay] = 5.5;
40 }
41}
CgemAlignment::LAYERNMAX
const int LAYERNMAX
Definition: CgemAlignment.h:13

◆ ~CgemMilleAlign()

CgemMilleAlign::~CgemMilleAlign ( )

Definition at line 43 of file CgemMilleAlign.cxx.

43 {
44 cout << "CgemMilleAlign::~CgemMilleAlign()" << endl;
45}

Member Function Documentation

◆ clear()

void CgemMilleAlign::clear ( )
virtual

Implements CgemAlignBase.

Definition at line 47 of file CgemMilleAlign.cxx.

47 {
48 delete m_hresAll_phi;
49 delete m_hresAll_v;
50 delete m_hresAll_z;
51 for(int lay=0; lay<LAYERNMAX; lay++)
52 {
53 delete m_hresLay_phi[lay];
54 delete m_hresLay_v[lay];
55 delete m_hresLay_z[lay];
56 };
57 delete m_hddoca;
58 delete m_pMilleAlign;
59}

◆ fillHist()

bool CgemMilleAlign::fillHist ( CgemAliEvent *  event,
CgemAlignPar *  alignPar 
)
virtual

Implements CgemAlignBase.

Definition at line 515 of file CgemMilleAlign.cxx.

515 {
516
517 IMessageSvc* msgSvc;
518 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
519 MsgStream log(msgSvc, "MilleAlign");
520 log << MSG::INFO << "MilleAlign::fillTree()" << endreq;
521
522 int recFlag;
523 int itrk;
524 int ihit;
525 int lay;
526 int Sheet_Meas;
527 int lay_virtual;
528 int Sheet_virtual;
529 double Phi_Meas;
530 double V_Meas;
531 double Z_Meas;
532 double resi_phi;
533 double resi_v;
534 double resi_phi_O;
535 double resi_v_O;
536 double resi_z;
537 double exp_phi;
538 double exp_v;
539 double exp_z;
540 double resi_comb;
541 double deri_Phi;
542 double deri_V;
543 int direction;
544 double hitSigm_Phi=0.0009;
545 double hitSigm_V=0.2239;
546 double lay_pos[6]; // 4 mis-alignment parameters for each layer
547
548 double trkpar[NTRKPAR];
549 double trkparms[NTRKPARALL]; // track parameters and errors
550
551 // numerical derivative
552 int ipar;
553 int iparGB;
554
555 CgemAliRecTrk* rectrk;
556 CgemAliRecHit* rechit;
557
558 HepPoint3D posUp, posDown;
559 double phiVUp[3], phiVDown[3];
560
561 HepPoint3D posUp_O, posDown_O;
562 double phiVUp_O[3], phiVDown_O[3];
563
564 int ntrk = event -> getNTrk();
565
566 log << MSG::INFO << "Number of tracks : "<<ntrk << endreq;
567
568 for(itrk=0; itrk<ntrk; itrk++){
569 rectrk = event->getRecTrk(itrk);
570 recFlag = rectrk->getStat();
571
572 trkpar[0] = rectrk -> getDr();
573 trkpar[1] = rectrk -> getPhi0();
574 trkpar[2] = rectrk -> getKappa();
575 trkpar[3] = rectrk -> getDz();
576 trkpar[4] = rectrk -> getTanLamda();
577 double chisq = rectrk -> getChisq();
578
579 int nHit = 0;
580 nHit = rectrk -> getNcluster(); //Why this number is two times larger than the inpur value ??? Guoaq -19June08
581
582 HepVector rechelix = rectrk->getHelix();
583 HepVector helix = rectrk->getHelix();
584 HepSymMatrix helixErr = rectrk->getHelixErr();
585
586 if(re_3lay&&nHit<12) continue; // only use the track can pass through all 3 layers of Cgem
587
588 int nHitUse = 0;
589 for(ihit=0; ihit<nHit/2; ihit++){
590 rechit = rectrk -> getRecHit(ihit);
591
592 lay = rechit -> getlayerid();
593 Phi_Meas = rechit -> getrecphi();
594 V_Meas = rechit -> getrecv();
595 Z_Meas = rechit -> getRecZ();
596 Sheet_Meas = rechit -> getsheetid();
597
598 if(lay==0) Sheet_virtual = (Phi_Meas<0)? 0:1;
599 else Sheet_virtual = Sheet_Meas;
600 lay_virtual = lay*2+Sheet_virtual;
601
602 lay_pos[0] = alignPar->getDx(lay_virtual);
603 lay_pos[1] = alignPar->getDy(lay_virtual);
604 lay_pos[2] = alignPar->getDz(lay_virtual);
605 lay_pos[3] = alignPar->getRx(lay_virtual);
606 lay_pos[4] = alignPar->getRy(lay_virtual);
607 lay_pos[5] = alignPar->getRz(lay_virtual);
608
609
610 if(m_GeoMethod==1)
611 {
612 // six or 3 in m_geoalign
613 m_geoalign->setDx(lay_virtual, lay_pos[0]);
614 m_geoalign->setDy(lay_virtual, lay_pos[1]);
615 m_geoalign->setDz(lay_virtual, lay_pos[2]);
616 m_geoalign->setRx(lay_virtual, lay_pos[3]);
617 m_geoalign->setRy(lay_virtual, lay_pos[4]);
618 m_geoalign->setRz(lay_virtual, lay_pos[5]);
619 m_middriftplane->setAlignment(m_geoalign);
620
621 StraightLine* trk = new StraightLine(trkpar[0], trkpar[1], trkpar[3], trkpar[4]);
622 m_middriftplane->getPointAligned(lay_virtual, trk, posUp, posDown, phiVUp, phiVDown);
623
624 phiVUp[0] = fmod((phiVUp[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
625 phiVDown[0] = fmod((phiVDown[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
626
627 phiVUp_O[0] = fmod((phiVUp_O[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
628 phiVDown_O[0] = fmod((phiVDown_O[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
629 Phi_Meas = fmod((Phi_Meas+2.0*TMath::Pi()),(2.0*TMath::Pi()));
630
631 resi_phi = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (Phi_Meas - phiVUp[0]) : (Phi_Meas - phiVDown[0]);
632 resi_v = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (V_Meas - phiVUp[1]) : (V_Meas - phiVDown[1]);
633
634 resi_phi_O = (fabs(Phi_Meas - phiVUp_O[0])<fabs(Phi_Meas - phiVDown_O[0]))? (Phi_Meas - phiVUp_O[0]) : (Phi_Meas - phiVDown_O[0]);
635 resi_v_O = (fabs(Phi_Meas - phiVUp_O[0])<fabs(Phi_Meas - phiVDown_O[0]))? (V_Meas - phiVUp_O[1]) : (V_Meas - phiVDown_O[1]);
636
637 resi_z = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (Z_Meas - posUp.z()) : (Z_Meas - posDown.z());
638
639 exp_phi = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (phiVUp[0]) : (phiVDown[0]);
640 exp_v = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (phiVUp[1]) : (phiVDown[1]);
641 exp_v = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (phiVUp[1]) : (phiVDown[1]);
642 exp_z = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (posUp.z()) : (posDown.z());
643
644 if(mul_R) resi_phi = resi_phi*m_r[lay_virtual];
645 if(mul_Ste_ang) resi_v = resi_v*sin(Ste_ang[lay_virtual]);
646
647 direction = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? -1 : 1;
648 delete trk;
649
650 }
651
652 // fill histograms
653 m_hresAll_phi->Fill(resi_phi);
654 m_hresAll_v->Fill(resi_v);
655 m_hresAll_z->Fill(resi_z);
656
657 m_hresLay_phi[lay_virtual]->Fill(resi_phi);
658 m_hresLay_v[lay_virtual]->Fill(resi_v);
659 m_hresLay_z[lay_virtual]->Fill(resi_z);
660
661 if (fabs(resi_phi) > m_resiCut_phi[lay_virtual]) continue;
662 if (fabs(resi_v) > m_resiCut_v[lay_virtual]) continue;
663
664 // reset the derivatives arrays
665 m_pMilleAlign -> ZerLoc(&m_derGB_Phi[0], &m_derLC_Phi[0], &m_derNonLin_Phi[0]);
666 m_pMilleAlign -> ZerLoc(&m_derGB_V[0], &m_derLC_V[0], &m_derNonLin_V[0]);
667
668 // derivatives of local parameters
669 for(ipar=0; ipar<m_nloc; ipar++){
670 if( ! getDeriLoc(nTrack_tot,ipar, lay_virtual, Sheet_Meas, direction, Phi_Meas, V_Meas, lay_pos, rechelix, helixErr, deri_Phi, deri_V) ){
671 cout << "getDeriLoc == false!" << setw(12) << itrk << setw(12) << ipar << endl;
672 return false;
673 }
674 if(debug)cout<<"deriLoc : ipar = "<<ipar<<", deri_Phi = "<<deri_Phi<<", deri_V = "<<deri_V<<endl;
675 m_derLC_Phi[ipar] = deri_Phi;
676 m_derLC_V[ipar] = deri_V;
677 }
678
679 // derivatives of global parameters
680 // ipar 0 - 3: dx, dy, dz, rz
681 for(ipar=0; ipar<m_nglo_on_lay; ipar++){
682 iparGB = getAlignParId(lay_virtual, ipar); // I need to use virtual layer here Guoaq/2020-09-08
683 if( ! getDeriGlo(nTrack_tot, ipar, iparGB, lay_virtual, Sheet_Meas, direction, Phi_Meas, V_Meas, lay_pos, helix, helixErr, deri_Phi, deri_V) )
684 {
685 cout << "getDeriGlo == false!" << setw(12) << itrk << setw(12) << ipar << endl;
686 return false;
687 }
688 m_derGB_Phi[iparGB] = deri_Phi;
689 m_derGB_V[iparGB] = deri_V;
690 if(debug)cout<<"deriGlo : ipar, iparGB, layer_virual = "<<ipar<<", "<<iparGB<<", "<<lay_virtual<<", deri_Phi = "<<deri_Phi<<", deri_V = "<<deri_V<<endl;
691 }
692 // Form one equation for each hit, what's that? How about this item? and this guy? guoaq-19-09-14
693
694 if(fabs(Phi_Meas-TMath::Pi())<0.025) continue;
695 if(fabs(Phi_Meas-0)<0.1) continue;
696 if(fabs(Phi_Meas-2.0*TMath::Pi())<0.025) continue;
697 if(fabs(fabs(trkpar[0])-m_r[lay_virtual])<10) continue;
698 if(fabs(fabs(Z_Meas)-m_z[lay_virtual])<5) continue;
699 if(m_derGB_V[getAlignParId(lay_virtual, 3)]<-1000||m_derGB_V[getAlignParId(lay_virtual, 3)]>1000) continue;
700 if(m_derGB_V[getAlignParId(lay_virtual, 4)]<-1000||m_derGB_V[getAlignParId(lay_virtual, 4)]>1000) continue;
701
702 double err_phi = 0;
703 double err_v = 0;
704
705 if(use_phi&&mul_R) err_phi = sig_phi[lay_virtual]*m_r[lay_virtual];
706 if(use_phi&&!mul_R) err_phi = sig_phi[lay_virtual];
707 if(use_v&&mul_Ste_ang) err_v = sig_v[lay_virtual]*sin(Ste_ang[lay_virtual]);
708 if(use_v&&!mul_Ste_ang) err_v = sig_v[lay_virtual];
709
710 m_pMilleAlign -> EquLoc(&m_derGB_Phi[0], &m_derLC_Phi[0], &m_derNonLin_Phi[0], resi_phi, err_phi);
711 m_pMilleAlign -> EquLoc(&m_derGB_V[0], &m_derLC_V[0], &m_derNonLin_V[0], resi_v, err_v);
712 nHitUse++;
713
714 } // loop of nhit
715
716 bool sc = m_pMilleAlign -> FitLoc(m_pMilleAlign->GetTrackNumber(), trkparms, 0);
717 if(sc) {m_pMilleAlign -> SetTrackNumber( m_pMilleAlign->GetTrackNumber()+1 );count[24]++;}
718
719 } // track loop
720
721 nTrack_tot++;
722 return true;
723}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
CgemAliRecTrk::getHelix
HepVector getHelix() const
Definition: CgemAliRecTrk.h:33
CgemAliRecTrk::getStat
int getStat() const
Definition: CgemAliRecTrk.h:27
CgemAliRecTrk::getHelixErr
HepSymMatrix getHelixErr() const
Definition: CgemAliRecTrk.h:34
CgemAlignPar::getDz
double getDz(int iEP)
Definition: CgemAlignPar.h:44
CgemAlignPar::getRy
double getRy(int iEP)
Definition: CgemAlignPar.h:46
CgemAlignPar::getDy
double getDy(int iEP)
Definition: CgemAlignPar.h:43
CgemAlignPar::getDx
double getDx(int iEP)
Definition: CgemAlignPar.h:42
CgemAlignPar::getRz
double getRz(int iEP)
Definition: CgemAlignPar.h:47
CgemAlignPar::getRx
double getRx(int iEP)
Definition: CgemAlignPar.h:45
CgemGeoAlign::setRy
void setRy(int layer_vir, double v)
Definition: CgemGeoAlign.h:88
CgemGeoAlign::setRz
void setRz(int layer_vir, double v)
Definition: CgemGeoAlign.h:89
CgemGeoAlign::setDz
void setDz(int layer_vir, double v)
Definition: CgemGeoAlign.h:86
CgemGeoAlign::setDy
void setDy(int layer_vir, double v)
Definition: CgemGeoAlign.h:85
CgemGeoAlign::setDx
void setDx(int layer_vir, double v)
Definition: CgemGeoAlign.h:84
CgemGeoAlign::setRx
void setRx(int layer_vir, double v)
Definition: CgemGeoAlign.h:87
CgemMidDriftPlane::setAlignment
void setAlignment(CgemGeoAlign *alignPtr)
Definition: CgemMidDriftPlane.h:14
CgemMidDriftPlane::getPointAligned
bool getPointAligned(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
Definition: CgemMidDriftPlane.cxx:164
CgemAlignment::NTRKPAR
const int NTRKPAR
Definition: CgemAlignment.h:17
CgemAlignment::NTRKPARALL
const int NTRKPARALL
Definition: CgemAlignment.h:18

◆ initialize()

void CgemMilleAlign::initialize ( TObjArray *  hlist,
ICgemGeomSvc *  cgemGeomSvc,
ICgemCalibFunSvc *  cgemFunSvc 
)
virtual

Implements CgemAlignBase.

Definition at line 62 of file CgemMilleAlign.cxx.

62 {
63 IMessageSvc* msgSvc;
64 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
65 MsgStream log(msgSvc, "MilleAlign");
66 log << MSG::INFO << "MilleAlign::initialize()" << endreq;
67
68 CgemAlignBase::initialize(hlist, cgemGeomSvc, cgemFunSvc);
69 m_hlist = hlist;
70 m_cgemGeomSvc = cgemGeomSvc;
71 m_cgemFunSvc = cgemFunSvc;
72
73 for(int il=0; il<6; il++)
74 {
75 int ilgeo = int(il/2);
76 m_r[il]=(m_cgemGeomSvc->getCgemLayer(ilgeo)->getMiddleROfGapD());
77 m_z[il]=(m_cgemGeomSvc->getCgemLayer(ilgeo)->getLengthOfCgemLayer()/2.0);
78 Ste_ang[il]=(m_cgemGeomSvc->getCgemLayer(ilgeo)->getAngleOfStereo()*TMath::Pi()/180.0);
79 }
80
81 int DataType = 1;
82
83 if(DataType==1)
84 {
85
86 sig_phi[0] = 0.13/m_r[0];
87 sig_phi[1] = 0.13/m_r[1];
88 sig_phi[2] = 0.13/m_r[2];
89 sig_phi[3] = 0.13/m_r[3];
90 sig_phi[4] = 0.13/m_r[4];
91 sig_phi[5] = 0.13/m_r[5];
92
93 sig_v[0] = 0.13;
94 sig_v[1] = 0.13;
95 sig_v[2] = 0.13;
96 sig_v[3] = 0.13;
97 sig_v[4] = 0.13;
98 sig_v[5] = 0.13;
99
100 }
101 else if(DataType==2)
102 {
103 // for DATA
104 sig_phi[0] = 0.03;
105 sig_phi[1] = 0.03;
106 sig_phi[2] = 0.03;
107 sig_phi[3] = 0.03;
108 sig_phi[4] = 0.03;
109 sig_phi[5] = 0.03;
110
111 sig_v[0] = 1.7;
112 sig_v[1] = 1.7;
113 sig_v[2] = 1.2;
114 sig_v[3] = 1.2;
115 sig_v[4] = 1.2;
116 sig_v[5] = 1.2;
117 }
118
119 for(int i=0; i<30; i++)
120 {
121 count[i] = 0 ;
122 }
123
124 // cuts of derivatives
125 cut_derLC_Phi[0][0] = -0.05;
126 cut_derLC_Phi[1][0] = 0.98;
127 cut_derLC_Phi[2][0] = -0.01;
128 cut_derLC_Phi[3][0] = -0.01;
129 cut_derLC_Phi[4][0] = -0.01;
130
131 cut_derLC_Phi[0][1] = 0.05;
132 cut_derLC_Phi[1][1] = 1.02;
133 cut_derLC_Phi[2][1] = 0.01;
134 cut_derLC_Phi[3][1] = 0.01;
135 cut_derLC_Phi[4][1] = 0.01;
136
137 cut_derLC_V[0][0] = -5.0;
138 cut_derLC_V[1][0] = 60.0;
139 cut_derLC_V[2][0] = -0.01;
140 cut_derLC_V[3][0] = 0.45;
141 cut_derLC_V[4][0] = -90.0;
142
143 cut_derLC_V[0][1] = 5.0;
144 cut_derLC_V[1][1] = 150;
145 cut_derLC_V[2][1] = 0.01;
146 cut_derLC_V[3][1] = 0.80;
147 cut_derLC_V[4][1] = 90.0;
148
149 cut_derGB_Phi[0][0] = -0.05;
150 cut_derGB_Phi[1][0] = -0.02;
151 cut_derGB_Phi[2][0] = -0.01;
152 cut_derGB_Phi[3][0] = -1.02;
153
154 cut_derGB_Phi[0][1] = 0.05;
155 cut_derGB_Phi[1][1] = 0.02;
156 cut_derGB_Phi[2][1] = 0.01;
157 cut_derGB_Phi[3][1] = -0.98;
158
159 cut_derGB_V[0][0] = -5.00;
160 cut_derGB_V[1][0] = -5.00;
161 cut_derGB_V[2][0] = -0.80;
162 cut_derGB_V[3][0] = -150;
163
164 cut_derGB_V[0][1] = 5.00;
165 cut_derGB_V[1][1] = 5.00;
166 cut_derGB_V[2][1] = -0.50;
167 cut_derGB_V[3][1] = -60;
168
169 // initialize hitograms
170 m_hresAll_phi = new TH1F("Resi_All_Phi", "", 200, -2.0, 2.0);
171 m_hlist->Add(m_hresAll_phi);
172
173 m_hresAll_v = new TH1F("Resi_All_v", "", 100, -30.0, 30.0);
174 m_hlist->Add(m_hresAll_v);
175
176 m_hresAll_z = new TH1F("Resi_All_z", "", 100, -30.0, 30.0);
177 m_hlist->Add(m_hresAll_z);
178
179 char hname_phi[200];
180 char hname_v[200];
181 char hname_z[200];
182
183 for(int lay=0; lay<LAYERNMAX; lay++){
184 sprintf(hname_phi, "Res_Phi_Layer%02d", lay);
185 m_hresLay_phi[lay] = new TH1F(hname_phi, "", 200, -2.0, 2.0);
186 m_hlist->Add(m_hresLay_phi[lay]);
187
188 sprintf(hname_v, "Res_V_Layer%02d", lay);
189 m_hresLay_v[lay] = new TH1F(hname_v, "", 100, -30.0, 30.0);
190 m_hlist->Add(m_hresLay_v[lay]);
191
192 sprintf(hname_z, "Res_Z_Layer%02d", lay);
193 m_hresLay_z[lay] = new TH1F(hname_z, "", 100, -30.0, 30.0);
194 m_hlist->Add(m_hresLay_z[lay]);
195 }
196
197 // for debug
198 m_hddoca = new TH1F("delt_doca", "", 200, -1.0, 1.0);
199 m_hlist->Add(m_hddoca);
200
201 char hname[200];
202 for(int lay=0; lay<LAYERNMAX; lay++){
203 sprintf(hname, "delt_docaLay%02d", lay);
204 m_hddocaLay[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
205 m_hlist->Add(m_hddocaLay[lay]);
206 }
207
208 m_GeoMethod = GEO_METHOD;
209 Plot_der = false;
210 debug = false;
211 re_3lay = false;
212 record_pos = true;
213 use_phi = true;
214 use_v = true;
215 mul_R = false;
216 mul_Ste_ang = false;
217
218 nTrack_tot = 0;
219
220 // initialize millepede
221
222 m_nlay = NLAYER;
223 m_nloc = NTRKPAR;
224 m_nglo_on_lay = NDOFALIGN; // NDOFALIGN : for each hit, 4 alignment parameters on each layer: dx, dy, dz, rz
225 m_npar = NDOFALIGN * m_nlay;// that is the global parameters, isn't it?
226
227 if(debug)cout<<"m_nlay, m_nloc, m_nglo_on_lay, m_npar = "<<m_nlay<<", "<<m_nloc<<", "<<m_nglo_on_lay<<", "<<m_npar<<endl;
228
229 m_middriftplane = m_cgemGeomSvc->getMidDriftPtr();
230 m_geoalign = m_cgemGeomSvc->getAlignPtr();
231
232 // get the initial alignment parameter and save it to a list
233 for(int i=0; i<NLAYER; i++)
234 {
235
236 ini_alig_par[i][0] = m_geoalign->getDx(i);
237 ini_alig_par[i][1] = m_geoalign->getDy(i);
238 ini_alig_par[i][2] = m_geoalign->getDz(i);
239 ini_alig_par[i][3] = m_geoalign->getRx(i);
240 ini_alig_par[i][4] = m_geoalign->getRy(i);
241 ini_alig_par[i][5] = m_geoalign->getRz(i);
242
243 }
244
245
246 int i;
247 for(i=0; i<NDOFALIGN; i++){
248 m_dofs[i] = g_dofs[i];
249 for (int j=i*m_nlay; j<(i+1)*m_nlay; j++)
250 {
251 m_sigm[j] = g_Sigm[j];
252 }
253 }
254
255 m_pMilleAlign = new Millepede();
256 m_pMilleAlign -> InitMille(&m_dofs[0], &m_sigm[0], m_nlay,m_nglo_on_lay, m_nloc,
257 g_start_chi_cut, 3, g_res_cut, g_res_cut_init);
258
259 m_derLC_Phi.resize(m_nloc);
260 m_derLC_V.resize(m_nloc);
261 m_derGB_Phi.resize(m_npar);
262 m_derGB_V.resize(m_npar);
263 m_derNonLin_Phi.resize(m_npar);
264 m_derNonLin_V.resize(m_npar);
265 m_par.resize(m_npar);
266 m_error.resize(m_npar);
267 m_pull.resize(m_npar);
268
269 log << MSG::INFO << "MilleAlign::initialize() finished!!" << endreq;
270
271 // contraints
272 std::vector<double> constDY1;
273 std::vector<double> constDY2;
274 std::vector<double> constDY3;
275 std::vector<double> constDY4;
276 std::vector<double> constDY5;
277 std::vector<double> constDY6;
278
279 std::vector<double> constDX5;
280 std::vector<double> constDZ5;
281 std::vector<double> constRX5;
282 std::vector<double> constRY5;
283 std::vector<double> constRZ5;
284
285 std::vector<double> constDX6;
286 std::vector<double> constDZ6;
287 std::vector<double> constRX6;
288 std::vector<double> constRY6;
289 std::vector<double> constRZ6;
290
291 std::vector<double> constDX1;
292 std::vector<double> constDZ1;
293 std::vector<double> constRX1;
294 std::vector<double> constRY1;
295 std::vector<double> constRZ1;
296
297 std::vector<double> constDX2;
298 std::vector<double> constDZ2;
299 std::vector<double> constRX2;
300 std::vector<double> constRY2;
301 std::vector<double> constRZ2;
302
303 std::vector<double> constDX3DX4;
304 std::vector<double> constDZ3DZ4;
305 std::vector<double> constRX3RX4;
306 std::vector<double> constRY3RY4;
307 std::vector<double> constRZ3RZ4;
308
309 std::vector<double> constRX3;
310 std::vector<double> constRY3;
311 std::vector<double> constRX4;
312 std::vector<double> constRY4;
313
314 // fix Dy of all layers
315 constDY1.resize(m_npar);
316 constDY2.resize(m_npar);
317 constDY3.resize(m_npar);
318 constDY4.resize(m_npar);
319 constDY5.resize(m_npar);
320 constDY6.resize(m_npar);
321
322 // fix Dx, Dz, Rz of the outmost layers
323 constDX5.resize(m_npar);
324 constDZ5.resize(m_npar);
325 constRX5.resize(m_npar);
326 constRY5.resize(m_npar);
327 constRZ5.resize(m_npar);
328
329 constDX6.resize(m_npar);
330 constDZ6.resize(m_npar);
331 constRX6.resize(m_npar);
332 constRY6.resize(m_npar);
333 constRZ6.resize(m_npar);
334
335 // fix Dx, Dz, Rz of the midle layers
336 constDX1.resize(m_npar);
337 constDZ1.resize(m_npar);
338 constRX1.resize(m_npar);
339 constRY1.resize(m_npar);
340 constRZ1.resize(m_npar);
341
342 constDX2.resize(m_npar);
343 constDZ2.resize(m_npar);
344 constRX2.resize(m_npar);
345 constRY2.resize(m_npar);
346 constRZ2.resize(m_npar);
347
348 constDX3DX4.resize(m_npar);
349 constDZ3DZ4.resize(m_npar);
350 constRX3RX4.resize(m_npar);
351 constRY3RY4.resize(m_npar);
352 constRZ3RZ4.resize(m_npar);
353
354 constRX3.resize(m_npar);
355 constRY3.resize(m_npar);
356 constRX4.resize(m_npar);
357 constRY4.resize(m_npar);
358
359 for(i=0; i<m_npar; i++){
360
361 constDY1[i] = 0.0;
362 constDY2[i] = 0.0;
363 constDY3[i] = 0.0;
364 constDY4[i] = 0.0;
365 constDY5[i] = 0.0;
366 constDY6[i] = 0.0;
367
368 constDX5[i] = 0.0;
369 constDZ5[i] = 0.0;
370 constRX5[i] = 0.0;
371 constRY5[i] = 0.0;
372 constRZ5[i] = 0.0;
373
374 constDX6[i] = 0.0;
375 constDZ6[i] = 0.0;
376 constRX6[i] = 0.0;
377 constRY6[i] = 0.0;
378 constRZ6[i] = 0.0;
379
380 constDX1[i] = 0.0;
381 constDZ1[i] = 0.0;
382 constRX1[i] = 0.0;
383 constRY1[i] = 0.0;
384 constRZ1[i] = 0.0;
385
386 constDX2[i] = 0.0;
387 constDZ2[i] = 0.0;
388 constRX2[i] = 0.0;
389 constRY2[i] = 0.0;
390 constRZ2[i] = 0.0;
391
392 constDX3DX4[i] = 0.0;
393 constDZ3DZ4[i] = 0.0;
394 constRX3RX4[i] = 0.0;
395 constRY3RY4[i] = 0.0;
396 constRZ3RZ4[i] = 0.0;
397
398 constRX3[i] = 0.0;
399 constRY3[i] = 0.0;
400 constRX4[i] = 0.0;
401 constRY4[i] = 0.0;
402
403 }
404
405 // Fix the position of 2nd and 3rd layer and dy in all 3 layers
406 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
407 // DX DY DZ RX RY RZ
408 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
409 // DX1 DX2 DX3 DX4 DX5 DX6 DY1 DY2 DY3 DY4 DY5 DY6 DZ1 DZ2 DZ3 DZ4 DZ5 DZ6 RX1 RX2 RX3 RX4 RX5 RX6 RY1 RY2 RY3 RY4 RY5 RY6 RZ1 RZ2 RZ3 RZ4 RZ5 RZ6
410 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
411 // index = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
412 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
413
414 // fix the DY on all layers
415 constDY1[6] = 1.0;
416 constDY2[7] = 1.0;
417 constDY3[8] = 1.0;
418 constDY4[9] = 1.0;
419 constDY5[10] = 1.0;
420 constDY6[11] = 1.0;
421
422 // fix the 3rd layer
423 constDX5[4] = 1.0;
424 constDZ5[16] = 1.0;
425 constRX5[22] = 1.0;
426 constRY5[28] = 1.0;
427 constRZ5[34] = 1.0;
428
429 constDX6[5] = 1.0;
430 constDZ6[17] = 1.0;
431 constRX6[23] = 1.0;
432 constRY6[29] = 1.0;
433 constRZ6[35] = 1.0;
434
435 // change to fix the 1st layer
436 constDX1[0] = 1.0;
437 constDZ1[12] = 1.0;
438 constRX1[18] = 1.0;
439 constRY1[24] = 1.0;
440 constRZ1[30] = 1.0;
441
442 constDX2[1] = 1.0;
443 constDZ2[13] = 1.0;
444 constRX2[19] = 1.0;
445 constRY2[25] = 1.0;
446 constRZ2[31] = 1.0;
447
448 // require the Dx3 = Dx4
449 constDX3DX4[2] = 1.0;
450 constDX3DX4[3] = -1.0;
451 // require the Dz3 = Dz4
452 constDZ3DZ4[14] = 1.0;
453 constDZ3DZ4[15] = -1.0;
454 // require the Rx3 = Rx4
455 constRX3RX4[20] = 1.0;
456 constRX3RX4[21] = -1.0;
457 constRX3[20] = 1.0;
458 constRX4[21] = 1.0;
459 // require the Ry3 = Ry4
460 constRY3RY4[26] = 1.0;
461 constRY3RY4[27] = -1.0;
462 constRY3[26] = 1.0;
463 constRY4[27] = 1.0;
464 // require the Rz3 = Rz4
465 constRZ3RZ4[32] = 1.0;
466 constRZ3RZ4[33] = -1.0;
467
468 // apply the constrains
469 m_pMilleAlign -> ConstF(&constDY1[0], 0.0);
470 m_pMilleAlign -> ConstF(&constDY2[0], 0.0);
471 m_pMilleAlign -> ConstF(&constDY3[0], 0.0);
472 m_pMilleAlign -> ConstF(&constDY4[0], 0.0);
473 m_pMilleAlign -> ConstF(&constDY5[0], 0.0);
474 m_pMilleAlign -> ConstF(&constDY6[0], 0.0);
475
476 m_pMilleAlign -> ConstF(&constDX5[0], 0.0);
477 m_pMilleAlign -> ConstF(&constDZ5[0], 0.0);
478 m_pMilleAlign -> ConstF(&constRX5[0], 0.0);
479 m_pMilleAlign -> ConstF(&constRY5[0], 0.0);
480 //m_pMilleAlign -> ConstF(&constRZ5[0], 0.0);
481
482 m_pMilleAlign -> ConstF(&constDX6[0], 0.0);
483 m_pMilleAlign -> ConstF(&constDZ6[0], 0.0);
484 m_pMilleAlign -> ConstF(&constRX6[0], 0.0);
485 m_pMilleAlign -> ConstF(&constRY6[0], 0.0);
486 //m_pMilleAlign -> ConstF(&constRZ6[0], 0.0);
487
488 m_pMilleAlign -> ConstF(&constDX1[0], 0.0);
489 m_pMilleAlign -> ConstF(&constDZ1[0], 0.0);
490 m_pMilleAlign -> ConstF(&constRX1[0], 0.0);
491 m_pMilleAlign -> ConstF(&constRY1[0], 0.0);
492 m_pMilleAlign -> ConstF(&constRZ1[0], 0.0);
493
494 m_pMilleAlign -> ConstF(&constDX2[0], 0.0);
495 m_pMilleAlign -> ConstF(&constDZ2[0], 0.0);
496 m_pMilleAlign -> ConstF(&constRX2[0], 0.0);
497 m_pMilleAlign -> ConstF(&constRY2[0], 0.0);
498 m_pMilleAlign -> ConstF(&constRZ2[0], 0.0);
499
500 m_pMilleAlign -> ConstF(&constDX3DX4[0], 0.0);
501 //m_pMilleAlign -> ConstF(&constDZ3DZ4[0], 0.0);
502
503 //m_pMilleAlign -> ConstF(&constRX3RX4[0], 0.0);
504 //m_pMilleAlign -> ConstF(&constRY3RY4[0], 0.0);
505
506 //m_pMilleAlign -> ConstF(&constRX3[0], 0.0);
507 //m_pMilleAlign -> ConstF(&constRX4[0], 0.0);
508 //m_pMilleAlign -> ConstF(&constRY3[0], 0.0);
509 //m_pMilleAlign -> ConstF(&constRY4[0], 0.0);
510
511 //m_pMilleAlign -> ConstF(&constRZ3RZ4[0], 0.0);
512
513}
CgemAlignBase::initialize
virtual void initialize(TObjArray *hlist, ICgemGeomSvc *cgemGeomSvc, ICgemCalibFunSvc *cgemFunSvc)=0
Definition: CgemAlignBase.cxx:13
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
CgemGeoLayer::getLengthOfCgemLayer
double getLengthOfCgemLayer() const
Definition: CgemGeoLayer.h:22
CgemGeoLayer::getAngleOfStereo
double getAngleOfStereo() const
Definition: CgemGeoLayer.h:218
ICgemGeomSvc::getCgemLayer
virtual CgemGeoLayer * getCgemLayer(int i) const =0
ICgemGeomSvc::getAlignPtr
virtual CgemGeoAlign * getAlignPtr() const =0
ICgemGeomSvc::getMidDriftPtr
virtual CgemMidDriftPlane * getMidDriftPtr() const =0
CgemAlignment::g_res_cut
const double g_res_cut
Definition: CgemAlignment.h:35
CgemAlignment::GEO_METHOD
const int GEO_METHOD
Definition: CgemAlignment.h:14
CgemAlignment::NLAYER
const int NLAYER
Definition: CgemAlignment.h:16
CgemAlignment::g_start_chi_cut
const double g_start_chi_cut
Definition: CgemAlignment.h:37
CgemAlignment::g_res_cut_init
const double g_res_cut_init
Definition: CgemAlignment.h:36
CgemAlignment::g_dofs
const bool g_dofs[6]
Definition: CgemAlignment.h:21
CgemAlignment::NDOFALIGN
const int NDOFALIGN
Definition: CgemAlignment.h:19
CgemAlignment::g_Sigm
const double g_Sigm[36]
Definition: CgemAlignment.h:25

◆ setParam()

void CgemMilleAlign::setParam ( CgemAliParams &  param)
inlinevirtual

Implements CgemAlignBase.

Definition at line 131 of file CgemMilleAlign.h.

131 {
132 CgemAlignBase::setParam(param);
133 m_param = param;
134}
CgemAlignBase::setParam
virtual void setParam(CgemAliParams &param)=0
Definition: CgemAlignBase.h:41

◆ updateAlignPar()

bool CgemMilleAlign::updateAlignPar ( CgemAlignPar *  alignPar)
virtual

Implements CgemAlignBase.

Definition at line 741 of file CgemMilleAlign.cxx.

741 {
742
743 IMessageSvc* msgSvc;
744 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
745 MsgStream log(msgSvc, "MilleAlign");
746 log << MSG::INFO << "MilleAlign::updateConst()" << endreq;
747
748 m_pMilleAlign -> MakeGlobalFit(&m_par[0], &m_error[0], &m_pull[0]);
749
750 int iLAYER;
751 int ipar;
752 double val;
753 double err;
754 for(int i=0; i<NDOFALIGN; i++){ //loop dx,dy,dz,rz,i
755 for(iLAYER=0; iLAYER<NLAYER; iLAYER++){ //loop 3 layer , iLAYER
756 ipar = i * NLAYER + iLAYER;
757
758 if(m_dofs[i]){
759 val = m_par[ipar];
760 err = m_error[ipar];
761 } else{
762 val = 0.0;
763 err = 0.0;
764 }
765 if(0 == i){
766 alignPar->setDelDx(iLAYER, val);
767 alignPar->setErrDx(iLAYER, err);
768 } else if(1 == i){
769 alignPar->setDelDy(iLAYER, val);
770 alignPar->setErrDy(iLAYER, err);
771 } else if(2 == i){
772 alignPar->setDelDz(iLAYER, val);
773 alignPar->setErrDz(iLAYER, err);
774 } else if(3 == i){
775 alignPar->setDelRx(iLAYER, val);
776 alignPar->setErrRx(iLAYER, err);
777 } else if(4 == i){
778 alignPar->setDelRy(iLAYER, val);
779 alignPar->setErrRy(iLAYER, err);
780 } else if(5 == i){
781 alignPar->setDelRz(iLAYER, val);
782 alignPar->setErrRz(iLAYER, err);
783 }
784
785 }
786 }
787 return true;
788}
CgemAlignPar::setErrRy
void setErrRy(int iEP, double val)
Definition: CgemAlignPar.h:39
CgemAlignPar::setDelDy
void setDelDy(int iEP, double val)
Definition: CgemAlignPar.h:29
CgemAlignPar::setErrDy
void setErrDy(int iEP, double val)
Definition: CgemAlignPar.h:36
CgemAlignPar::setErrRx
void setErrRx(int iEP, double val)
Definition: CgemAlignPar.h:38
CgemAlignPar::setDelRy
void setDelRy(int iEP, double val)
Definition: CgemAlignPar.h:32
CgemAlignPar::setDelDz
void setDelDz(int iEP, double val)
Definition: CgemAlignPar.h:30
CgemAlignPar::setDelRz
void setDelRz(int iEP, double val)
Definition: CgemAlignPar.h:33
CgemAlignPar::setErrDz
void setErrDz(int iEP, double val)
Definition: CgemAlignPar.h:37
CgemAlignPar::setErrDx
void setErrDx(int iEP, double val)
Definition: CgemAlignPar.h:35
CgemAlignPar::setErrRz
void setErrRz(int iEP, double val)
Definition: CgemAlignPar.h:40
CgemAlignPar::setDelDx
void setDelDx(int iEP, double val)
Definition: CgemAlignPar.h:28
CgemAlignPar::setDelRx
void setDelRx(int iEP, double val)
Definition: CgemAlignPar.h:31
The documentation for this class was generated from the following files: