FairRoot/PandaRoot
PndSttTrackFinderIdeal.cxx
Go to the documentation of this file.
1 
2 #include "PndSttTrackFinderIdeal.h"
3 
4 #include "PndSttHit.h"
5 #include "PndSttPoint.h"
6 #include "PndSttHelixHit.h"
7 #include "PndSttHoughDefines.h"
8 #include "PndDetectorList.h"
9 #include "PndTrackCand.h"
10 #include "PndTrack.h"
11 #include "PndSttTube.h"
12 
13 #include "FairMCPoint.h"
14 #include "FairTrackParP.h"
15 #include "FairRootManager.h"
16 
17 // ROOT includes
18 #include "TClonesArray.h"
19 #include "TCanvas.h"
20 #include "TH2F.h"
21 #include "TArc.h"
22 #include "TDatabasePDG.h"
23 #include "TRandom.h"
24 
25 // C++ includes
26 #include <iostream>
27 #include <map>
28 #include <cmath>
29 
30 using std::cout;
31 using std::cin;
32 using std::endl;
33 using std::map;
34 
35 // ----- Default constructor -------------------------------------------
36 PndSttTrackFinderIdeal::PndSttTrackFinderIdeal()
37 {
38  rootoutput = kFALSE;
39 
40  fMCTrackArray = NULL;
41  fVerbose = 1;
42  fHelixHitProduction = kFALSE;
43 }
44 // -------------------------------------------------------------------------
45 
46 
47 
48 // ----- Standard constructor ------------------------------------------
49 PndSttTrackFinderIdeal::PndSttTrackFinderIdeal(Int_t verbose)
50 {
51  fMCTrackArray = NULL;
52  fVerbose = verbose;
53  fHelixHitProduction = kFALSE;
54 
55  if (verbose >= 4) rootoutput = kTRUE;
56  else rootoutput = kFALSE; // stt1
57 
58 }
59 // -------------------------------------------------------------------------
60 
61 
62 
63 // ----- Destructor ----------------------------------------------------
64 PndSttTrackFinderIdeal::~PndSttTrackFinderIdeal()
65 {
66 }
67 // -------------------------------------------------------------------------
68 
69 
70 
71 // ----- Public method Init --------------------------------------------
72 void PndSttTrackFinderIdeal::Init()
73 {
74  // Get and check FairRootManager
75  FairRootManager* ioman = FairRootManager::Instance();
76 
77  if (!ioman)
78  {
79  cout << "-E- PndSttTrackFinderIdeal::Init: "
80  << "RootManager not instantised!" << endl;
81  return;
82  }
83 
84  // Get MCTrack array
85  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
86  if ( ! fMCTrackArray)
87  {
88  cout << "-E- PndSttTrackFinderIdeal::Init: No MCTrack array!"
89  << endl;
90  return;
91  }
92 }
93 
94 
95 Int_t PndSttTrackFinderIdeal::DoFind( TClonesArray* trackCandArray, TClonesArray *trackArray, TClonesArray* helixHitArray)
96 {
97  // Check pointers
98  if ( !fMCTrackArray )
99  {
100  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
101  << "MCTrack array missing! " << endl;
102  return -1;
103  }
104 
105  if (fHitCollectionList.GetEntries() == 0)
106  {
107  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
108  << "No hit arrays present, call AddHitCollection() first (at least once)! " << endl;
109  return -1;
110  }
111 
112  if (fPointCollectionList.GetEntries() == 0)
113  {
114  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
115  << "No point arrays present, call AddHitCollection() first (at least once)! " << endl;
116  return -1;
117  }
118 
119  if ( !trackCandArray )
120  {
121  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
122  << "Track cand array missing! " << endl;
123  return -1;
124  }
125 
126  if ( !trackArray )
127  {
128  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
129  << "Track array missing! " << endl;
130  return -1;
131  }
132 
133  TCanvas
134  *finderCanvas;
135 
136  if (rootoutput)
137  {
138  TH2F
139  myRange("range", "range", 100, -50., 50., 100, -50., 50);
140 
141  finderCanvas = new TCanvas("findercanvas", "findercanvas", 800, 600);
142  myRange.DrawCopy();
143  plotAllStraws();
144  }
145 
146 
147  // Initialise control counters
148  Int_t nNoMCTrack = 0;
149  Int_t nNoTrack = 0;
150  Int_t nNoSttPoint = 0;
151  Int_t nNoSttHit = 0;
152 
153  // Create pointers to hit and SttPoint
154  PndSttHit* pMhit = NULL;
155  FairMCPoint* pMCpt = NULL;
156  PndMCTrack* pMCtr = NULL;
157  PndTrackCand* pTrckCand = NULL;
158  PndTrack* pTrck = NULL;
159 
160  // Number of STT hits
161  Int_t
162  nHits = 0;
163 
164  for (Int_t hitListCounter = 0; hitListCounter < fHitCollectionList.GetEntries(); hitListCounter++)
165  {
166  nHits += ((TClonesArray *)fHitCollectionList.At(hitListCounter))->GetEntriesFast();
167  }
168 
169  // Declare some variables outside the loops
170  Int_t ptIndex = 0; // MCPoint index
171  Int_t mcTrackIndex = 0; // MCTrack index
172  Int_t trackIndex = 0; // STTTrack index
173 
174  // Create STL map from MCtrack index to number of valid SttHits
175  map<Int_t, map<Double_t, Int_t> >
176  hitMap;
177 
178  // Loop over hits
179  for (Int_t iHit = 0; iHit < nHits; iHit++)
180  {
181  pMhit = GetHitFromCollections(iHit);
182 
183  if ( ! pMhit )
184  continue;
185 
186  // tubeID CHECK added
187  Int_t tubeID = pMhit->GetTubeID();
188  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
189 
190  ptIndex = pMhit->GetRefIndex();
191 
192  if (rootoutput)
193  {
194  TArc
195  myArc;
196 
197  myArc.SetFillStyle(0);
198  myArc.SetLineColor(1);
199 
200  if (pMhit->GetIsochrone() == 0)
201  myArc.DrawArc(tube->GetPosition().X(), tube->GetPosition().Y(), 1.0);
202  else
203  myArc.DrawArc(tube->GetPosition().X(), tube->GetPosition().Y(), pMhit->GetIsochrone());
204  }
205 
206  if (ptIndex < 0)
207  continue; // fake or background hit
208 
209  pMCpt = GetPointFromCollections(iHit);
210 
211  if ( ! pMCpt )
212  continue;
213 
214  mcTrackIndex = pMCpt->GetTrackID();
215 
216  Double_t
217  wireX = tube->GetPosition().X(),
218  wireY = tube->GetPosition().Y();
219  (hitMap[mcTrackIndex])[wireX * wireX + wireY * wireY]++;
220  }
221 
222  // Create STL map from MCTrack index to SttTrack index
223  map<Int_t, Int_t>
224  correlationMap,
225  trackMap;
226 
227  // Create STTTracks for reconstructable MCTracks
228  Int_t nMCacc = 0; // accepted MCTracks (more than 3 points)
229  Int_t nTracks = 0; // reconstructable MCTracks
230  Int_t nMCTracks = fMCTrackArray->GetEntriesFast();
231 
232  for (Int_t iMCTrack = 0; iMCTrack < nMCTracks; iMCTrack++)
233  {
234  pMCtr = (PndMCTrack*) fMCTrackArray->At(iMCTrack);
235  if ( ! pMCtr )
236  continue;
237 
238  if (hitMap[iMCTrack].size() < 3)
239  continue;
240 
241  // temporary hack, fix this in monte carlo ....
242  if (TDatabasePDG::Instance()->GetParticle(pMCtr->GetPdgCode()) == NULL)
243  continue;
244 
245  if (!(fabs(TDatabasePDG::Instance()->GetParticle(pMCtr->GetPdgCode())->Charge() / 3.0) > 0.))
246  continue;
247  nMCacc++;
248 
249  new((*trackCandArray)[nTracks]) PndTrackCand();
250 
251  if (fVerbose >= 2) cout << "-I- PndSttTrackFinderIdeal: STTTrack "
252  << nTracks << " created from MCTrack "
253  << iMCTrack << " (" << pMCtr->GetNPoints(kSTT)
254  << " STTPoints)" << endl;
255 
256  correlationMap[nTracks] = iMCTrack;
257  trackMap[iMCTrack] = nTracks++;
258 
259 
260 
261  }
262 
263  // Loop over hits. Get corresponding MCPoint and MCTrack index
264  for (Int_t iHit = 0; iHit < nHits; iHit++)
265  {
266  pMhit = GetHitFromCollections(iHit);
267  if ( ! pMhit )
268  {
269  cout << "-E- PndSttTrackFinderIdeal::DoFind: Empty slot "
270  << "in HitArray at position " << iHit << endl;
271  nNoSttHit++;
272  continue;
273  }
274 
275  // tubeID CHECK added
276  Int_t tubeID = pMhit->GetTubeID();
277  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
278 
279  if (pMhit->GetIsochrone() == 0.)
280  {
281  if (fVerbose == 3) cout << "take the center: " << tube->GetPosition().X() << " " << tube->GetPosition().Y() << " " << tube->GetPosition().Z() << endl;
282  }
283 
284  ptIndex = pMhit->GetRefIndex();
285  if (ptIndex < 0)
286  continue; // fake or background hit
287 
288  pMCpt = GetPointFromCollections(iHit);
289 
290  if ( ! pMCpt )
291  {
292  nNoSttPoint++;
293  continue;
294  }
295 
296  mcTrackIndex = pMCpt->GetTrackID();
297 
298  if (mcTrackIndex < 0 || mcTrackIndex > nMCTracks)
299  {
300  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
301  << "MCTrack index out of range. " << mcTrackIndex << " "
302  << nMCTracks << endl;
303  nNoMCTrack++;
304  continue;
305  }
306 
307  if (trackMap.find(mcTrackIndex) == trackMap.end())
308  continue;
309 
310  trackIndex = trackMap[mcTrackIndex];
311  pTrckCand = (PndTrackCand*) trackCandArray->At(trackIndex);
312  if ( ! pTrckCand )
313  {
314  cout << "-E- PndSttTrackFinderIdeal::DoFind: "
315  << "No SttTrack pointer. " << iHit << " " << ptIndex
316  << " " << mcTrackIndex << " " << trackIndex << endl;
317  nNoTrack++;
318  continue;
319  }
320 
321  // ---------------------------------------------------------------------
322 
323 
324  TVector3 MCmom;
325  pMCpt->Momentum(MCmom);
326 
327  // tubeID CHECK added
328  tubeID = pMhit->GetTubeID();
329  tube = (PndSttTube*) fTubeArray->At(tubeID);
330 
331  Double_t wireRad = tube->GetPosition().Perp();
332 
333  if(MCmom.Mag() < 0.3) {
334  // for low momentum particles, hits added to pTrck by R could not be
335  // in the correct order due to spiralizing tracks
336  // cout << "LOW MOMENTUM --> AddHit by HitID " << endl;
337  if(pTrckCand->GetNHits() > 25) continue;
338  // CHECK: test iHit and how to organize sorting (here...)
339  // pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("STTHit"), iHit, iHit);
340  pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("STTHit"), iHit, iHit);
341 
342  }
343  else {
344  // cout << "HIGH MOMENTUM --> AddHit by R " << endl;
345  // CHECK: test iHit and how to organize sorting (... and here)
346  // pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("STTHit"), iHit, wireRad);
347  pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("STTHit"), iHit, wireRad);
348 
349  }
350 
351  // ---------------------------------------------------------------------
352  if(fHelixHitProduction) {
353 
354  // HelixHit Production after PR
355  TClonesArray& clref = *helixHitArray;
356  Int_t size = clref.GetEntriesFast();
357 
358  // CHECK remember to SET the errors on position!!!!!!!
359  TVector3 pos;
360  pMCpt->Position(pos);
361  TVector3 posErr(0, 0, 0);
362  PndSttHelixHit *helixhit = new(clref[size]) PndSttHelixHit(pMhit->GetDetectorID(), pMhit->GetTubeID(), iHit, pMhit->GetRefIndex(),
363  pos, posErr,
364  pMhit->GetIsochrone(), pMhit->GetIsochroneError(),
365  0.0);
366  // dedx
367  // calculate it from MC
368  double InOut[6];
369  memset(InOut, 0, sizeof(InOut));
370  InOut[0] = ((PndSttPoint *) pMCpt)->GetXInLocal();
371  InOut[1] = ((PndSttPoint *) pMCpt)->GetYInLocal();
372  InOut[2] = ((PndSttPoint *) pMCpt)->GetZInLocal();
373  InOut[3] = ((PndSttPoint *) pMCpt)->GetXOutLocal();
374  InOut[4] = ((PndSttPoint *) pMCpt)->GetYOutLocal();
375  InOut[5] = ((PndSttPoint *) pMCpt)->GetZOutLocal();
376  TVector3 diff3(InOut[0] - InOut[3], InOut[1] - InOut[4], InOut[2] - InOut[5]);
377  double distance = diff3.Mag(); //
378  Double_t dedx = 999;
379  Double_t eloss = ((PndSttPoint *) pMCpt)->GetEnergyLoss();
380  if (distance != 0) dedx = eloss/(distance); // in GeV/cm (I guess) CHECK
381  helixhit->SetdEdx(dedx);
382 
383  }
384  // ---------------------------------------------------------------------
385 
386  if (rootoutput)
387  {
388  TArc
389  myArc;
390 
391  myArc.SetFillStyle(0);
392  myArc.SetLineColor(trackIndex + 2);
393 
394  if (pMhit->GetIsochrone() == 0)
395  myArc.DrawArc(pMhit->GetX(), pMhit->GetY(), 1.0);
396  else
397  myArc.DrawArc(pMhit->GetX(), pMhit->GetY(), pMhit->GetIsochrone());
398  }
399 
400  if (fVerbose > 2) cout << "Hit " << iHit << " from STTPoint "
401  << ptIndex << " (MCTrack "
402  << mcTrackIndex << ") added to STTTrack "
403  << trackIndex << endl;
404  }
405 
406  for (Int_t trackTeller = 0; trackTeller < nTracks; trackTeller++)
407  {
408  // loop over
409  pTrckCand = (PndTrackCand*) trackCandArray->At(trackTeller);
410 
411  if ( pTrckCand != NULL)
412  {
413  pTrckCand->setMcTrackId(correlationMap[trackTeller]);
414 
415  Double_t
416  dSeed,
417  rSeed,
418  phiSeed,
419  zSeed,
420  tanLamSeed;
421 
422  GetTrack(dSeed, phiSeed, rSeed, zSeed, tanLamSeed, correlationMap[trackTeller]);
423 
424  Double_t
425  xSeed = (dSeed + rSeed) * cos(phiSeed),
426  ySeed = (dSeed + rSeed) * sin(phiSeed);
427 
428  //Double_t //[R.K. 01/2017] unused variable
429  //xSeed_old = xSeed, //[R.K. 01/2017] unused variable?
430  //ySeed_old = ySeed, //[R.K. 01/2017] unused variable?
431  //rSeed_old = rSeed; //[R.K. 01/2017] unused variable?
432 
433  // ZoomTrack(dSeed, phiSeed, rSeed, pTrck); // CHECK!!!
434 
435  xSeed = (dSeed + rSeed) * cos(phiSeed);
436  ySeed = (dSeed + rSeed) * sin(phiSeed);
437 
438  if (rootoutput)
439  {
440  TArc
441  myArc;
442 
443  myArc.SetFillStyle(0);
444  //myArc.SetLineColor(trackTeller + 2);
445  //myArc.DrawArc(xSeed_old, ySeed_old, rSeed_old);
446 
447  myArc.SetLineColor(trackTeller + 2);
448  myArc.DrawArc(xSeed, ySeed, rSeed);
449  }
450 
451  // ----------------------------------
452  // check: seeds directly from MC
453  Double_t vxSeed = dSeed * TMath::Cos(phiSeed);
454  Double_t vySeed = dSeed * TMath::Sin(phiSeed);
455  Double_t vzSeed = zSeed;
456  TVector3 posSeed(vxSeed, vySeed, vzSeed); // starting point
457 
458  PndMCTrack
459  *mcTrack2 = (PndMCTrack*) fMCTrackArray->At(correlationMap[trackTeller]);
460 
461  Int_t pdg = mcTrack2->GetPdgCode();
462  TDatabasePDG *fdbPDG = TDatabasePDG::Instance();
463  TParticlePDG *fParticle = fdbPDG->GetParticle(pdg);
464  Double_t chargeSeed = fParticle->Charge()/3.;
465 
466  TVector3 dirSeed(mcTrack2->GetMomentum().X(),
467  mcTrack2->GetMomentum().Y(),
468  mcTrack2->GetMomentum().Z()); // momentum direction in starting point
469 
470  Double_t qop = chargeSeed/dirSeed.Mag(); // q over p
471  dirSeed.SetMag(1.);
472 
473  // set up PndTrack ***********************************************************************
474  TVector3 momSeed(mcTrack2->GetMomentum().X(),
475  mcTrack2->GetMomentum().Y(),
476  mcTrack2->GetMomentum().Z()); // momentum direction in starting point
477 
478  // get 1st and last hits
479  Int_t hitcounter = pTrckCand->GetNHits();
480  PndTrackCandHit candhit = pTrckCand->GetSortedHit(0);
481  Int_t iHit = candhit.GetHitId();
482  PndSttPoint *firstpnt = (PndSttPoint*) GetPointFromCollections(iHit);
483  if(!firstpnt) { cout << "PndSttTrackFinderIdeal::DoFind ERROR: 1st pnt " << iHit << endl; continue; }
484  TVector3 MCfirstPos(firstpnt->GetX(), firstpnt->GetY(), firstpnt->GetZ());
485  TVector3 MCfirstMom(firstpnt->GetPx(), firstpnt->GetPy(), firstpnt->GetPz());
486  PndSttTube *firsttube = (PndSttTube*) fTubeArray->At(firstpnt->GetTubeID());
487  TVector3 difirst = MCfirstMom; difirst.SetMag(1.);
488  TVector3 djfirst = firsttube->GetWireDirection(); djfirst.SetMag(1.);
489  TVector3 dkfirst = difirst.Cross(djfirst); dkfirst.SetMag(1.);
490 
491  candhit = pTrckCand->GetSortedHit(hitcounter - 1);
492  iHit = candhit.GetHitId();
493  PndSttPoint *lastpnt = (PndSttPoint*) GetPointFromCollections(iHit);
494  if(!lastpnt) { cout << "PndSttTrackFinderIdeal::DoFind ERROR last pnt " << iHit << endl; continue; }
495  TVector3 MClastPos(lastpnt->GetX(), lastpnt->GetY(), lastpnt->GetZ());
496  TVector3 MClastMom(lastpnt->GetPx(), lastpnt->GetPy(), lastpnt->GetPz());
497  PndSttTube *lasttube = (PndSttTube*) fTubeArray->At(lastpnt->GetTubeID());
498  TVector3 dilast = MClastMom; dilast.SetMag(1.);
499  TVector3 djlast = lasttube->GetWireDirection(); djlast.SetMag(1.);
500  TVector3 dklast = dilast.Cross(djlast); dklast.SetMag(1.);
501 
502 
503  FairTrackParP first(MCfirstPos, MCfirstMom,
504  TVector3(0., 0., 0.), TVector3(0., 0., 0.), ((int) chargeSeed),
505  MCfirstPos, djfirst, dkfirst);
506 
507 
508  FairTrackParP last(MClastPos, MClastMom,
509  TVector3(0., 0., 0.), TVector3(0., 0., 0.), ((int) chargeSeed),
510  MClastPos, djlast, dklast);
511 
512  // pTrck = new((*trackArray)[trackTeller]) PndTrack(first, last, *pTrckCand, 0, -1., 0, 0, trackTeller, -1);
513  pTrck = new((*trackArray)[trackTeller]) PndTrack(first, last, *pTrckCand);
514  pTrck->SetRefIndex("STTTrackCand", trackTeller);
515 // pTrck->SetLink("STTTrackCand", trackTeller);
516 
517 
518  // *****************************************************************************************
519 
520 
521  // CHECK to be deleted: tests!
522  // position: MY posSeed, PREV dSeed
523  if(1==0) { // to pilot the print outs
524  if(fabs(dSeed - posSeed.Mag()) > 1e-3) {
525  cout << "pos" << endl;
526  cout << dSeed << " " << posSeed.Mag() << endl;
527  posSeed.Print() ;
528  }
529 
530  // direction: MY dirSeed, PREV perp to radius <-> (0, 0, 0)
531  TVector3 dtest(-TMath::Sin(phiSeed) * 0.006 * rSeed, TMath::Cos(phiSeed) * 0.006 * rSeed, mcTrack2->GetMomentum().Z());
532  dtest.SetMag(1.);
533  if(fabs(fabs(dtest.X()) - fabs(dirSeed.X())) > 1e-3) {
534  cout << "dir x " << endl;
535  dirSeed.Print();
536  dtest.Print();
537  }
538 
539  if(fabs(fabs(dtest.Y()) - fabs(dirSeed.Y())) > 1e-3) {
540  cout << "dir y " << endl;
541  dirSeed.Print();
542  dtest.Print();
543  }
544 
545  if(fabs(fabs(dtest.Z()) - fabs(dirSeed.Z())) > 1e-3) {
546  cout << "dir z " << endl;
547  dirSeed.Print();
548  dtest.Print();
549  }
550 
551  if(fabs(qop - (chargeSeed/((0.006 * rSeed) * (0.006 * rSeed) + mcTrack2->GetMomentum().Z() * mcTrack2->GetMomentum().Z()))) > 1e-3)
552  cout << "qop " << qop << " " << (chargeSeed/((0.006 * rSeed) * (0.006 * rSeed) + mcTrack2->GetMomentum().Z() * mcTrack2->GetMomentum().Z())) << " " << fabs(qop - (chargeSeed/((0.006 * rSeed) * (0.006 * rSeed) + mcTrack2->GetMomentum().Z() * mcTrack2->GetMomentum().Z()))) << endl;
553 
554  }
555 
556 
557  }
558  }
559 
560 
561 
562  if (rootoutput)
563  {
564  finderCanvas->Update();
565  finderCanvas->Show();
566 
567  char
568  waitchar;
569 
570  cout << "press any key to continue." << endl;
571  cin >> waitchar;
572 
573  delete finderCanvas;
574  }
575 
576  if (fVerbose == 3)
577  {
578  cout << endl;
579  cout << "-------------------------------------------------------"
580  << endl;
581  cout << "-I- Ideal STT track finding -I-"
582  << endl;
583  cout << "Hits: " << nHits << endl;
584  cout << "MCTracks: total " << nMCTracks << ", accepted " << nMCacc
585  << ", reconstructable: " << nTracks << endl;
586 
587  cout << "SttHits not found : "
588  << nNoSttHit << endl;
589  cout << "SttPoints not found : "
590  << nNoSttPoint << endl;
591  cout << "MCTracks not found : "
592  << nNoMCTrack << endl;
593  cout << "SttTracks not found : "
594  << nNoTrack << endl;
595  cout << "-------------------------------------------------------"
596  << endl;
597  }
598  else if(fVerbose) cout << "-I- PndSttTrackFinderIdeal: all " << nMCTracks
599  << ", acc. " << nMCacc << ", rec. " << nTracks << endl;
600 
601  return nTracks;
602 }
603 
604 void PndSttTrackFinderIdeal::GetTrackletCircular(Double_t firstX, Double_t firstY, Double_t firstR,
605  Double_t secondX, Double_t secondY, Double_t secondR,
606  Double_t thirdX, Double_t thirdY, Double_t thirdR,
607  Double_t *circleRadii, Double_t *circleCentersX,
608  Double_t *circleCentersY) const
609 {
610  Int_t
611  trackletCounter = 0;
612 
613  Double_t
614  small_limit = 0.0001;
615 
616  for (Int_t sign1 = -1; sign1 < 3; sign1 += 2)
617  {
618  for (Int_t sign2 = -1; sign2 < 3; sign2 += 2)
619  {
620  for (Int_t sign3 = -1; sign3 < 3; sign3 += 2)
621  {
622  // if three points are collinear, shift middle point by 50 microns
623  if ((firstX - secondX == 0) && (secondX - thirdX == 0))
624  {
625  secondX += 0.005;
626  secondY += 0.005;
627  }
628  else if (!((fabs(firstX - secondX) < small_limit) || (fabs(secondX - thirdX) < small_limit)))
629  {
630  if (fabs((firstY - secondY) / (firstX - secondX) - (secondY - thirdY) / (secondX - thirdX)) < small_limit)
631  {
632  secondX += 0.005;
633  secondY += 0.005;
634  }
635  }
636 
637  Double_t
638  a = -2. * (firstX - secondX),
639  b = -2. * (firstY - secondY),
640  c = 2. * (sign1 * firstR - sign2 * secondR),
641  d = ((firstX * firstX + firstY * firstY - firstR * firstR) -
642  (secondX * secondX + secondY * secondY - secondR * secondR)),
643  e = -2. * (firstX - thirdX),
644  f = -2. * (firstY - thirdY),
645  g = 2. * (sign1 * firstR - sign3 * thirdR),
646  h = ((firstX * firstX + firstY * firstY - firstR * firstR) -
647  (thirdX * thirdX + thirdY * thirdY - thirdR * thirdR));
648 
649  Double_t
650  A = -1. * (f * d - b * h) / (a * f - b * e),
651  B = -1. * (f * c - b * g) / (a * f - b * e),
652  C = (e * d - a * h) / (a * f - e * b),
653  D = (e * c - a * g) / (a * f - e * b),
654 
655  I = B * B + D * D - 1.,
656  II = 2. * A * B - 2. * B * firstX + 2. * C * D - 2. * D * firstY + sign1 * 2. * firstR,
657  III = A * A - 2. * A * firstX + firstX * firstX + C * C - 2. * C * firstY + firstY * firstY - firstR * firstR;
658 
659  if ((fabs(I) > small_limit) && ((II * II - 4. * I * III) > 0.))
660  {
661  Double_t
662  r = (-1. * II - sqrt(II * II - 4. * I * III)) / (2. * I),
663 
664  x = A + B * r,
665  y = C + D * r;
666 
667  circleRadii[trackletCounter] = sqrt(r * r);
668  circleCentersX[trackletCounter] = x;
669  circleCentersY[trackletCounter] = y;
670  trackletCounter++;
671  }
672  else
673  {
674  circleRadii[trackletCounter] = -1.;
675  circleCentersX[trackletCounter] = 0.;
676  circleCentersY[trackletCounter] = 0.;
677  trackletCounter++;
678  }
679  }
680  }
681  }
682 }
683 
684 // void PndSttTrackFinderIdeal::ZoomTrack(Double_t &dSeed, Double_t &phiSeed,
685 // Double_t &rSeed, PndSttTrack *track)
686 // {
687 // track->SortHits();
688 
689 // if (track->GetNofHits() < 3)
690 // {
691 // cout << "-W- PndSttTrackFinderHough::GetSeed2Circular: less than three hits found, not possible to fit!"
692 // << endl;
693 // return;
694 // }
695 // else
696 // {
697 // PndSttHoughAccumulatorNew
698 // *accumulator = new PndSttHoughAccumulatorNew(2, kTRUE,
699 // 150, dSeed - 5., dSeed + 5.,
700 // 150, phiSeed - 0.5, phiSeed + 0.5,
701 // 100, rSeed - 200., rSeed + 200.);
702 
703 // Int_t
704 // segmentLength1 = track->GetNofHits() / 3,
705 // segmentLength2 = track->GetNofHits() / 3,
706 // segmentLength3 = track->GetNofHits() - segmentLength2 - segmentLength1;
707 
708 // Int_t
709 // segmentStart1 = 0,
710 // segmentStart2 = segmentLength1,
711 // segmentStart3 = segmentLength1 + segmentLength2;
712 
713 // // select n random combinations
714 // for (Int_t teller = 0; teller < 100; teller++)
715 // {
716 // Int_t
717 // firstSample = gRandom->Integer(segmentLength1),
718 // secondSample = gRandom->Integer(segmentLength2),
719 // thirdSample = gRandom->Integer(segmentLength3);
720 
721 // Int_t
722 // i = segmentStart1 + firstSample,
723 // ii = segmentStart2 + secondSample,
724 // iii = segmentStart3 + thirdSample;
725 
726 // PndSttHit
727 // *iFirstHit = GetHitFromCollections(track->GetHitIndex(i)),
728 // *iSecondHit = GetHitFromCollections(track->GetHitIndex(ii)),
729 // *iThirdHit = GetHitFromCollections(track->GetHitIndex(iii));
730 
731 // if ((!iFirstHit) || (!iSecondHit) || (!iThirdHit))
732 // continue;
733 
734 // Double_t
735 // firstX = iFirstHit->GetX(),
736 // firstY = iFirstHit->GetY(),
737 // firstR = iFirstHit->GetIsochrone(),
738 // secondX = iSecondHit->GetX(),
739 // secondY = iSecondHit->GetY(),
740 // secondR = iSecondHit->GetIsochrone(),
741 // thirdX = iThirdHit->GetX(),
742 // thirdY = iThirdHit->GetY(),
743 // thirdR = iThirdHit->GetIsochrone();
744 
745 // Double_t
746 // circleRadii[8],
747 // circleCentersX[8],
748 // circleCentersY[8];
749 
750 // GetTrackletCircular(firstX, firstY, firstR,
751 // secondX, secondY, secondR,
752 // thirdX, thirdY, thirdR,
753 // circleRadii, circleCentersX, circleCentersY);
754 
755 // // loop over all tracklets between this and next
756 // for (Int_t trackletteller = 0; trackletteller < 8; trackletteller++)
757 // {
758 // if (circleRadii[trackletteller] > 0.)
759 // {
760 // Double_t
761 // phi = atan(circleCentersY[trackletteller] / circleCentersX[trackletteller]),
762 // dist = sqrt(circleCentersX[trackletteller] * circleCentersX[trackletteller] +
763 // circleCentersY[trackletteller] * circleCentersY[trackletteller])
764 // - circleRadii[trackletteller];
765 
766 // if (circleCentersX[trackletteller] < 0.)
767 // {
768 // if (circleCentersY[trackletteller] < 0.)
769 // phi -= dPi;
770 // else
771 // phi += dPi;
772 // }
773 
774 // accumulator->AddHit(dist, i,
775 // phi, ii,
776 // circleRadii[trackletteller], iii);
777 // }
778 // }
779 // }
780 
781 // accumulator->Cluster();
782 
783 // dSeed = accumulator->GetHighestPeakX();
784 // phiSeed = accumulator->GetHighestPeakY();
785 // rSeed = accumulator->GetHighestPeakZ();
786 
787 // delete accumulator;
788 // }
789 // }
790 
791 
792 
793 
794 void PndSttTrackFinderIdeal::GetTrack(Double_t &dSeed, Double_t &phiSeed, Double_t &rSeed,
795  Double_t &zSeed, Double_t &tanLamSeed, Int_t mcTrackNo)
796 {
797  PndMCTrack
798  *mcTrack = (PndMCTrack*) fMCTrackArray->At(mcTrackNo);
799 
800  // TODO: read field from container
801  rSeed = sqrt(mcTrack->GetMomentum().X() * mcTrack->GetMomentum().X() +
802  mcTrack->GetMomentum().Y() * mcTrack->GetMomentum().Y()) / 0.006;
803 
804  Double_t
805  phiStart = atan(mcTrack->GetMomentum().Y() / mcTrack->GetMomentum().X()),
806  xStart = mcTrack->GetStartVertex().X(),
807  yStart = mcTrack->GetStartVertex().Y();
808 
809  if (mcTrack->GetMomentum().X() < 0.)
810  {
811  if (mcTrack->GetMomentum().Y() < 0.)
812  phiStart -= dPi;
813  else
814  phiStart += dPi;
815  }
816 
817  Double_t
818  sign = 1.;
819 
820  if (TDatabasePDG::Instance()->GetParticle(mcTrack->GetPdgCode()) != NULL)
821  {
822  if (TDatabasePDG::Instance()->GetParticle(mcTrack->GetPdgCode())->Charge() < 0.)
823  {
824  sign = -1.;
825  }
826  }
827 
828  Double_t
829  xCircleCenter = xStart + sign * (rSeed * sin(phiStart)),
830  yCircleCenter = yStart - sign * (rSeed * cos(phiStart));
831 
832  dSeed = sqrt(xCircleCenter * xCircleCenter + yCircleCenter * yCircleCenter) - rSeed;
833  phiSeed = atan(yCircleCenter / xCircleCenter);
834 
835 
836  if (xCircleCenter < 0.)
837  {
838  if (yCircleCenter < 0.)
839  phiSeed -= dPi;
840  else
841  phiSeed += dPi;
842  }
843 
844  // z - tan(lam)
845  tanLamSeed = mcTrack->GetMomentum().Z() / sqrt(mcTrack->GetMomentum().X() * mcTrack->GetMomentum().X() +
846  mcTrack->GetMomentum().Y() * mcTrack->GetMomentum().Y());
847 
848  zSeed = mcTrack->GetStartVertex().Z();
849 
850 
851 }
852 
853 // -------------------------------------------------------------------------
854 Bool_t PndSttTrackFinderIdeal::putStraw(Double_t xpos, Double_t ypos, Double_t radius)
855 {
856  Double_t
857  pipeDiam = 0.42,
858  tubeOuterDiam = 1.032,
859  CoverThickness = 0.1,
860  outerRadius = 42.,
861  innerRadius = 15.;
862 
863  if ((sqrt(xpos * xpos + ypos * ypos) < outerRadius - CoverThickness - (tubeOuterDiam / 2.)) &&
864  (sqrt(xpos * xpos + ypos * ypos) > innerRadius + CoverThickness + (tubeOuterDiam / 2.)) &&
865  (sqrt(ypos * ypos) > (pipeDiam + tubeOuterDiam / 2.)))
866  {
867  TArc
868  myArc;
869 
870  myArc.SetFillStyle(0);
871  myArc.SetLineColor(17);
872  myArc.DrawArc(xpos, ypos, radius);
873  }
874  else
875  {
876  return false;
877  }
878 
879  return true;
880 }
881 
882 void PndSttTrackFinderIdeal::plotAllStraws()
883 {
884  Double_t
885  tubeOuterDiam = 1.032,
886  sttCenterX = 0.,
887  sttCenterY = 0.;
888 
889  Int_t
890  ringmax = -1;
891 
892  Bool_t
893  started = kFALSE,
894  goOn = kTRUE;
895 
896  Int_t
897  ringteller = 18;
898 
899  while (goOn)
900  {
901  goOn = kFALSE;
902 
903  Double_t
904  sqrt3 = sqrt(3.),
905  radius = tubeOuterDiam / 2.;
906 
907  Double_t
908  zpos = radius;
909 
910  // place first
911  Double_t
912  xpos = sttCenterX + ((ringteller) * 2 * radius),
913  ypos = sttCenterY;
914 
915  for(Int_t i = 0; i < ringteller; i++)
916  {
917  xpos -= radius;
918  ypos += sqrt3 * radius;
919  if (putStraw(xpos, ypos, zpos))
920  goOn = kTRUE;
921  }
922  for(Int_t i = 0; i < ringteller; i++)
923  {
924  xpos -= 2 * radius;
925  if (putStraw(xpos, ypos, zpos))
926  goOn = kTRUE;
927  }
928  for(Int_t i = 0; i < ringteller; i++)
929  {
930  xpos -= radius;
931  ypos -= sqrt3 * radius;
932  if (putStraw(xpos, ypos, zpos))
933  goOn = kTRUE;
934  }
935  for(Int_t i = 0; i < ringteller; i++)
936  {
937  xpos += radius;
938  ypos -= sqrt3 * radius;
939  if (putStraw(xpos, ypos, zpos))
940  goOn = kTRUE;
941  }
942  for(Int_t i = 0; i < ringteller; i++)
943  {
944  xpos += 2 * radius;
945  if (putStraw(xpos, ypos, zpos))
946  goOn = kTRUE;
947  }
948  for(Int_t i = 0; i < ringteller; i++)
949  {
950  xpos += radius;
951  ypos += sqrt3 * radius;
952  if (putStraw(xpos, ypos, zpos))
953  goOn = kTRUE;
954  }
955 
956  if (goOn)
957  started = kTRUE;
958 
959  if (!started)
960  goOn = kTRUE;
961 
962  ringteller++;
963  if ((ringmax > 0) && (ringteller == ringmax))
964  goOn = kFALSE;
965  }
966 }
967 
968 
969 PndSttHit* PndSttTrackFinderIdeal::GetHitFromCollections(Int_t hitCounter)
970 {
971  PndSttHit
972  *retval = NULL;
973 
974  Int_t
975  relativeCounter = hitCounter;
976 
977  for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++)
978  {
979  Int_t
980  size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast();
981 
982  if (relativeCounter < size)
983  {
984  retval = (PndSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter);
985  break;
986  }
987  else
988  {
989  relativeCounter -= size;
990  }
991  }
992  return retval;
993 }
994 
995 FairMCPoint* PndSttTrackFinderIdeal::GetPointFromCollections(Int_t hitCounter)
996 {
997  FairMCPoint
998  *retval = NULL;
999 
1000  Int_t
1001  relativeCounter = hitCounter;
1002 
1003  for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++)
1004  {
1005  Int_t
1006  size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast();
1007 
1008  if (relativeCounter < size)
1009  {
1010  Int_t
1011  tmpHit = ((PndSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter))->GetRefIndex();
1012 
1013  retval = (FairMCPoint*) ((TClonesArray *)fPointCollectionList.At(collectionCounter))->At(tmpHit);
1014 
1015  break;
1016  }
1017  else
1018  {
1019  relativeCounter -= size;
1020  }
1021  }
1022  return retval;
1023 }
1024 
1025 
1026 
1027 
1028 ClassImp(PndSttTrackFinderIdeal)
1029 
1030 
1031 
1032 
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndTrack::SetRefIndex
void SetRefIndex(TString branch, Int_t i)
Definition: PndTrack.h:41
trackArray
TClonesArray * trackArray
Definition: NHitsPerEvent.C:14
r
double r
Definition: RiemannTest.C:14
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
PndSttPoint::GetTubeID
Int_t GetTubeID()
Definition: PndSttPoint.h:77
i
Int_t i
Definition: run_full.C:25
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
PndSttTrackFinderIdeal::GetHitFromCollections
PndSttHit * GetHitFromCollections(Int_t hitCounter)
Definition: PndSttTrackFinderIdeal.cxx:969
map
PndTransMap * map
Definition: sim_emc_apd.C:99
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndMCTrack::GetNPoints
Int_t GetNPoints(DetectorId detId) const
Definition: PndMCTrack.cxx:120
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
sttCenterX
#define sttCenterX
Definition: createSTT.C:46
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndTrackCand::GetSortedHit
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
PndTrackCand::setMcTrackId
void setMcTrackId(int i)
Definition: PndTrackCand.h:72
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
g
TFile * g
Definition: GammaSpectraAnalysis.C:22
eloss
double eloss
Definition: anaLmdSim.C:34
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
GetEntriesFast
cout<< "POINTs for new FwEndCap == "<< tsim-> GetEntriesFast()
Definition: anal_point_fwendcap.C:90
PndSttTrackFinderIdeal::DoFind
virtual Int_t DoFind(TClonesArray *trackCandArray, TClonesArray *trackArray, TClonesArray *helixHitArray)
Definition: PndSttTrackFinderIdeal.cxx:95
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndSttTrackFinderIdeal::putStraw
Bool_t putStraw(Double_t xpos, Double_t ypos, Double_t radius)
Definition: PndSttTrackFinderIdeal.cxx:854
At
cout<< "blue = Monte Carlo "<< endl;cout<< "red = Helix Hit "<< endl;cout<< "green = Center Of Tubes "<< endl;for(Int_t k=0;k< track->GetEntriesFast();k++){PndSttTrack *stttrack=(PndSttTrack *) track-> At(k)
Definition: checkhelixhit.C:56
kSTT
Definition: production/tdrct/eta_c/track_check.C:4
C
int Pic_FED Eff_lEE C()
Definition: Pic_FEP-Eff_lE2.C:1
fParticle
Int_t * fParticle
Definition: run_full.C:24
a
Int_t a
Definition: anaLmdDigi.C:126
PndTrackCand::AddHit
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Definition: PndTrackCand.cxx:36
nHits
int nHits
Definition: RiemannTest.C:16
PndSttHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
zpos
const Double_t zpos
Definition: createRootFscGeometryFile.C:11
PndSttTube::GetPosition
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
pipeDiam
#define pipeDiam
Definition: createSTT.C:55
f
TFile * f
Definition: bump_analys.C:12
PndSttHit::GetTubeID
Int_t GetTubeID() const
Definition: PndSttHit.h:75
PndSttHit::GetIsochroneError
Double_t GetIsochroneError() const
Definition: PndSttHit.h:63
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
sttCenterY
#define sttCenterY
Definition: createSTT.C:47
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndTrackCand::GetNHits
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
dPi
#define dPi
Definition: PndSttHoughDefines.h:4
ClassImp
ClassImp(PndAnaContFact)
sign
int sign(T val)
Definition: PndCADef.h:48
y
Double_t y
Definition: createRootFscGeometryFile.C:159
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndSttTrackFinderIdeal::GetTrackletCircular
void GetTrackletCircular(Double_t firstX, Double_t firstY, Double_t firstR, Double_t secondX, Double_t secondY, Double_t secondR, Double_t thirdX, Double_t thirdY, Double_t thirdR, Double_t *circleRadii, Double_t *circleCentersX, Double_t *circleCentersY) const
Definition: PndSttTrackFinderIdeal.cxx:604
PndTrackCandHit::GetHitId
Int_t GetHitId() const
Definition: PndTrackCandHit.h:55
PndSttHelixHit::SetdEdx
void SetdEdx(Double_t dedx)
Definition: PndSttHelixHit.h:65
PndSttTrackFinderIdeal::GetTrack
void GetTrack(Double_t &dSeed, Double_t &phiSeed, Double_t &rSeed, Double_t &zSeed, Double_t &tanLamSeed, Int_t mcTrackNo)
Definition: PndSttTrackFinderIdeal.cxx:794
tubeOuterDiam
#define tubeOuterDiam
Definition: createSTT.C:38
PndSttTrackFinderIdeal::GetPointFromCollections
FairMCPoint * GetPointFromCollections(Int_t hitCounter)
Definition: PndSttTrackFinderIdeal.cxx:995
Generated on Sun Aug 18 2019 03:01:50 for FairRoot/PandaRoot by   doxygen 1.8.8