FairRoot/PandaRoot
SQM.h
Go to the documentation of this file.
1 void PndFTSCAGBTracker::Refit_1( FTSCANPletV& triplet, const FTSCAHits& hits )
2 {
3 
4 //#include "Tang.h"
5 
6 
7 Double_t det;
8 TMatrixD H(4,4);
9 TMatrixD HY(4,1);
10 TMatrixD Cov(4,4);
11 TMatrixD HC(4,1);
12 
13 float f1, f2, f3, f4, f5, Y, chi2;
14 
15 for (Int_t i=0; i<4; i++){
16  HY(i,0)=0;
17  HC(i,0)=0;
18  for (Int_t j=0; j<4; j++){
19  H(i,j)=0;
20  }
21 }
22 
23  const int N = triplet.N();
24 
25 float Tx, Ty, Bx, By;
26 
27 //#include "SQM_0.h"
28 
29 cout << " **************** Param before Refit_1 ************ " << endl;
30 
31  PndFTSCATrackParamVector &param = triplet.Param();
32 
33  cout << " param.Tx()[0] " << param.Tx()[0]
34  << " param.Ty()[0] " << param.Ty()[0]
35  << " param.QP()[0] " << param.QP()[0]
36  << " param.X()[0] " << param.X()[0]
37  << " param.Y()[0] " << param.Y()[0]
38  << " param.Z()[0] " << param.Z()[0]
39  << " param.Chi2()[0] " << param.Chi2()[0]
40  << " param.NDF()[0] " << param.NDF()[0] << endl;
41 
42 Tx = param.Tx()[0];
43 Ty = param.Ty()[0];
44 Bx = param.X()[0] - param.Tx()[0]*param.Z()[0];
45 By = param.Y()[0] - param.Ty()[0]*param.Z()[0];
46 
47 float Tx_K = Tx;
48 float Ty_K = Ty;
49 float Bx_K = Bx;
50 float By_K = By;
51 
52 cout << " **************** Before Refit_1 ************ " << endl;
53 cout << " Tx " << Tx << endl;
54 cout << " Ty " << Ty << endl;
55 cout << " Bx " << Bx << endl;
56 cout << " By " << By << endl;
57 
58 // float_m active = triplet.IsValid();
59 
60  for ( int it = 0; it < 3; it++ ) {
61 
62 for (Int_t i=0; i<4; i++){
63  HY(i,0)=0;
64  HC(i,0)=0;
65  for (Int_t j=0; j<4; j++){
66  H(i,j)=0;
67  }
68 }
69 
70  for ( int ihit = 0; ihit < N; ihit++ ) {
71  const TESV& index = triplet.IHit(ihit);
72 cout << "index.s station " << index.s << " index.e " << index.e << endl;
73  const FTSCAElementsOnStation<FTSCAHit> &hits_st = hits.OnStationConst(index.s[0]);
74 
75  const FTSCAHit &hit = hits_st[index.e[0]];
76 
77  int ISta = (int) hit.IStation();
78  cout << "chtck hit.IStation()" << ISta << endl;
79 
80 //float fX = param.X()[0] - param.Tx()[0]*(param.Z()[0] - hit.X0());
81 //float fY = param.Y()[0] - param.Ty()[0]*(param.Z()[0] - hit.X0());
82 
83 
84 
85  const PndFTSCAGBHit &hit_1 = fHits[hit.Id()];
86 
87 cout << " ************************************************* " << endl;
88 cout << " MC_Tx " << hit_1.point_Px/hit_1.point_Pz << endl;
89 cout << " MC_Ty " << hit_1.point_Py/hit_1.point_Pz << endl;
90 cout << " MC_Bx " << ( hit_1.point_X - (hit_1.point_Px/hit_1.point_Pz)*hit_1.point_Z) << endl;
91 cout << " MC_By " << (hit_1.point_Y - (hit_1.point_Py/hit_1.point_Pz)*hit_1.point_Z) << endl;
92 
93 cout << " hit.R() " << hit.R() << " hit.ErrR() " << sqrt(hit.Err2R()) << endl;
94 cout << " X " << hit.X1() << " Y " << hit.X2() << " Z " << hit.X0() << endl;
95 
96 
97 
98 const FTSCAStation &station = GetParameters().Station( ISta );
99  const float s = station.f.sin;
100 
101  const float c = station.f.cos;
102 
103 cout << " s " << s << " c " << c << endl;
104 
105 float Tx_MC = hit_1.point_Px/hit_1.point_Pz;
106 float Ty_MC = hit_1.point_Py/hit_1.point_Pz;
107 float Bx_MC = ( hit_1.point_X - (hit_1.point_Px/hit_1.point_Pz)*hit_1.point_Z);
108 float By_MC = (hit_1.point_Y - (hit_1.point_Py/hit_1.point_Pz)*hit_1.point_Z);
109  float R_MC = (c*(Tx_MC*hit.X0()+ Bx_MC - hit.X1())
110 + s*(Ty_MC*hit.X0()+ By_MC - hit.X2()) )/sqrt(1.f + (c*Tx_MC + s*Ty_MC)*(c*Tx_MC + s*Ty_MC));
111 
112 if ( s==0)
113 {
114 
115 float D2_mc_hit = sqrt ((hit_1.point_X - hit.X1())*(hit_1.point_X - hit.X1()) + (hit_1.point_Z - hit.X0())*(hit_1.point_Z - hit.X0()));
116 float Diff_D_R = D2_mc_hit -hit.R();
117 cout << " D2_mc_hit " << D2_mc_hit << " hit.R() " << hit.R() <<
118 " Diff_D_R " << Diff_D_R << endl;
119 
120 }
121 
122 float diff_R;
123 
124 if ( R_MC > 0.f)
125 diff_R = R_MC - hit.R();
126 else
127 diff_R = R_MC + hit.R();
128 
129 cout << " diff_R " << diff_R << endl;
130 
131 float R_MC_1 = (c*(Tx_MC*hit_1.point_Z + Bx_MC - hit.X1())
132  + s*(Ty_MC*hit_1.point_Z+ By_MC - hit.X2()) )/sqrt(1.f + (c*Tx_MC + s*Ty_MC)*(c*Tx_MC + s*Ty_MC));
133 
134 float diff_R_1;
135 if ( R_MC_1 > 0.f)
136  diff_R_1 = R_MC_1 - hit.R();
137 else
138  diff_R_1 = R_MC_1 + hit.R();
139 
140 cout << " diff_R_1 " << diff_R_1 << endl;
141 
142 float R_MC_2 = (c*(hit_1.point_X - hit.X1())
143  + s*(hit_1.point_Y - hit.X2()) )/sqrt(1.f + (c*Tx_MC + s*Ty_MC)*(c*Tx_MC + s*Ty_MC));
144 
145 cout << " R_MC_2 " << R_MC_2 << " hit.R()" << hit.R() << endl;
146 
147 float diff_R_2;
148 if (R_MC_2<0.f)
149  diff_R_2 = R_MC_2 + hit.R();
150 else
151 diff_R_2 = R_MC_2 - hit.R();
152 
153 cout << " diff_R_2 " << diff_R_2 << endl;
154 
155 
156 
157 // dbg
158 /*
159 if ( ihit==0 ) {
160 Tx = hit_1.point_Px/hit_1.point_Pz;
161 Ty = hit_1.point_Py/hit_1.point_Pz;
162 Bx = hit_1.point_X - (hit_1.point_Px/hit_1.point_Pz)*hit_1.point_Z;
163 By = hit_1.point_Y - (hit_1.point_Py/hit_1.point_Pz)*hit_1.point_Z;
164 }
165 
166 */
167 
168 float Y_point_of_wire=0.f;
169 //dbg
170 
171 
172 float fX = Tx*hit.X0() + Bx;
173 float fY = Ty*hit.X0() + By;
174 
175 cout << " Tx " << Tx << " Bx " << Bx << " Ty " << Ty << " By " << By << endl;
176 
177 
178 float tx = c*Tx + s*Ty;
179 float rCorrection = sqrt(1.f+tx*tx);
180 
181 float Diff = c*(fX - hit.X1()) + s*(fY - hit.X2());
182 //float Diff = c*(fX - hit.X1()) + s*(fY - Y_point_of_wire);
183 
184 float Diff_mc = c*(hit_1.point_X -hit.X1()) + s*(hit_1.point_Y - hit.X2());
185 
186 //float Diff_mc = c*(hit_1.point_X -hit.X1()) + s*(hit_1.point_Y - Y_point_of_wire);
187 
188 cout << " Diff " << Diff << " Diff_mc " << Diff_mc << endl;
189 
190 cout << " ************************************************* " << endl;
191 
192 float R = hit.R();
193 
194 if (Diff < 0.f ) R=(-1.f)*R;
195 
196 //if (Diff_mc < 0.f ) R=(-1.f)*R; //dbg
197 
198 //float err = (sqrt(hit.Err2R())+ s*sqrt(hit.Err2X2()))*rCorrection;
199 
200 
201 float err = sqrt(hit.Err2R());
202 //float err1= 2*hit.R()*err;
203 
204 Y=(R-(Diff/rCorrection))/err;
205 //Y = (R*R - Diff*Diff/(rCorrection*rCorrection) )/err;
206 
207 f1 = (c*hit.X0()/rCorrection - c*Diff*tx/(rCorrection*rCorrection*rCorrection))/err;
208 
209 
210 
211 
212 //float rCorrection2= rCorrection*rCorrection;
213 //f1 = (2*Diff*c*hit.X0()/rCorrection2 - (Diff*Diff*2*tx*c)/(rCorrection2*rCorrection2))/err;
214 
215 
216 f2 = c/(rCorrection*err);
217 //f2 = 2*Diff*c/(rCorrection*rCorrection*err);
218 
219 f3= (s*hit.X0()/rCorrection - s*Diff*tx/(rCorrection*rCorrection*rCorrection))/err;
220 //f3 = (2*Diff*s*hit.X0()/rCorrection2- (Diff*Diff*2*tx*s)/(rCorrection2*rCorrection2))/err;
221 
222 
223 f4 = s/(rCorrection*err);
224 
225 //f4 =2*Diff*s/(rCorrection*rCorrection*err);
226 
227 //f4=2*Diff*s/(rCorrection*rCorrection*err);
228 
229 
230 
231  H(0,0) = H(0,0) + f1*f1;
232  H(0,1) = H(0,1) + f1*f2;
233  H(0,2) = H(0,2) + f1*f3;
234  H(0,3) = H(0,3) + f1*f4;
235  H(1,0) = H(0,1);
236  H(1,1) = H(1,1) + f2*f2;
237  H(1,2) = H(1,2) + f2*f3;
238  H(1,3) = H(1,3) + f2*f4;
239  H(2,0) = H(0,2);
240  H(2,1) = H(1,2);
241  H(2,2) = H(2,2) + f3*f3;
242  H(2,3) = H(2,3) + f3*f4;
243  H(3,0) = H(0,3);
244  H(3,1) = H(1,3);
245  H(3,2) = H(2,3);
246  H(3,3) = H(2,3) + f4*f4;
247 
248  HY(0,0) = HY(0,0) + Y*f1;
249  HY(1,0) = HY(1,0) + Y*f2;
250  HY(2,0) = HY(2,0) + Y*f3;
251  HY(3,0) = HY(3,0) + Y*f4;
252 
253 
254 } //ihit
255 
256  H.Print();
257  H.Invert(&det);
258  Cov=H;
259  H.Print();
260  Cov.Print();
261  cout << "det= " << det << endl;
262 
263  HC=H*HY;
264  HC.Print();
265 
266 
267 
268 cout << "DTx= " << HC(0,0) << " +- " << sqrt(Cov(0,0)) << endl;
269 cout << "DBx= " << HC(1,0) << " +- " << sqrt(Cov(1,1)) << endl;
270 cout << "DTy= " << HC(2,0) << " +- " << sqrt(Cov(2,2)) << endl;
271 cout << "DBy= " << HC(3,0) << " +- " << sqrt(Cov(3,3)) << endl;
272 
273 
274  Tx= Tx + HC(0,0);
275  Bx= Bx + HC(1,0);
276  Ty= Ty + HC(2,0);
277  By= By + HC(3,0);
278 
279 cout << "Tx= " << Tx << endl;
280 cout << "Ty= " << Ty << endl;
281 cout << "Bx= " << Bx << endl;
282 cout << "By= " << By << endl;
283 
284 
285 } //it
286 
287 param.SetTx(Tx);
288 param.SetTy(Ty);
289 param.SetX(Tx*param.Z()[0] + Bx);
290 param.SetY(Ty*param.Z()[0] + By);
291 param.SetBx(Bx);
292 //float Err2_Tx=(float)Cov(0,0);
293 //param.SetCov( 5, Err2_Tx);
294 
295 
296 
297 float Tx_1 = param.Tx()[0];
298 float Ty_1 = param.Ty()[0];
299 float Bx_1 = param.Bx()[0];
300 float By_1 = By;
301 
302 
303  chi2=0.0;
304 float chi2_k=0.0;
305  for ( int ihit = 0; ihit < N; ihit++ ) {
306  const TESV& index = triplet.IHit(ihit);
307 cout << "index.s station " << index.s << " index.e " << index.e << endl;
308  const FTSCAElementsOnStation<FTSCAHit> &hits_st = hits.OnStationConst(index.s[0]);
309 
310  const FTSCAHit &hit = hits_st[index.e[0]];
311 
312  int ISta = (int) hit.IStation();
313  cout << "chtck hit.IStation()" << ISta << endl;
314 
315 const FTSCAStation &station = GetParameters().Station( ISta );
316  const float s_1 = station.f.sin;
317 
318  const float c_1 = station.f.cos;
319 
320 cout << " s_1 " << s_1 << " c_1 " << c_1 << endl;
321 
322 float tx_1 = c_1*Tx_1 + s_1*Ty_1;
323 float rCorrection_1 = sqrt(1.f+tx_1*tx_1);
324 float fX_1 = Tx_1*hit.X0() + Bx_1;
325 float fY_1 = Ty_1*hit.X0() + By_1;
326 
327 float Diff_1 = c_1*(fX_1 - hit.X1()) + s_1*(fY_1 - hit.X2());
328 
329 float tx_k = c_1*Tx_K + s_1*Ty_K;
330 float rCorrection_k = sqrt(1.f+tx_k*tx_k);
331 float fX_k = Tx_K*hit.X0() + Bx_K;
332 float fY_k = Ty_K*hit.X0() + By_K;
333 
334 float Diff_k = c_1*(fX_k - hit.X1()) + s_1*(fY_k - hit.X2());
335 float R_k = hit.R();
336 if (Diff_k < 0.f ) R_k=(-1.f)*R_k;
337 //float err_k = (sqrt(hit.Err2R())+ s_1*sqrt(hit.Err2X2()) )*rCorrection_k;
338 
339 
340 
341 float Y_point_of_wire_1=0.f;
342 //float Diff_1 = c_1*(fX_1 - hit.X1()) + s_1*(fY_1 - Y_point_of_wire_1);
343 float R_1 = hit.R();
344 if (Diff_1 < 0.f ) R_1=(-1.f)*R_1;
345 
346 //float err_1 = (sqrt(hit.Err2R())+ s_1*sqrt(hit.Err2X2()) )*rCorrection_1;
347 float err_1 = sqrt(hit.Err2R());
348 float err_k =err_1;
349 
350 float Mes, Mes_k;
351 if (ISta > 50 )
352 {
353 Mes = (Diff_1/rCorrection_1 - R_1)/(sqrt(hit.Err2R())+ s_1*sqrt(hit.Err2X2()) + c_1*sqrt(hit.Err2X1()) );
354 Mes_k = (Diff_k/rCorrection_k - R_k)/(sqrt(hit.Err2R())+ s_1*sqrt(hit.Err2X2()) + c_1*sqrt(hit.Err2X1()) );
355 }
356 else
357 {
358 // Mes = (Diff_1/rCorrection_1 - R_1)/(sqrt(hit.Err2R())+ s_1*sqrt(hit.Err2X2()) );
359 //Mes_k = (Diff_k/rCorrection_k - R_k)/(sqrt(hit.Err2R())+ s_1*sqrt(hit.Err2X2()) );
360 Mes = (Diff_1/rCorrection_1 - R_1)/sqrt(hit.Err2R());
361 Mes_k = (Diff_k/rCorrection_k - R_k)/sqrt(hit.Err2R());
362 
363 }
364 float Mes1 = (Diff_1 - R_1*rCorrection_1)/(err_1*rCorrection_1);
365 if ( Mes > 10.f)
366 cout << " change sighn for hit " << ihit << endl;
367 
368 cout << " Mes " << Mes << " Mes1 " << Mes1 << endl;
369 chi2=chi2 + Mes*Mes;
370 chi2_k=chi2_k + Mes_k*Mes_k;
371 cout << " Mes*Mes " << Mes*Mes << endl;
372 
373 }
374 cout << " chi2 " << chi2 << " chi2_k " << chi2_k << endl;
375 
376 param.SetChi2(chi2);
377 
378 }
FTSCAHit::X2
float X2() const
Definition: FTSCAHits.h:43
TESV::s
int_v s
Definition: FTSCATES.h:42
FTSCANPletV::IHit
const TESV & IHit(int IH) const
Definition: FTSCANPletsV.h:53
PndFTSCAParam::Station
const FTSCAStation & Station(short i) const
Definition: PndFTSCAParam.h:45
i
Int_t i
Definition: run_full.C:25
FTSCAStation::f
FTSCAStripInfo f
Definition: FTSCAStation.h:32
FTSCAHit::IStation
char IStation() const
Definition: FTSCAHits.h:33
f1
TF1 * f1
Definition: reco_analys2.C:50
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
s
TLorentzVector s
Definition: Pnd2DStar.C:50
fX
Double_t fX
Definition: PndCaloDraw.cxx:34
TESV
Definition: FTSCATES.h:28
f4
TFile * f4
Definition: ReadMesDecCorrRes.C:9
Y
double Y
Definition: anaLmdDigi.C:68
FTSCAHits::OnStationConst
const FTSCAElementsOnStation< T > & OnStationConst(char i) const
Definition: FTSCAHits.h:134
FTSCAHit::X0
float X0() const
Definition: FTSCAHits.h:41
FTSCAHit::X1
float X1() const
Definition: FTSCAHits.h:42
FTSCANPletV::N
int N() const
Definition: FTSCANPletsV.h:51
PndFTSCAGBTracker::fHits
PndFTSResizableArray< PndFTSCAGBHit > fHits
Definition: PndFTSCAGBTracker.h:181
FTSCAHit::Err2X2
float Err2X2() const
Definition: FTSCAHits.h:52
f3
TFile * f3
Definition: ReadMesDecCorrRes.C:8
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndFTSCAGBTracker::GetParameters
const PndFTSCAParam & GetParameters() const
Definition: PndFTSCAGBTracker.h:73
f
TFile * f
Definition: bump_analys.C:12
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
TESV::e
uint_v e
Definition: FTSCATES.h:43
fY
Double_t fY
Definition: PndCaloDraw.cxx:34
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
FTSCAHit::Id
int Id() const
Definition: FTSCAHits.h:35
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
FTSCAHit::Err2X1
float Err2X1() const
Definition: FTSCAHits.h:50
FTSCANPletV::Param
const PndFTSCATrackParamVector & Param() const
Definition: FTSCANPletsV.h:55
R
Double_t R
Definition: checkhelixhit.C:61
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
PndFTSCAGBTracker::Refit_1
void Refit_1(FTSCANPletV &triplet, const FTSCAHits &hits)
Definition: SQM.h:1
Generated on Wed Apr 3 2019 10:03:02 for FairRoot/PandaRoot by   doxygen 1.8.5