FairRoot/PandaRoot
CbmGeaneTrKalStt.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- FairGeaneTrKalStt source file -----
3 // ----- Created 28/09/06 by M. Al-Turany -----
4 // -------------------------------------------------------------------------
5 
6 
7 #include "TClonesArray.h"
8 #include "FairRootManager.h"
9 #include "FairGeaneTrKalStt.h"
10 #include "TGeant3TGeo.h"
11 #include "TGeant3.h"
12 #include "TVector3.h"
13 #include "CbmPlanePoint.h"
14 #include "TTree.h"
15 #include "TDatabasePDG.h"
16 #include "PndSttTrack.h"
17 #include "PndSttHit.h"
18 #include "PndSttPoint.h"
19 #include "DetPlane.h"
20 #include "FairRunAna.h"
21 #include "PndSttHelixTrackFitter.h"
22 #include "FairTrackParH.h"
23 #include "FairTrackParP.h"
24 #include "TCanvas.h"
25 #include "TH1F.h"
26 
27 #include <iostream>
28 
29 using namespace std;
30 
31 // TCanvas *can2;
32 // TH1F *dstvtx, *momvtx;
33 Int_t eventnum3;
34 // ----- Default constructor -------------------------------------------
35 FairGeaneTrKalStt::FairGeaneTrKalStt() :
36  FairTask("Test") { }
37 // -------------------------------------------------------------------------
38 
39 
40 
41 // ----- Destructor ----------------------------------------------------
42 FairGeaneTrKalStt::~FairGeaneTrKalStt() { }
43 // -------------------------------------------------------------------------
44 
45 
46 // ----- Public method Init --------------------------------------------
47 InitStatus FairGeaneTrKalStt::Init() {
48 
49  eventnum3 = 0;
50 
51  // Get RootManager
52  FairRootManager* ioman = FairRootManager::Instance();
53  if ( ! ioman ) {
54  cout << "-E- FairGeaneTrKalStt::Init: "
55  << "RootManager not instantised!" << endl;
56  return kFATAL;
57  }
58 
59  welldone = 0; notdone = 0;
60  notconsidered = 0;
61  notbackprocessedhit = 0; notprocessedhit = 0;
62  tothits = 0; total = 0;
63 
64  // can2 = new TCanvas();
65  // can2->Divide(1,2);
66  // dstvtx = new TH1F("dstvtx","dstvtx",100,0,10);
67  // momvtx = new TH1F("momvtx","momvtx",100,0,0.5);
68 
69 
70  fTrackParGeane = new TClonesArray("FairTrackParP");
71  ioman->Register("GeaneTrackPar","Geane", fTrackParGeane, kTRUE);
72 
73  fTrackParIni = new TClonesArray("FairTrackParP");
74  ioman->Register("GeaneTrackIni","Geane", fTrackParIni, kTRUE);
75 
76  fTrackParFinal = new TClonesArray("FairTrackParP");
77  ioman->Register("GeaneTrackFinal","Geane", fTrackParFinal, kTRUE);
78 
79  // Get hit Array
80  fHitArray = (TClonesArray*) ioman->GetObject("STTHit");
81  if (!fHitArray) cout << "-W- FairGeaneTrKalStt::Init: No Hit array!" << endl;
82  // Get point Array
83  fPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
84  if (!fPointArray) cout << "-W- FairGeaneTrKalStt::Init: No Point array!" << endl;
85 
86  // Get SttTrack array
87  fTrackArray = (TClonesArray*) ioman->GetObject("STTTrack");
88  if (!fTrackArray) {
89  cout << "-E- FairGeaneTrKalStt::Init: No SttTrack array!" << endl;
90  return kERROR;
91  }
92 
93  fPro = new FairGeanePro();
94  fUtil = new FairGeaneUtil();
95  // fPro->PropagateToVolume("mL3",0,1);
96 
97 
98  return kSUCCESS;
99 
100 }
101 // -------------------------------------------------------------------------
102 
103 void FairGeaneTrKalStt::FinishTask() {
104  cout << "FairGeaneTrKalStt::FinishTask" << " total " << total << " total hits " << tothits << endl;
105  cout << "WELL DONE " << welldone << " NOT DONE " << notdone << endl;
106  cout << "notconsidered " << notconsidered << " notprocessedhit " << notprocessedhit << " notbackprocessedhit " << notbackprocessedhit << endl;
107 }
108 
109 // ----- Public method Exec --------------------------------------------
110 void FairGeaneTrKalStt::Exec(Option_t*) {
111 
112  // cout << "EVENT " << eventnum3 << endl;
113  eventnum3++;
114 
115  if(total > 900) {
116  cout << "FairGeaneTrKalStt::Exec" << " total " << total << " total hits " << tothits << endl;
117  cout << "WELL DONE " << welldone << " NOT DONE " << notdone << endl;
118  cout << "notconsidered " << notconsidered << " notprocessedhit " << notprocessedhit << " notbackprocessedhit " << notbackprocessedhit << endl;
119  }
120 
121  fTrackParGeane->Delete();
122  fTrackParIni->Delete();
123  fTrackParFinal->Delete();
124 
125  if (!fTrackArray ) return;
126 
127 
128  Int_t nTracks = fTrackArray->GetEntriesFast();
129 
130  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
131  {
132 
133  PndSttTrack* pTrack = (PndSttTrack*) fTrackArray->At(iTrack);
134  if(!pTrack) continue;
135 
136  total++;
137 
138 
139  if(pTrack->GetFlag() < 3) {
140  notconsidered++;
141  continue;
142  } // per ora voglio solo le tracce perfette
143 
144  // STARTING POINT
145 
146  /*
147  Int_t iHit1 = pTrack->GetHitIndex(0);
148  PndSttHit *currenthit1 = (PndSttHit*) fHitArray->At(iHit1);
149  if(!currenthit1) continue;
150  if(currenthit1->GetXint() == -999 || currenthit1->GetYint() == -999) continue;
151  Int_t refindex1 = currenthit1->GetRefIndex();
152  // get point
153  PndSttPoint *fPoint1 = (PndSttPoint*) fPointArray->At(refindex1);
154  StartPos = TVector3 (fPoint1->GetXtot(),fPoint1->GetYtot(),fPoint1->GetZtot());
155  StartPosErr = TVector3(0.0,0.0,0.0);
156  StartMom = TVector3 (fPoint1->GetPXtot(),fPoint1->GetPYtot(),fPoint1->GetPZtot());
157  StartMomErr = TVector3(0.0,0.0,0.0);
158  */
159 
160  // FROM VERTEX
161  Bool_t vtx = RetrieveVertex(pTrack);
162 
163  // traccio fino all' ultimo hit della traccia PCA al wire
164  Int_t hitcounter = pTrack->GetNofHits();
165  tothits += hitcounter;
166 
167  Int_t iHit2 = pTrack->GetHitIndex(0);
168  PndSttHit *currenthit2 = (PndSttHit*) fHitArray->At(iHit2);
169  if(!currenthit2) continue;
170  if(currenthit2->GetXint() == -999 || currenthit2->GetYint() == -999) continue;
171  Int_t refindex2 = currenthit2->GetRefIndex();
172  // get point
173  PndSttPoint *iPoint2 = (PndSttPoint*) fPointArray->At(refindex2);
174  // if(TMath::Sqrt(iPoint2->GetXInLocal()*iPoint2->GetXInLocal() + iPoint2->GetYInLocal()*iPoint2->GetYInLocal()) < 0.4998) continue;
175  // cout << "************************************" << endl;
176  // cout << "HIT CENTER " << currenthit2->GetX() << " " << currenthit2->GetY() << " " << currenthit2->GetZ() << endl;
177 
178  // out
179  // let' s take (approx!) the medium point
180  // position --> (posout - posin) /2
181  // momentum --> (momout - momin) / 2
182  TVector3 EndPos = TVector3(iPoint2->GetXtot(), iPoint2->GetYtot(), iPoint2->GetZtot());
183  TVector3 EndPosErr = TVector3(0,0,0);
184  TVector3 EndMom = TVector3(iPoint2->GetPXtot(), iPoint2->GetPYtot(), iPoint2->GetPZtot());
185  TVector3 EndMomErr = TVector3(0,0,0);
186 
187  PDGCode= 13;
188  TDatabasePDG *fdbPDG= TDatabasePDG::Instance();
189  TParticlePDG *fParticle= fdbPDG->GetParticle(PDGCode);
190  Double_t fCharge= fParticle->Charge();
191 
192  TClonesArray& clref1 = *fTrackParIni;
193  Int_t size1 = clref1.GetEntriesFast();
194  FairTrackParP *fStart = new (clref1[size1]) FairTrackParP(StartPos, StartMom, StartPosErr, StartMomErr, fCharge, StartPos, TVector3(0.,1.,0.), TVector3(0.,0.,1.));
195 
196  TClonesArray& clref = *fTrackParGeane;
197  Int_t size = clref.GetEntriesFast();
198  FairTrackParP *fRes = new(clref[size]) FairTrackParP();
199 
200  TClonesArray& clref2 = *fTrackParFinal;
201  Int_t size2 = clref2.GetEntriesFast();
202  FairTrackParP *fFinal = new(clref2[size2]) FairTrackParP();
203  FairTrackParH *fFinalH = new FairTrackParH(EndPos, EndMom, EndPosErr, EndMomErr, fCharge);
204 
205  Int_t processedhit = 0;
206  Int_t backprocessedhit = 0;
207 
208  FairTrackParP *fRunningRes = new FairTrackParP();
209 
210  for(Int_t iteration = 0; iteration < 3; iteration++)
211  {
212  // if(StartMom.Mag() < 0.100) continue;
213 
214  // cout << "ITERATION " << iteration << endl;
215  // running start point: ogni volta viene riaggiornato col punto di kalman
216  FairTrackParP *fRunningStart = new FairTrackParP(StartPos, StartMom, StartPosErr, StartMomErr, fCharge, TVector3(0.,0.,0.), TVector3(0.,1.,0.), TVector3(0.,0.,1.));
217  if(!fRunningStart) continue;
218  fRunningRes->Reset();
219 
220  // cout << "START " << backprocessedhit << " " << hitcounter << endl;
221  // propagationo ok
222  Int_t donecounter = 0;
223  processedhit = 0;
224  // loop on hits
225  // loop on hits
226  Int_t start = 0;
227  if(iteration > 0) start = backprocessedhit + 1;
228  for(int i = start; i < hitcounter; i++)
229  // for(int i = backprocessedhit; i < hitcounter; i++)
230  // for(int i = 1; i < hitcounter; i++)
231  {
232 
233  // cout << "i " << i << endl;
234 
235  Bool_t rc_prochit = ProcessHit(pTrack, i, fRunningStart, fRunningRes, "FWD");
236  if(rc_prochit == kFALSE) {
237  notprocessedhit++;
238  continue;
239  }
240  processedhit = i;
241  donecounter++;
242  }
243  if(donecounter == 0) {
244  notdone++;
245  continue; // se nessun estrap/kalman e' andato bene butto la traccia (o tengo il prefit? CHECK!)
246  }
247  // BACK TRACKING AL PRIMO PUNTO
248  // cout << "STARTING BACK TRACKING" << donecounter << endl;
249  FairTrackParP *fRunningKal = new FairTrackParP();
250  Double_t KD2[15];
251  fRunningRes->GetCov(KD2);
252  fRunningKal->SetTrackPar(fRunningRes->GetV(), fRunningRes->GetW(),
253  fRunningRes->GetTV(), fRunningRes->GetTW(),
254  fRunningRes->GetQp(), KD2, fRunningRes->GetOrigin(),
255  fRunningRes->GetIVer(), fRunningRes->GetJVer(), fRunningRes->GetKVer(),
256  fRunningRes->GetSPU());
257  Int_t backcounter = 0;
258  backprocessedhit = 0;
259  // loop on hits
260  for(int i = (processedhit-1); i >= 0; i--)
261  // for(int i = (donecounter-2); i >= 0; i--)
262  {
263  // cout << "i " << i << endl;
264  Bool_t rc_back = ProcessHit(pTrack, i, fRunningKal, fRunningRes, "BKW");
265  if(rc_back == kFALSE) {
266  notbackprocessedhit++;
267  continue;
268  }
269  backprocessedhit = i;
270  backcounter++;
271  StartPos = TVector3(fRunningRes->GetX(), fRunningRes->GetY(), fRunningRes->GetZ());
272  StartMom = TVector3(fRunningRes->GetPx(), fRunningRes->GetPy(), fRunningRes->GetPz());
273  StartPosErr = TVector3(fRunningRes->GetDX(), fRunningRes->GetDY(), fRunningRes->GetDZ());
274  StartMomErr = TVector3(fRunningRes->GetDPx(), fRunningRes->GetDPy(), fRunningRes->GetDPz());
275  }
276  if(backcounter == 0) {
277  notdone++;
278  continue;
279  }
280 
281 
282 
283  // cout << "MOME " << StartMom.Mag() << endl;
284  }
285 
286  //------------------------------------------------------
287  // BACK TRACK TO VERTEX
288  Bool_t rc_vtx = BackToVertex(fRunningRes, fRes);
289  if(rc_vtx == kFALSE) {
290  notdone++;
291  continue;
292  }
293  welldone++;
294 
295 
296  if(!fRes) continue;
297 
298  Double_t PD[3], RD[15], H[3], CH, SPU, DJ[3], DK[3], PC[3], RC[15];
299  Int_t IERR;
300 
301  PC[0] = fFinalH->GetQp()/fFinalH->GetQ();
302  PC[1] = fFinalH->GetLambda();
303  PC[2] = fFinalH->GetPhi();
304 
305  fFinalH->GetCov(RC);
306 
307  H[0] = 0; H[1] = 0; H[2] = 20;
308 
309  if(fRunningRes->GetQ()!=0) CH = fRunningRes->GetQ()/TMath::Abs(fRunningRes->GetQ());
310  else continue;
311 
312  DJ[0] = fRunningRes->GetJVer().X();
313  DJ[1] = fRunningRes->GetJVer().Y();
314  DJ[2] = fRunningRes->GetJVer().Z();
315  DK[0] = fRunningRes->GetKVer().X();
316  DK[1] = fRunningRes->GetKVer().Y();
317  DK[2] = fRunningRes->GetKVer().Z();
318 
319  fUtil->FromSCToSD(PC, RC, H, CH, DJ, DK, IERR, SPU, PD, RD);
320 
321  TVector3 positionsd = fUtil->FromMARSToSDCoord(TVector3(fFinalH->GetX(), fFinalH->GetY(), fFinalH->GetZ()), fRunningRes->GetOrigin(), fRunningRes->GetIVer(), fRunningRes->GetJVer(), fRunningRes->GetKVer());
322  Double_t u = positionsd.X(); // CHECK
323  Double_t v = positionsd.Y(); // CHECK
324  Double_t w = positionsd.Z(); // CHECK
325 
326  double spu;
327  fFinal->SetTrackPar(v, w, PD[1], PD[2], CH*PD[0], RD, fRunningRes->GetOrigin(), fRunningRes->GetIVer(), fRunningRes->GetJVer(),
328  fRunningRes->GetKVer(),spu);
329  // -----------------------------------------------------
330 
331  }
332 
333 
334  fTrackArray->Delete();
335  fHitArray->Delete();
336  fPointArray->Delete();
337 
338 }
339 
340 Bool_t FairGeaneTrKalStt::ProcessHit(PndSttTrack *pTrack, Int_t i, FairTrackParP *fRunningStart, FairTrackParP *fRunningRes, TString fb)
341 {
342  if(!fRunningStart) return kFALSE;
343  Int_t iHit = pTrack->GetHitIndex(i);
344  PndSttHit *currenthit = (PndSttHit*) fHitArray->At(iHit);
345  if(!currenthit) { cout << "NON HIT" << endl; return kFALSE;}
346  if(currenthit->GetXint() == -999 || currenthit->GetYint() == -999) { cout << "NON HIT INT" << endl; return kFALSE; }
347  Int_t refindex = currenthit->GetRefIndex();
348  // get point
349  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
350  // if(TMath::Sqrt(iPoint->GetXInLocal()*iPoint->GetXInLocal()+iPoint->GetYInLocal()*iPoint->GetYInLocal()) < 0.4998){cout << "OUT OF TUBE " << TMath::Sqrt(iPoint->GetXInLocal()*iPoint->GetXInLocal()+iPoint->GetYInLocal()*iPoint->GetYInLocal()) << endl; return kFALSE;}
351 
352  // check if points are too close
353  // TVector3 p1 = TVector3(fRunningStart->GetX(), fRunningStart->GetY(), 0.);
354  // TVector3 p2 = TVector3(currenthit->GetX(), currenthit->GetY(), 0.);
355  // TVector3 p3 = p1 - p2;
356  // if(p3.Mag() < 1.5) return kFALSE;
357 
358 
359 // cout << "START POINT " << fRunningStart->GetX()<<" "<< fRunningStart->GetY() << " " << fRunningStart->GetZ() << endl;
360 // cout << "START POINT ERR " << fRunningStart->GetDX()<<" " << fRunningStart->GetDY()<<" "<<fRunningStart->GetDZ() << endl;
361 // cout <<"START MOM " << fRunningStart->GetPx()<<" "<< fRunningStart->GetPy()<<" " << fRunningStart->GetPz() << endl;
362 // cout << "START MOM ERR " << fRunningStart->GetDPx()<<" "<<fRunningStart->GetDPy()<<" "<<fRunningStart->GetDPz() << endl;
363  // cout << "START POINT " << fRunningStart->GetV()<<" "<< fRunningStart->GetW() << endl;
364  // cout << "START ERR " << fRunningStart->GetDV()<<" "<< fRunningStart->GetDW() << endl;
365 
366  Bool_t rc_pro = Propagation(currenthit, fRunningStart, fRunningRes, fb);
367 
368  if(rc_pro == kFALSE) {
369  // cout << "RESET" << endl;
370  // fRunningRes->Reset();
371  return kFALSE;
372  }
373 
374  // cout << "EXTRAP POINT " << fRunningRes->GetX()<<" "<< fRunningRes->GetY() << " " << fRunningRes->GetZ() << endl;
375  // cout << "EXTRAP POINT ERR " << fRunningRes->GetDX()<<" " << fRunningRes->GetDY()<<" "<<fRunningRes->GetDZ() << endl;
376  // cout << "EXTRAP MOM " << fRunningRes->GetPx()<<" "<< fRunningRes->GetPy()<<" " << fRunningRes->GetPz() << endl;
377  // cout << "EXTRAP MOM ERR " << fRunningRes->GetDPx()<<" "<<fRunningRes->GetDPy()<<" "<<fRunningRes->GetDPz() << endl;
378  // cout << "EXTRAP POINT " << fRunningRes->GetV()<<" "<< fRunningRes->GetW() << endl;
379  // cout << "EXTRAP ERR " << fRunningRes->GetDV()<<" "<< fRunningRes->GetDW() << endl;
380 
381 
382  // kalman call
383  Bool_t rc_kal = Kalman(currenthit, fRunningRes, fRunningStart);
384  // cout << "KAL " << rc_kal << endl;
385  if(rc_kal == kFALSE) return kFALSE;
386  //......................................
387  // per controllare che il segno delle componenti sia corretto (da cancellare)
388  // if(premomx * fRunningStart->GetPx() < 0) cout << "AAAAAAAAIUTOOOOOOOOO X" << endl;
389  // if(premomy * fRunningStart->GetPy() < 0) cout << "AAAAAAAAIUTOOOOOOOOO Y" << endl;
390  // if(premomz * fRunningStart->GetPz() < 0) cout << "AAAAAAAAIUTOOOOOOOOO Z" << endl;
391  return kTRUE;
392 
393 }
394 
395 Bool_t FairGeaneTrKalStt::Propagation(PndSttHit *currenthit, FairTrackParP *fRunningStart, FairTrackParP *fRunningRes, TString fb)
396 {
397  // cout << "................." << endl;
398  // cout << "O " << fRunningStart->GetOrigin().X()<<" "<< fRunningStart->GetOrigin().Y() << " " << fRunningStart->GetOrigin().Z() << endl;
399  // cout << "I " << fRunningStart->GetIVer().X()<<" "<< fRunningStart->GetIVer().Y() << " " << fRunningStart->GetIVer().Z() << endl;
400  // cout << "J " << fRunningStart->GetJVer().X()<<" "<< fRunningStart->GetJVer().Y() << " " << fRunningStart->GetJVer().Z() << endl;
401  // cout << "K " << fRunningStart->GetKVer().X()<<" "<< fRunningStart->GetKVer().Y() << " " << fRunningStart->GetKVer().Z() << endl;
402  // cout << "R1START POINT " << fRunningStart->GetX()<<" "<< fRunningStart->GetY() << " " << fRunningStart->GetZ() << endl;
403  // cout << "R1START POINT ERR " << fRunningStart->GetDX()<<" " << fRunningStart->GetDY()<<" "<<fRunningStart->GetDZ() << endl;
404  // cout << "R1START MOM " << fRunningStart->GetPx()<<" "<< fRunningStart->GetPy()<<" " << fRunningStart->GetPz() << endl;
405  // cout << "R1START MOM ERR " << fRunningStart->GetDPx()<<" "<<fRunningStart->GetDPy()<<" "<<fRunningStart->GetDPz() << endl;
406 
407  // per controllare che il segno delle componenti sia corretto (da cancellare)
408  // double premomx = fRunningStart->GetPx();
409  // double premomy = fRunningStart->GetPy();
410  // double premomz = fRunningStart->GetPz();
411 
412  // ----- propagation: uso la propagate to closest ------
413  TVector3 v0 = TVector3(0., 0., 0.);
414  fPro->SetPoint(v0);
415 
416 
417  TVector3 center = TVector3(currenthit->GetX(), currenthit->GetY(), currenthit->GetZ());
418  TVector3 wiredir = currenthit->GetWireDirection();
419  wiredir *= 75.;
420  TVector3 second = (TVector3)(center + wiredir);
421  TVector3 first = (TVector3)(center - wiredir);
422 
423  fPro->SetWire(first, second);
424  if(fb == "FWD") fPro->PropagateToVirtualPlaneAtPCA(2); // 1 if point; 2 if wire
425  else if(fb == "BKW") fPro->BackTrackToVirtualPlaneAtPCA(2); // 1 if point; 2 if wire
426 
427  Bool_t rc = kFALSE;
428  rc = fPro->Propagate(fRunningStart, fRunningRes, PDGCode);
429  // cout << "PROPAGATION RC " << rc << " " << currenthit->GetX() << " " << currenthit->GetY() << " " << currenthit->GetZ() << endl;
430 
431  // if propagation failed
432  // if(fabs(fPro->GetPCAOnWire().X()) > 50. || fabs(fPro->GetPCAOnWire().Y()) > 50. || fabs(fPro->GetPCAOnWire().Z()) > 75.) rc = kFALSE;
433  if(fabs(fRunningRes->GetX()) > 50. || fabs(fRunningRes->GetY()) > 50. || fabs(fRunningRes->GetZ()) > 50) rc = kFALSE;
434 
435 
436 
437  return rc;
438 
439 }
440 Bool_t FairGeaneTrKalStt::Kalman(PndSttHit *currenthit, FairTrackParP *fRunningRes, FairTrackParP *fRunningStart)
441 {
442  // KALMAN
443  // punto estrapolato
444  TMatrixT<double> extrap(5,1);
445  extrap[0][0] = fRunningRes->GetQp()/fRunningRes->GetQ();
446  extrap[1][0] = fRunningRes->GetTV();
447  extrap[2][0] = fRunningRes->GetTW();
448  extrap[3][0] = fRunningRes->GetV();
449  extrap[4][0] = fRunningRes->GetW();
450  // cout << "extrap " << extrap[0][0] << " " << extrap[1][0] << " " << extrap[2][0] << " " << extrap[3][0] << " " << extrap[4][0] << endl;
451 
452  // punto estrapolato matrice errori
453  TMatrixT<double> sigmaex(5,5);
454  Double_t SIGMAEX[15], SIGMAEX55[5][5];
455  fRunningRes->GetCov(SIGMAEX);
456  fUtil->FromVec15ToMat25(SIGMAEX, SIGMAEX55);
457  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) sigmaex[i][j] = SIGMAEX55[i][j];
458  TMatrixT<double> weightex(5,5);
459  weightex = sigmaex;
460  Double_t det1 = weightex.Determinant();
461  if(det1 != 0) weightex.Invert();
462  else {
463  cout << "det1 = 0!!!!!" << endl;
464  return kFALSE;
465  }
466 
467  // punto misurato
468  Double_t rdrift = currenthit->GetIsochrone();
469  TMatrixT<double> meas(5,1);
470  meas[0][0] = 0.;
471  meas[1][0] = 0.;
472  meas[2][0] = 0.;
473  meas[3][0] = rdrift;
474  // meas[4][0] = currenthit->GetZint();
475  if(currenthit->GetWireDirection() == TVector3(0.,0.,1.)) meas[4][0] = 0.;
476  else if(currenthit->GetZint() != -999) meas[4][0] = currenthit->GetZint();
477  else meas[4][0] = 0.;
478 
479 
480  // punto misurato matrice pesi
481  TMatrixT<double> weightmi(5,5);
482  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) weightmi[i][j] = 0.;
483  weightmi[3][3] = 1./(0.0150 * 0.0150); // CHECK
484  if(meas[4][0] == 0) weightmi[4][4] = 1./(1.5 * 1.5); // CHECK
485  // cout << "MEAS POINT " << meas[3][0] <<" "<< meas[4][0] << endl;
486  // cout << "MEAS ERR " << sqrt(1./weightmi[3][3]);
487  // if(weightmi[4][4]!=0) cout << " " << sqrt(1./weightmi[4][4]) << endl;
488  // else cout << " 0" << endl;
489 
490  // cout << "meas " << meas[0][0] << " " << meas[1][0] << " " << meas[2][0] << " " << meas[3][0] << " " << meas[4][0] << endl;
491 
492  // somma pesiex + pesimi --> sigmakal = errori di kalman
493  TMatrixT<double> weightsum(5,5);
494  weightsum = weightex + weightmi;
495  // inverto
496  TMatrixT<double> sigmakal(5,5);
497  sigmakal = weightsum;
498  Double_t det2 = sigmakal.Determinant();
499  if(det2 != 0) sigmakal.Invert();
500  else {
501  cout << "det2 = 0!!!!!" << endl;
502  return kFALSE;
503  }
504 
505  // somma pesiex * extr + pesimi * meas
506  // pesiex * extr
507  TMatrixT<double> wextrap(weightex,TMatrixT<double>::kMult,extrap);
508  // pesimi * meas
509  TMatrixT<double> wmeas(weightmi,TMatrixT<double>::kMult,meas);
510  // sommo
511  TMatrixT<double> sum(5,1);
512  // for(int i = 0; i < 5; i++) wextrap[i][0] = fabs(wextrap[i][0]);
513  // for(int i = 0; i < 5; i++) wmeas[i][0] = fabs(wmeas[i][0]);
514 
515 
516 
517  sum = wextrap + wmeas;
518 
519  // cout << "wextrap " << wextrap[0][0] << " " << wextrap[1][0] << " " << wextrap[2][0] << " " << wextrap[3][0] << " " << wextrap[4][0] << endl;
520  // cout << "wmeas " << wmeas[0][0] << " " << wmeas[1][0] << " " << wmeas[2][0] << " " << wmeas[3][0] << " " << wmeas[4][0] << endl;
521 
522 
523 
524 
525 
526  // kalman point
527  TMatrixT<double> kalman(sigmakal,TMatrixT<double>::kMult,sum);
528 
529  // cout << "kalman " << kalman[0][0] << " " << kalman[1][0] << " " << kalman[2][0] << " " << kalman[3][0] << " " << kalman[4][0] << endl;
530 
531  Double_t SIGMAKAL[5][5], SIGMAKAL15[15];
532  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) SIGMAKAL[i][j] = sigmakal[i][j];
533  fUtil->FromMat25ToVec15(SIGMAKAL, SIGMAKAL15);
534 
535  fRunningStart->SetTrackPar(kalman[3][0], kalman[4][0],
536  kalman[1][0], kalman[2][0],
537  kalman[0][0] * fRunningRes->GetQ(), SIGMAKAL15, fRunningRes->GetOrigin(),
538  fRunningRes->GetIVer(), fRunningRes->GetJVer(), fRunningRes->GetKVer(),
539  fRunningRes->GetSPU());
540 
541  // cout << "................." << endl;
542  // cout << "O " << fRunningStart->GetOrigin().X()<<" "<< fRunningStart->GetOrigin().Y() << " " << fRunningStart->GetOrigin().Z() << endl;
543  // cout << "I " << fRunningStart->GetIVer().X()<<" "<< fRunningStart->GetIVer().Y() << " " << fRunningStart->GetIVer().Z() << endl;
544  // cout << "J " << fRunningStart->GetJVer().X()<<" "<< fRunningStart->GetJVer().Y() << " " << fRunningStart->GetJVer().Z() << endl;
545  // cout << "K " << fRunningStart->GetKVer().X()<<" "<< fRunningStart->GetKVer().Y() << " " << fRunningStart->GetKVer().Z() << endl;
546 
547  // cout << "KALMAN POINT " << fRunningStart->GetX()<<" "<< fRunningStart->GetY() << " " << fRunningStart->GetZ() << endl;
548  // cout << "KALMAN POINT ERR " << fRunningStart->GetDX()<<" " << fRunningStart->GetDY()<<" "<<fRunningStart->GetDZ() << endl;
549  // cout << "KALMAN MOM " << fRunningStart->GetPx()<<" "<< fRunningStart->GetPy()<<" " << fRunningStart->GetPz() << endl;
550  // cout << "KALMAN MOM ERR " << fRunningStart->GetDPx()<<" "<<fRunningStart->GetDPy()<<" "<<fRunningStart->GetDPz() << endl;
551  // cout << "KALMAN POINT " << fRunningStart->GetV()<<" "<< fRunningStart->GetW() << endl;
552  // cout << "KALMAN ERR " << fRunningStart->GetDV()<<" "<< fRunningStart->GetDW() << endl;
553 
554 
555  // cout << "KALMANMOMENTUM " << TMath::Sqrt(fRunningStart->GetPx() * fRunningStart->GetPx() + fRunningStart->GetPy() * fRunningStart->GetPy() + fRunningStart->GetPz() * fRunningStart->GetPz()) << endl;
556 
557 
558  return kTRUE;
559 
560 
561 }
562 
563 //-------------------
564 Bool_t FairGeaneTrKalStt::RetrieveVertex(PndSttTrack *pTrack)
565 {
566 
567  PndSttHelixTrackFitter *fFitter;
568  TVector3 *vertex = new TVector3(0.,0.,0.); // CHECK
569  TVector3 *recovtx = fFitter->PCAToPoint(pTrack, vertex);
570  TVector3 *momeatvtx;
571  momeatvtx = fFitter->MomentumAtPoint(pTrack, recovtx);
572 
573  // Double_t pterr = pTrack->TransverseMomentumError();
574  // Double_t plerr = pTrack->LongitudinalMomentumError();
575 
576  Double_t phi = pTrack->GetParamLast()->GetY();
577  Double_t d = pTrack->GetParamLast()->GetX();
578  Double_t cosp = TMath::Cos(phi);
579  Double_t sinp = TMath::Sin(phi);
580 // Double_t errd = pTrack->GetD0Error();
581 // Double_t errphi = pTrack->GetPhi0Error();
582  // Double_t errxp = sqrt(cosp * cosp * errd * errd + d * d * sinp * sinp * errphi * errphi);
583  // Double_t erryp = sqrt(sinp * sinp * errd * errd + d * d * cosp * cosp * errphi * errphi);
584  // Double_t errzp = pTrack->GetZ0Error();
585 
586  // ERRORI SUL VERTICE
587  // errors * 10;
588  // errxp *= 10.; erryp *= 10.; errzp *= 10.;
589  // pterr *= 10.; plerr *= 10.;
590  StartPos = *recovtx;
591  Double_t pterr, plerr;
592  Double_t errxp, erryp, errzp;
593  errxp = 10 * 0.02;
594  erryp = errxp;
595  errzp = 10 * 0.15;
596  StartPosErr = TVector3(errxp, erryp, errzp);
597  StartMom = *momeatvtx;
598 
599  Double_t momerr; // = 0.04 * 1.5;
600  momerr = 0.04 * 2.5;
601 
602  // pterr = 10 * 0.02;
603  // plerr = 10 * 0.02;
604  pterr = momerr;
605  plerr = momerr;
606 
607  StartMomErr = TVector3(pterr, pterr, plerr);
608  // ..................................................
609 
610  //cout << "vertex pos " << StartPos.X() << " " << StartPos.Y() << " " << StartPos.Z() << endl;
611  //cout << "vertex pos err " << StartPosErr.X() << " " << StartPosErr.Y() << " " << StartPosErr.Z() << endl;
612  //cout << "vertex mom " << StartMom.X() << " " << StartMom.Y() << " " << StartMom.Z() << endl;
613  //cout << "vertex mom err " << StartMomErr.X() << " " << StartMomErr.Y() << " " << StartMomErr.Z() << endl;
614 
615  return kTRUE;
616 
617 }
618 
619 Bool_t FairGeaneTrKalStt::BackToVertex(FairTrackParP *fRunningRes, FairTrackParP *fRunningKal)
620 {
621 
622  // cout << "back to vertex " << fRunningRes->GetQp() << " " << fRunningRes->GetTV() << " " << fRunningRes->GetTW() << " " << fRunningRes->GetV() << " " << fRunningRes->GetW() << endl;
623 
624 
625  if(fRunningRes->GetQp() == 0) return kFALSE;
626  // BACK TRACK TO VERTEX
627  TVector3 v0 = TVector3(0., 0., 0.);
628  fPro->SetPoint(v0);
629  fPro->SetWire(v0,v0);
630  fPro->BackTrackToVertex();
631  Bool_t rc_vtx = kFALSE;
632  TVector3 lastPos = TVector3(fRunningRes->GetX(),
633  fRunningRes->GetY(),
634  fRunningRes->GetZ());
635  TVector3 lastMom = TVector3(fRunningRes->GetPx(),
636  fRunningRes->GetPy(),
637  fRunningRes->GetPz());
638  TVector3 lastPosErr = TVector3(fRunningRes->GetDX(),
639  fRunningRes->GetDY(),
640  fRunningRes->GetDZ());
641  TVector3 lastMomErr = TVector3(fRunningRes->GetDPx(),
642  fRunningRes->GetDPy(),
643  fRunningRes->GetDPz());
644 
645  FairTrackParH *fKalH1 = new FairTrackParH(lastPos, lastMom,
646  lastPosErr, lastMomErr,
647  fRunningRes->GetQ());
648 
649 
650  // cout << "fKalH1 " << fKalH1->GetQp() << " " << fKalH1->GetX() << " " << fKalH1->GetY() << " " << fKalH1->GetZ() << endl;
651 
652  FairTrackParH *fKalH2 = new FairTrackParH();
653  rc_vtx = fPro->Propagate(fKalH1, fKalH2, PDGCode);
654  // cout << "fKalH2 " << fKalH2->GetQp() << " " << fKalH2->GetX() << " " << fKalH2->GetY() << " " << fKalH2->GetZ() << endl;
655 
656  if(rc_vtx == kFALSE) return kFALSE;
657  Double_t PC[3], RC[15], H[3], CH, DJ[3], DK[3], SPU, PD[3], RD[15];
658  Int_t IERR;
659  PC[0] = fKalH2->GetQp()/fKalH2->GetQ(); PC[1] = fKalH2->GetLambda(); PC[2] = fKalH2->GetPhi();
660  fKalH2->GetCov(RC);
661  H[0] = 0.; H[1] = 0.; H[2] = 20.;
662  CH = fKalH2->GetQ();
663  DJ[0] = 0.; DJ[1] = 1.; DJ[2] = 0;
664  DK[0] = 0.; DK[1] = 0.; DK[2] = 1;
665  fUtil->FromSCToSD(PC,RC, H,CH, DJ, DK, IERR, SPU,PD,RD);
666 
667  fRunningKal->SetTrackPar(fKalH2->GetX(), fKalH2->GetY(), fKalH2->GetZ(),
668  fKalH2->GetPx(), fKalH2->GetPy(), fKalH2->GetPz(),
669  fKalH2->GetQ(), RD,
670  TVector3(0.,0.,0.), TVector3(1.,0.,0.), TVector3(0.,1.,0.), TVector3(0.,0.,1.));
671  // cout << "fRunningKal " << fRunningKal->GetQp() << " " << fRunningKal->GetTV() << " " << fRunningKal->GetTW() << " " << fRunningKal->GetV() << " " << fRunningKal->GetW() << endl;
672  // cout << "VTX POINT " << fRunningKal->GetX()<<" "<< fRunningKal->GetY() << " " << fRunningKal->GetZ() << endl;
673  // cout << "VTX POINT ERR " << fRunningKal->GetDX()<<" " << fRunningKal->GetDY()<<" "<<fRunningKal->GetDZ() << endl;
674  // cout << "VTX MOM " << fRunningKal->GetPx()<<" "<< fRunningKal->GetPy()<<" " << fRunningKal->GetPz() << endl;
675  // cout << "VTX MOM ERR " << fRunningKal->GetDPx()<<" "<<fRunningKal->GetDPy()<<" "<<fRunningKal->GetDPz() << endl;
676  // cout << "VTX MOMENTUM " << TMath::Sqrt(fRunningKal->GetPx() * fRunningKal->GetPx() + fRunningKal->GetPy() * fRunningKal->GetPy() + fRunningKal->GetPz() * fRunningKal->GetPz()) << endl;
677 }
678 
679 Bool_t FairGeaneTrKalStt::BackToVertex2(FairTrackParP *fRunningRes, FairTrackParP *fRunningKal)
680 {
681 
682  if(!fRunningRes) return kFALSE;
683  // cout << "back to vertex " << fRunningRes->GetQp() << " " << fRunningRes->GetTV() << " " << fRunningRes->GetTW() << " " << fRunningRes->GetV() << " " << fRunningRes->GetW() << endl;
684 
685  // ----- propagation: uso la propagate to closest ------
686  TVector3 v0 = TVector3(0., 0., 0.);
687  fPro->SetPoint(v0);
688  fPro->SetWire(v0,v0);
689 
690  fPro->BackTrackToVirtualPlaneAtPCA(1); // 1 if point; 2 if wire
691 
692  Bool_t rc = kFALSE;
693  rc = fPro->Propagate(fRunningRes, fRunningKal, PDGCode);
694 
695  //cout << "fRunningKal " << fRunningKal->GetQp() << " " << fRunningKal->GetTV() << " " << fRunningKal->GetTW() << " " << fRunningKal->GetV() << " " << fRunningKal->GetW() << endl;
696  // cout << "VTX POINT " << fRunningKal->GetX()<<" "<< fRunningKal->GetY() << " " << fRunningKal->GetZ() << endl;
697  //cout << "VTX POINT ERR " << fRunningKal->GetDX()<<" " << fRunningKal->GetDY()<<" "<<fRunningKal->GetDZ() << endl;
698  // cout << "VTX MOM " << fRunningKal->GetPx()<<" "<< fRunningKal->GetPy()<<" " << fRunningKal->GetPz() << endl;
699  //cout << "VTX MOM ERR " << fRunningKal->GetDPx()<<" "<<fRunningKal->GetDPy()<<" "<<fRunningKal->GetDPz() << endl;
700  // cout << "VTX MOMENTUM " << TMath::Sqrt(fRunningKal->GetPx() * fRunningKal->GetPx() + fRunningKal->GetPy() * fRunningKal->GetPy() + fRunningKal->GetPz() * fRunningKal->GetPz()) << endl;
701 }
702 
703 
704 ClassImp(FairGeaneTrKalStt)
705 
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
i
Int_t i
Definition: run_full.C:25
FairGeaneTrKalStt::fPro
FairGeanePro * fPro
Definition: CbmGeaneTrKalStt.h:61
FairGeaneTrKalStt::Exec
virtual void Exec(Option_t *opt)
Definition: CbmGeaneTrKalStt.cxx:110
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
FairGeaneTrKalStt::ProcessHit
Bool_t ProcessHit(PndSttTrack *pTrack, Int_t k, FairTrackParP *fRunningStart, FairTrackParP *fRunningRes, TString fb)
Definition: CbmGeaneTrKalStt.cxx:340
kalman
KalmanTask * kalman
Definition: eta_c/run_reco_tpc_evt.C:80
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
v
__m128 v
Definition: P4_F32vec4.h:4
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
FairGeaneTrKalStt::fTrackArray
TClonesArray * fTrackArray
Definition: CbmGeaneTrKalStt.h:49
fParticle
Int_t * fParticle
Definition: run_full.C:24
PndSttPoint::GetYtot
Double_t GetYtot() const
Definition: PndSttPoint.h:55
FairGeaneTrKalStt::fUtil
FairGeaneUtil * fUtil
Definition: CbmGeaneTrKalStt.h:62
PndSttPoint::GetZtot
Double_t GetZtot() const
Definition: PndSttPoint.h:56
FairGeaneTrKalStt::Init
virtual InitStatus Init()
Definition: CbmGeaneTrKalStt.cxx:47
FairGeaneTrKalStt::BackToVertex
Bool_t BackToVertex(FairTrackParP *fRunningRes, FairTrackParP *fRes)
Definition: CbmGeaneTrKalStt.cxx:619
FairGeaneTrKalStt::fTrackParFinal
TClonesArray * fTrackParFinal
Definition: CbmGeaneTrKalStt.h:56
PndSttHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
FairGeaneTrKalStt::fTrackParGeane
TClonesArray * fTrackParGeane
Definition: CbmGeaneTrKalStt.h:55
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
FairGeaneTrKalStt::RetrieveVertex
Bool_t RetrieveVertex(PndSttTrack *pTrack)
Definition: CbmGeaneTrKalStt.cxx:564
phi
double phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
FairGeaneTrKalStt::BackToVertex2
Bool_t BackToVertex2(FairTrackParP *fRunningRes, FairTrackParP *fRes)
Definition: CbmGeaneTrKalStt.cxx:679
FairGeaneTrKalStt::Propagation
Bool_t Propagation(PndSttHit *currenthit, FairTrackParP *fRunningStart, FairTrackParP *fRunningRes, TString fb)
Definition: CbmGeaneTrKalStt.cxx:395
FairGeaneTrKalStt::fPointArray
TClonesArray * fPointArray
Definition: CbmGeaneTrKalStt.h:48
ClassImp
ClassImp(PndAnaContFact)
fb
TFile * fb
Definition: modifyRootFscGeometryFile.C:270
FairGeaneTrKalStt::Kalman
Bool_t Kalman(PndSttHit *currenthit, FairTrackParP *fRunningRes, FairTrackParP *fRunningStart)
Definition: CbmGeaneTrKalStt.cxx:440
PndSttPoint::GetXtot
Double_t GetXtot() const
Definition: PndSttPoint.h:54
eventnum3
Int_t eventnum3
Definition: CbmGeaneTrKalStt.cxx:33
PndSttTrack::GetFlag
Int_t GetFlag() const
Definition: PndSttTrack.h:51
FairGeaneTrKalStt::fHitArray
TClonesArray * fHitArray
Definition: CbmGeaneTrKalStt.h:47
FairGeaneTrKalStt::fTrackParIni
TClonesArray * fTrackParIni
Definition: CbmGeaneTrKalStt.h:54
Generated on Wed Apr 3 2019 08:43:49 for FairRoot/PandaRoot by   doxygen 1.8.5