FairRoot/PandaRoot
PndMvdSttGemRiemannTrackFinder.cxx
Go to the documentation of this file.
1 #include "PndMvdSttGemRiemannTrackFinder.h"
2 #include <iostream>
3 #include <math.h>
4 #include "PndSdsHit.h"
5 #include "PndGemHit.h"
6 #include "PndSdsHit.h"
7 #include "PndSttHit.h"
8 #include "TString.h"
9 #include "TMath.h"
10 
11 using namespace std;
12 
13  ClassImp(PndMvdSttGemRiemannTrackFinder);
14 
15 PndMvdSttGemRiemannTrackFinder::PndMvdSttGemRiemannTrackFinder():PndRiemannTrackFinder(), fSttHits(0), fGemHits(0),fZClosePar(0.1), fLastLayerId(0)
16 {
17  fLayers.resize(22);
18  if (fUseZeroPos)
19  fLayers[0].push_back(0);
20  fGeoH = PndGeoHandling::Instance();
21  fSttHitsInSectors.resize(6);
22  fSttSkewedHitsInSectors.resize(6);
23  InitLayerMap();
24 }
25 
26 PndMvdSttGemRiemannTrackFinder::~PndMvdSttGemRiemannTrackFinder(){}
27 
28 void PndMvdSttGemRiemannTrackFinder::InitLayerMap()
29 {
30  InitLayerMapMvd();
31  InitLayerMapGem();
32 }
33 
34 void PndMvdSttGemRiemannTrackFinder::InitLayerMapMvd()
35 {
36  fLayerMap["PixeloBlo1"] = fLastLayerId + 1;
37  fLayerMap["PixeloSdko(Silicon)_1"] = fLastLayerId + 2; //naming after MVD2.2
38  fLayerMap["PixeloSdkoco(Silicon)_1"] = fLastLayerId + 2;
39  fLayerMap["PixeloSdko(Silicon)_2"] = fLastLayerId + 3; //naming after MVD2.2
40  fLayerMap["PixeloSdkoco(Silicon)_2"] = fLastLayerId + 3;
41  fLayerMap["PixeloSdko(Silicon)_3"] = fLastLayerId + 4; //naming after MVD2.2
42  fLayerMap["PixeloSdkoco(Silicon)_3"] = fLastLayerId + 4;
43  fLayerMap["PixeloSdko(Silicon)_4"] = fLastLayerId + 5; //naming after MVD2.2
44  fLayerMap["PixeloSdkoco(Silicon)_4"] = fLastLayerId + 5;
45  fLayerMap["PixeloBlo2"] = fLastLayerId + 6;
46  fLayerMap["PixeloLdkoio(Silicon)_1"] = fLastLayerId + 7; //naming after MVD2.2
47  fLayerMap["PixeloLdkoco(Silicon)_1"] = fLastLayerId + 7;
48  fLayerMap["PixeloLdkoiio(Silicon)_1"] = fLastLayerId + 8; //naming after MVD2.2
49  fLayerMap["PixeloLdkoco(Silicon)_2"] = fLastLayerId + 8;
50  fLayerMap["PixeloLdkoiiio(Silicon)_1"] = fLastLayerId + 9; //naming after MVD2.2
51  fLayerMap["PixeloLdkoco(Silicon)_3"] = fLastLayerId + 9;
52  fLayerMap["PixeloLdkoiiio(Silicon)_2"] = fLastLayerId + 10; //naming after MVD2.2
53  fLayerMap["PixeloLdkoco(Silicon)_4"] = fLastLayerId + 10;
54  fLayerMap["PixeloLdkoiio(Silicon)_2"] = fLastLayerId + 11; //naming after MVD2.2
55  fLayerMap["PixeloLdkoco(Silicon)_5"] = fLastLayerId + 11;
56  fLayerMap["PixeloLdkoio(Silicon)_2"] = fLastLayerId + 12; //naming after MVD2.2
57  fLayerMap["PixeloLdkoco(Silicon)_6"] = fLastLayerId + 12;
58  fLayerMap["PixeloLdkoiiio(Silicon)_3"] = fLastLayerId + 13; //naming after MVD2.2
59  fLayerMap["PixeloLdkoco(Silicon)_7"] = fLastLayerId + 13;
60  fLayerMap["PixeloLdkoiiio(Silicon)_4"] = fLastLayerId + 14; //naming after MVD2.2
61  fLayerMap["PixeloLdkoco(Silicon)_8"] = fLastLayerId + 14;
62 
63  fLayerMap["StripoBl3o(Silicon)"] = fLastLayerId + 15;
64  fLayerMap["Fwdo(Silicon)_1"] = fLastLayerId + 16;
65  fLayerMap["StripoBl4o(Silicon)"] = fLastLayerId + 17;
66  //fLayerMap["Fwdo(Silicon)_2")){Layer=10;flag=false;}
67 
68  fLayerMap["StripoLdkoTrapSoRingAoSilicon_1"] = fLastLayerId + 18;
69  fLayerMap["StripoLdkoTrapSoRingAoSilicon_2"] = fLastLayerId + 19;
70  fLayerMap["StripoLdkoTrapSoRingBoSilicon_1"] = fLastLayerId + 18;
71  fLayerMap["StripoLdkoTrapSoRingBoSilicon_2"] = fLastLayerId + 19;
72  fLayerMap["StripoLdko5-6oTrapSo(Silicon)_1"] = fLastLayerId + 20;
73 
74  fLayerMap["LambdaDisk_1"] = fLastLayerId + 21;
75  fLayerMap["LambdaDisk_2"] = fLastLayerId + 22;
76  fLastLayerId += 22;
77 }
78 
79 void PndMvdSttGemRiemannTrackFinder::InitLayerMapGem()
80 {
81  fLayerMap["Gem_1_1"] = fLastLayerId + 1;
82  fLayerMap["Gem_1_2"] = fLastLayerId + 2;
83 
84  fLayerMap["Gem_2_1"] = fLastLayerId + 3;
85  fLayerMap["Gem_2_2"] = fLastLayerId + 4;
86 
87  fLayerMap["Gem_3_1"] = fLastLayerId + 5;
88  fLayerMap["Gem_3_2"] = fLastLayerId + 6;
89 
90  fLastLayerId += 6;
91 
92 }
93 
94 
95 void PndMvdSttGemRiemannTrackFinder::AddHits(TClonesArray* hits, Int_t branchId)
96 {
97 // SetVerbose(3);
98  FairRootManager* man = FairRootManager::Instance();
99  // PndRiemannTrackFinder::AddHits(hits, branchId);
100  if (branchId == man->GetBranchId("STTHit")){
101  fSttHits = hits;
102  for (int i = 0; i < fSttHits->GetEntries(); i++){
103  PndSttHit* myHit = (PndSttHit*)(fSttHits->At(i));
104  myHit->SetEntryNr(FairLink(man->GetBranchId("SttHit"), i));
105  fHits.push_back(myHit);
106  myHit->SetPositionError(TVector3(0.5/TMath::Sqrt(12), 0.5/TMath::Sqrt(12), 200));
107  TVector3 sttPos;
108  myHit->Position(sttPos);
109  FairLink myID;
110 
111  if (myHit->GetEntryNr().GetIndex() < 0 || myHit->GetEntryNr().GetType() < 0){
112  myID = FairLink(branchId, i);
113  myHit->SetEntryNr(myID);
114  }
115  else
116  myID = myHit->GetEntryNr();
117  double phi = 0;
118 // if (sttPos.Phi() < 0)
119 // phi = sttPos.Phi() + TMath::Pi() * 2;
120 // else
121 // phi = sttPos.Phi();
122  if (fVerbose > 1) std::cout << "PndMvdSttGemRiemannTrackFinder::AddHits StrawMap Sector: " << myHit->GetEntryNr() << " " << myHit->GetTubeID() << ": Phi: " << sttPos.x() << "/" << sttPos.y() << " " << phi << " Sector: " << fStrawMap.GetSector(myHit->GetTubeID())<< std::endl;
123  fSttHitsInSectors[fStrawMap.GetSector(myHit->GetTubeID())].push_back(myHit);
124 
125  fMapHitToID[fHits.size()-1]=myID;
126  fMapIDtoHit[myID] = fHits.size()-1;
127  }
128  return;
129  }
130 
131  if(branchId == man->GetBranchId("STTCombinedSkewedHits")){
132  if (hits == 0){
133  std::cout << "PndMvdSttGemRiemannTrackFinder::AddHits: no SkewedHits Branch" << std::endl;
134  return;
135  }
136  std::cout << "PndMvdSttGemRiemannTrackFinder::AddHits: SkewedHits " << hits->GetEntriesFast() << std::endl;
137  fSkewedSttHits = hits;
138  for (int i = 0; i < fSkewedSttHits->GetEntries(); i++){
139  FairHit* myHit = (FairHit*)(fSkewedSttHits->At(i));
140  myHit->SetEntryNr(FairLink(man->GetBranchId("STTCombinedSkewedHits"), i));
141  fHits.push_back(myHit);
142  myHit->SetPositionError(TVector3(0.5/TMath::Sqrt(12), 0.5/TMath::Sqrt(12), 200));
143  TVector3 sttPos;
144  myHit->Position(sttPos);
145  FairLink myID;
146 
147  if (myHit->GetEntryNr().GetIndex() < 0 || myHit->GetEntryNr().GetType() < 0){
148  myID = FairLink(branchId, i);
149  myHit->SetEntryNr(myID);
150  }
151  else
152  myID = myHit->GetEntryNr();
153  double phi = 0;
154  if (sttPos.Phi() < 0)
155  phi = sttPos.Phi() + TMath::Pi() * 2;
156  else
157  phi = sttPos.Phi();
158 // if (fVerbose > 1)
159  std::cout << "PndMvdSttGemRiemannTrackFinder::AddHits StrawMap Sector: " << myHit->GetEntryNr() << ": Phi: " << sttPos.x() << "/" << sttPos.y() << " " << phi << " Sector: " << fStrawMap.FindPhiSector(phi)<< std::endl;
160  fSttSkewedHitsInSectors[fStrawMap.FindPhiSector(phi)].push_back(myHit);
161 
162  fMapHitToID[fHits.size()-1]=myID;
163  fMapIDtoHit[myID] = fHits.size()-1;
164  }
165  return;
166  }
167 
168  if (branchId == man->GetBranchId("GEMHit")){
169  fGemHits = hits;
170  if (fVerbose > 1) std::cout << "PndMvdSttGemRiemannTrackFinder::AddHits GEMHitsAdded: " << hits->GetEntries() << std::endl;
171  for (int i = 0; i < hits->GetEntries(); i++){
172  FairHit* myHit = (FairHit*)(hits->At(i));
173  fHits.push_back(myHit);
174  FairLink myID;
175  if (myHit->GetEntryNr().GetIndex() < 0){
176  myID = FairLink(branchId, i);
177  myHit->SetEntryNr(myID);
178  }
179  else
180  myID = myHit->GetEntryNr();
181  fMapHitToID[fHits.size()-1]=myID;
182  fMapIDtoHit[myID] = fHits.size()-1;
183  }
184  return;
185  }
186 
187  for (int i = 0; i < hits->GetEntries(); i++){
188  FairHit* myHit = (FairHit*)(hits->At(i));
189  fHits.push_back(myHit);
190  FairLink myID;
191  if (myHit->GetEntryNr().GetIndex() < 0){
192  myID = FairLink(branchId, i);
193  myHit->SetEntryNr(myID);
194  }
195  else
196  myID = myHit->GetEntryNr();
197  fMapHitToID[fHits.size()-1]=myID;
198  fMapIDtoHit[myID] = fHits.size()-1;
199 
200  int Layer = 0;
201  if (branchId == man->GetBranchId("MVDHitsStrip") || branchId == man->GetBranchId("MVDHitsPixel")){
202  Layer = GetLayerMvd(myHit);
203  } else if (branchId == man->GetBranchId("GEMHit")) {
204  Layer = GetLayerGem(myHit);
205  }
206 
207  fNLayers = fLayers.size();
208  while (fNLayers < Layer+1){
209  std::vector<int> dummy;
210  fLayers.push_back(dummy);
211  fNLayers = fLayers.size();
212  }
213 
214 
215  if (fVerbose > 1) std::cout << "fMapHitToId: " << fHits.size() -1 << " : " << myID << " "
216  << myHit->GetX() << "/" << myHit->GetY() << "/" << myHit->GetZ()
217  << " Layer: " << Layer << std::endl;
218 
219  if (Layer > 0)
220  fLayers[Layer].push_back(fHits.size()-1); //putting hit in layers array
221  }
222  fNLayers = fLayers.size();
223 }
224 
225 
226 int PndMvdSttGemRiemannTrackFinder::GetLayer(TString identifier)
227 {
228  std::map<TString, int>::iterator layerIter;
229  for (layerIter = fLayerMap.begin(); layerIter != fLayerMap.end(); layerIter++){
230  if(identifier.Contains(layerIter->first)){
231  return layerIter->second;
232  }
233  }
234  return 0;
235 }
236 
237 int PndMvdSttGemRiemannTrackFinder::GetLayerGem(FairHit* hit)
238 {
239  PndGemHit* gemHit = (PndGemHit*)hit;
240 
241  TString prefix("Gem_");
242  prefix+=(gemHit->GetStationNr());
243  prefix+=("_");
244  prefix+=(gemHit->GetSensorNr());
245 
246  return GetLayer(prefix);
247 
248 }
249 
250 int PndMvdSttGemRiemannTrackFinder::GetLayerMvd(FairHit* hit)
251 {
252  PndSdsHit* tempHit = (PndSdsHit*) (hit);
253  TString geoPath = fGeoH->GetPath(tempHit->GetSensorID());
254 
255  return GetLayer(geoPath);
256 
257 }
258 
259 void PndMvdSttGemRiemannTrackFinder::FindTracks()
260 {
261 
262 // SetVerbose(3);
263 
264  std::vector<std::set<Int_t> > Tracks = GetStartTracks(); //Get the possible track seeds
265 
266  std::set<int> tooClose;
267  std::set<int>::iterator iter;
268 
269  fTracks.clear();
270  fTrackCand.clear();
271  for (unsigned int trackId = 0; trackId < Tracks.size(); trackId++){ //Go through all track seeds and search for additional points
272 
273  if (Tracks[trackId].size() != 3)
274  if (fVerbose > 1) std::cout << "-E- PndMVDRiemannTrackFinder::FindTracks: Start points: " << Tracks[trackId].size()
275  << " in Track: " << trackId << std::endl;
276 
277  std::set<Int_t> StartTrack = Tracks[trackId];
278 
279  if (fVerbose > 1){
280  if (fVerbose > 1)
281  {
282  std::cout << "------------------------------------" << std::endl;
283 
284  iter = StartTrack.begin();
285  std::cout << "Start Plane from Points: " << fMapHitToID[*iter] << " "
286  << fMapHitToID[*(++iter)] << " " << fMapHitToID[*(++iter)] << std::endl;
287  iter = StartTrack.begin();
288  std::cout << "Start Plane from Points: " << *iter << " " << *(iter++) << " " << *(iter++) << std::endl;
289  }
290  }
291 
292  if (TrackExists(StartTrack) == true){
293  if (fVerbose > 1) std::cout << "Track exists already!" << std::endl;
294  continue;
295  }
296 
297  PndRiemannTrack actTrack = CreateRiemannTrack(StartTrack);
298  int startLayer=0;
299  int startHit=*(StartTrack.begin());
300 
301  bool flag=false;
302  for(int i=1;i<fNLayers;i++){
303  for(unsigned int j=0;j<fLayers[i].size();j++){
304  if(fLayers[i][j]==startHit){
305  startLayer=i;
306  flag=true;
307  break;
308  }
309  }
310  if (flag) break;
311  }
312 
313  int testHit;
314  for (int Layer=startLayer;Layer<fNLayers;Layer++){
315  //if ((fHits[startHit]->GetZ()<(-fZClosePar)) && (Layer==2 or Layer==3 or Layer==5 or Layer==6 or Layer==8 or Layer==10)) continue; /// < in case of backward tracks disk layers can't contain hits
316  for(unsigned int testHitInLayer=0; testHitInLayer<fLayers[Layer].size();testHitInLayer++){
317  testHit=fLayers[Layer][testHitInLayer];
318 
319  if (fVerbose > 2) if (fVerbose > 1) std::cout << "Layer: " << Layer << " hitInLayer: " << testHitInLayer << " hitID " << testHit << " ";
320  if (fVerbose > 1) std::cout << "Point " << fMapHitToID[testHit] << " ";
321 
322  if (CheckHitInTrack(StartTrack, testHit)) continue;
323  iter = StartTrack.begin();
324  if (CheckHitDistance(*(iter) , testHit)!=true) continue;
325  if (CheckHitDistance(*(++iter), testHit)!=true) continue;
326  if (CheckHitDistance(*(++iter), testHit)!=true) continue;
327  if (CheckZeroPassing(StartTrack, testHit)== true) continue;
328 
329 
330  //if ((fHits[startHit]->GetZ())*(fHits[testHit]->GetZ())<0 && fabs(fHits[startHit]->GetZ())>fZClosePar && fabs(fHits[testHit]->GetZ())>fZClosePar) continue; //check the same direction on z axis
331 
332  PndRiemannHit actHit(fHits[testHit], testHit);
333  if (CheckRiemannHit(&actTrack, &actHit, fHits[testHit]) != true) continue;
334 
335  StartTrack.insert(testHit);
336 
337  actTrack.addHit(actHit);
338  actTrack.refit(false);
339  actTrack.szFit(false);
340 
341  TVectorD orig = actTrack.orig();
342  if (fVerbose > 1) std::cout << "actHit added: " << testHit << " r: " << actTrack.r() << " orig: " << orig[0] << " " << orig[1] << std::endl;
343  tooClose=GetTooCloseHitsInLayer(Layer,testHit);
344  fHitsTooClose[trackId].insert(tooClose.begin(),tooClose.end());
345  //break;
346  }
347  }
348 
349  if (fVerbose > 1) {
350  std::cout << std::endl;
351  std::cout << "PndMvdSttGemRiemannTrackFinder::FindTracks ActTrack for Assignment: ";
352  for (iter = StartTrack.begin(); iter != StartTrack.end(); iter++){
353  int hitId = *iter;
354  std::cout << fMapHitToID[hitId] << " ";
355  }
356  std::cout << std::endl;
357  }
358 
359  if (fVerbose > 1) std::cout << "Track before STT: " << actTrack << std::endl;
360 
361  if (fSttHits != NULL){ //assign STTHits to Mvd track
362  if (fVerbose > 1 )
363  std::cout << "PndMvdSttGemRiemannTrackFinder::FindTracks AssignSttHits" << std::endl;
364  AssignSttHits(actTrack, StartTrack);
365  if (fVerbose > 1 )
366  std::cout << "ActTrack size with Stt: " << actTrack.getNumHits() << " " << StartTrack.size() << std::endl;
367  actTrack.refit(false);
368  actTrack.szFit(false);
369  actTrack.correctSttHits();
370  actTrack.refit(false);
371  actTrack.szFit(false);
372  actTrack.correctSttHits();
373  actTrack.refit(false);
374  actTrack.szFit(false);
375 
376 
377  }
378 
379  if (fSkewedSttHits != NULL){
380  if (fVerbose > 1 )
381  std::cout << "PndMvdSttGemRiemannTrackFinder::FindTracks AssignSkewedSttHits" << std::endl;
382  AssignSkewedSttHits(actTrack, StartTrack);
383  if (fVerbose > 1 )
384  std::cout << "ActTrack size with Stt: " << actTrack.getNumHits() << " " << StartTrack.size() << std::endl;
385  actTrack.refit(false);
386  actTrack.szFit(false);
387  }
388 
389  if (fVerbose > 1) std::cout << "Track after STT: " << actTrack << std::endl;
390 
391  if (fGemHits != NULL){
392  if (fVerbose > 1 )
393  if (fVerbose > 1) std::cout << "PndMvdSttGemRiemannTrackFinder::FindTracks AssignGemHits" << std::endl;
394  AssignGemHits(actTrack, StartTrack);
395  if (fVerbose > 1 )
396  if (fVerbose > 1) std::cout << "ActTrack size with Gem: " << actTrack.getNumHits() << " " << StartTrack.size() << std::endl;
397  actTrack.refit(false);
398  actTrack.szFit(false);
399 
400  }
401 
402  if (fVerbose > 1) std::cout << "Track after GEM: " << actTrack << std::endl;
403 
404 
405  if ((int)actTrack.getNumHits() > (fMinNumberOfHits-1)) //if you have a track match check hits which were to close
406  {
407  std::set<int> hits = fHitsTooClose[trackId];
408  if (fVerbose > 1) std::cout << "HitsTooClose Test for Track: " << actTrack << std::endl;
409  for (iter = hits.begin(); iter != hits.end(); iter++){
410  if (fVerbose > 1) std::cout << "Too Close Point " << *iter << ": " << fMapHitToID[*iter];
411  if (CheckHitInTrack(StartTrack, *iter)) continue;
412  PndRiemannHit actHit(fHits[*iter]);
413  if (CheckRiemannHit(&actTrack, &actHit, fHits[*iter])!= true) continue;
414  StartTrack.insert(*iter);
415  actTrack.addHit(actHit);
416 
417  }
418 
419 // actTrack.refit(false);
420 // actTrack.szFit(false);
421 
422  if (fVerbose > 1) std::cout << "Hits in Track: " << StartTrack.size() << std::endl;
423 
424  fTracks.push_back(actTrack);
425  fHitsInTracks.push_back(StartTrack);
426 
427  if (fVerbose > 1){
428  TVectorD orig = actTrack.orig();
429  if (fVerbose > 1)
430  iter = StartTrack.begin();
431  std::cout << "Track added! " << *iter << " " << *(++iter)
432  << " " << *(++iter) << " r: " << actTrack.r()
433  << " orig: " << orig[0] << " " << orig[1]
434  << " sz-m: " << actTrack.getSZm() << " sz-t: " << actTrack.getSZt()
435  << " dip: " << actTrack.dip()
436  << std::endl;
437  }
438  if (fVerbose > 1) {
439  for (iter = StartTrack.begin(); iter != StartTrack.end(); iter++)
440  {
441  std::cout << " " << fMapHitToID[*iter];
442  }
443  }
444  }
445  }
446 
447  // if (fVerbose > 1) std::cout << "In PndMvdRiemannTrackFinder are: " << fTracks.size() << " tracks" << std::endl;
448  for (unsigned int n = 0; n < fTracks.size(); n++){
449  PndTrackCand myTrackCand;
450  bool trackFound = false;
451  std::vector<PndRiemannHit> TrackHits = fTracks[n].getHits();
452 
453  for (unsigned int p = 0; p < TrackHits.size(); p++){
454  if (fVerbose > 1) std::cout << "HitID: " << TrackHits[p].hitID() << std::endl;
455  if (TrackHits[p].hitID() > -1){
456  myTrackCand.AddHit(fMapHitToID[TrackHits[p].hitID()], TrackHits[p].s());
457  }
458 
459  }
460  if (fVerbose > 1) std::cout << "NPoints TrackCand: " << myTrackCand.GetNHits() << std::endl;
461  if (fVerbose > 1) myTrackCand.Print();
462  FairMultiLinkedData gemHits = myTrackCand.GetLinksWithType(FairRootManager::Instance()->GetBranchId("GEMHit"));
463  FairMultiLinkedData sttHits = myTrackCand.GetLinksWithType(FairRootManager::Instance()->GetBranchId("STTHit"));
464  if (gemHits.GetNLinks() > 3){
465  fTrackCand.push_back(myTrackCand);
466  trackFound = true;
467  }
468  else if (sttHits.GetNLinks() > 3 && CheckBoarderHitsStt(&myTrackCand)) {
469  fTrackCand.push_back(myTrackCand);
470  trackFound = true;
471 // } else if (gemHits.GetNLinks() == 0 && sttHits.GetNLinks() == 0) {
472 // fTrackCand.push_back(myTrackCand);
473 // trackFound = true;
474  } else {
475  fTracks.erase(fTracks.begin()+n);
476  n--;
477  }
478  if (trackFound) {
479 // fTracks[n].refit(false);
480 // fTracks[n].szFit(false);
481  std::pair<double,double> CurvDip(1/fTracks[n].r(),fTracks[n].dip());
482  fCurvAndDipOfCand.push_back(CurvDip);
483 
484  if (fVerbose > 0){
485  std::cout << "PndMvdSttGemRiemannTrackFinder::FindTracks TrackFound: " << myTrackCand << std::endl;
486  }
487  if (fVerbose > 1) std::cout << " numHits: " << fTracks[n].getNumHits() << std::endl;
488 // if (fVerbose > 1)
489  {
490  TVectorD myOrig = fTracks[n].orig();
491  std::cout << "Track: " << fTracks[n] << std::endl;
492  std::cout << " curv: " << 1/fTracks[n].r() << "+/-" << fTracks[n].dR()/(fTracks[n].r() * fTracks[n].r())
493  << " dip: " << fTracks[n].dip() << "+/-" << fTracks[n].dDip()
494  << " orig: " << myOrig[0] << "+/-" << fTracks[n].dX()
495  << " " << myOrig[1] << "+/-" << fTracks[n].dY() << std::endl;
496  std::cout << "PlaneChiSquare: " << fTracks[n].calcChi2Plane() << std::endl;
497  }
498  }
499  }
500 
501  // if (fVerbose > 1) std::cout << "Hits before merging: " << fTrackCand.size() << std::endl;
502  // MergeTracks();
503  // fTrackCand.clear();
504  // fTrackCand = fMergedTrackCand;
505  // if (fVerbose > 0) {
506  // if (fVerbose > 1) std::cout << "Tracks after merging:" << fMergedTrackCand.size() << std::endl;
507  // for (unsigned int p = 0; p < fTrackCand.size(); p++){
508  // fTrackCand[p].Print();
509  // } ///todo Create RiemannTracks out of PndTrackCands
510  // }
511 
512 }
513 
514 std::vector< std::set<Int_t> > PndMvdSttGemRiemannTrackFinder::GetStartTracks()
515 {
516  std::set<Int_t> actCandidates;
517  std::vector<std::set<Int_t> > Tracks;
518  std::set<int> tooCloseFirst;
519  std::set<int> tooCloseSecond;
520  std::set<int> tooCloseThird;
521  if (fHits.size() > 3) {
522  int shift = 0;
523  if (fUseZeroPos)
524  shift = 1;
525  for (int FirstLayer = 1 - shift; FirstLayer < fNLayers - 2; FirstLayer++) {
526  for (int SecondLayer = FirstLayer + 1; SecondLayer < fNLayers - 1; SecondLayer++) {
527 // for (int ThirdLayer = SecondLayer + 1; ThirdLayer < fNLayers; ThirdLayer++) {
528  for (int ThirdLayer = fNLayers - 1; ThirdLayer > SecondLayer; ThirdLayer--) {
529 
530  for (unsigned int firstInLayer = 0; firstInLayer < fLayers[FirstLayer].size(); firstInLayer++) {
531  int first = fLayers[FirstLayer][firstInLayer];
532  // if ((fHits[first]->GetZ()<(-fZClosePar)) && (SecondLayer==2 or SecondLayer==3 or SecondLayer==5 or SecondLayer==6 or SecondLayer==8 or SecondLayer==10)) continue;
533  for (unsigned int secondInLayer = 0; secondInLayer < fLayers[SecondLayer].size(); secondInLayer++) {
534  int second = fLayers[SecondLayer][secondInLayer];
535  // if ((fHits[second]->GetZ()<(-fZClosePar)) && (ThirdLayer==2 or ThirdLayer==3 or ThirdLayer==5 or ThirdLayer==6 or ThirdLayer==8 or ThirdLayer==10)) continue;
536  // if ((fHits[first]->GetZ())*(fHits[second]->GetZ())<0 && fabs(fHits[first]->GetZ())>fZClosePar && fabs(fHits[second]->GetZ())>fZClosePar) {/*printf("Diff Sign of Points z1=%e z2=%e \n",fHits[first]->GetZ(),fHits[second]->GetZ()); */continue;} //my// check the same direction on z axis
537 
538  for (unsigned int thirdInLayer = 0; thirdInLayer < fLayers[ThirdLayer].size(); thirdInLayer++) {
539  int third = fLayers[ThirdLayer][thirdInLayer];
540  //if ((fHits[second]->GetZ())*(fHits[third]->GetZ())<0 && fabs(fHits[second]->GetZ())>fZClosePar && fabs(fHits[third]->GetZ())>fZClosePar) {/*printf("Diff Sign of Points z1=%e z2=%e \n",fHits[first]->GetZ(),fHits[second]->GetZ());*/ continue;} //my// check the same direction on z axis
541  if (fVerbose > 2)
542  std::cout
543  << "Checking Points for Start Triplet: "
544  << first << " " << second << " "
545  << third << std::endl;
546  if (fVerbose > 2)
547  std::cout
548  << "Checking Points for Start Triplet: "
549  << fMapHitToID[first] << " "
550  << fMapHitToID[second] << " "
551  << fMapHitToID[third] << std::endl;
552 
553  if (CheckHitDistance(first, third) != true) {
554  tooCloseFirst.insert(third);
555  continue;
556  }
557  if (CheckHitDistance(second, third) != true) {
558  tooCloseSecond.insert(third);
559  continue;
560  }
561  if (CheckHitInSameSensor(first, third) == true)
562  continue;
563  if (CheckHitInSameSensor(second, third) == true)
564  continue;
565 
566  actCandidates.clear();
567  actCandidates.insert(first);
568  actCandidates.insert(second);
569 
570  if (CheckZeroPassing(actCandidates, third)== true) continue;
571 
572  actCandidates.insert(third);
573 
574  PndRiemannTrack actTrack; // = new PndRiemannTrack();
575  PndRiemannHit hit1(fHits[first]);
576  PndRiemannHit hit2(fHits[second]);
577  PndRiemannHit hit3(fHits[third]);
578  actTrack.addHit(hit1);
579  actTrack.addHit(hit2);
580  actTrack.addHit(hit3);
581  actTrack.refit(false);
582  if (CheckSZ(actTrack) != true)
583  continue;
584  TVectorT<double> orig = actTrack.orig();
585 // if (fVerbose > 1)
586  {
587  std::set<Int_t>::iterator iter = actCandidates.begin();
588  std::cout << "Base plane from Points: "
589  << *(iter) << " " << *(iter++) << " "
590  << *(iter++) << std::endl;
591  std::cout << "TrackParameter: " << actTrack << std::endl;
592  }
593  Tracks.push_back(actCandidates);
594  tooCloseFirst = GetTooCloseHitsInLayer(
595  FirstLayer, first);
596  tooCloseSecond = GetTooCloseHitsInLayer(
597  SecondLayer, second);
598  tooCloseThird = GetTooCloseHitsInLayer(
599  ThirdLayer, third);
600  std::set<int> combHits = tooCloseFirst;
601  combHits.insert(tooCloseSecond.begin(), tooCloseSecond.end());
602  combHits.insert(tooCloseThird.begin(), tooCloseThird.end());
603  fHitsTooClose.push_back(combHits);
604 
605  }
606 
607  }
608 
609  }
610 
611  }
612  }
613  }
614  }
615 
616  if (fVerbose > 1) {
617  std::cout << "Start Tracks are: " << std::endl;
618  for (unsigned int i = 0; i < Tracks.size(); i++) {
619  std::set<int> aTrack = Tracks[i];
620  for (std::set<int>::iterator iter = aTrack.begin(); iter != aTrack.end(); iter++) {
621  std::cout << *iter << " ";
622  }
623  std::cout << std::endl;
624  for(std::set<int>::iterator iter = aTrack.begin(); iter != aTrack.end(); iter++) {
625  std::cout << fMapHitToID[*iter] << " ";
626  }
627  std::cout << std::endl;
628  }
629  std::cout << std::endl;
630  }
631  return Tracks;
632 }
633 
634 
635 
636 void PndMvdSttGemRiemannTrackFinder::AssignSttHits(PndRiemannTrack& actTrack, std::set<Int_t>& startTrack)
637 {
638  int startSector = GetStrawSector(actTrack);
639  int searchSector = startSector;
640  bool continueSearch = true;
641  int nextSectorDirection = 0;
642  int nHitsStt = -1;
643 
644  if (fVerbose > 2)
645  std::cout << "StartSector: " << searchSector << std::endl;
646 
647  while (continueSearch){
648  if (nHitsStt == 0){ //check if you have found any Stt hits in the selected startSector. If this is not the case take the left one. If there are still no Stt hits take the right one
649  if (startSector == searchSector){
650  searchSector = fStrawMap.GetRightSector(startSector);
651  }
652  else {
653  searchSector = fStrawMap.GetLeftSector(startSector);
654  }
655  }
656 
657  std::vector<PndSttHit*> sttHits = fSttHitsInSectors[searchSector];
658  continueSearch = false;
659  if (fVerbose > 2)
660  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: Actual searchSector: " << searchSector << std::endl;
661 
662  for (size_t i = 0; i < sttHits.size(); i++){
663 
664  if (nHitsStt < 0) nHitsStt = 0;
665 
666  PndSttHit* fairSttHit = sttHits[i];
667 
668  if (fStrawMap.IsSkewedStraw(fairSttHit->GetTubeID())){
669 // actTrack.correctSttSkewedHit(fairSttHit, fStrawMap.GetTube(fairSttHit->GetTubeID()));
670  continue;
671  }
672 
673 
674 // PndSttHit* correctedSttHit = new PndSttHit(CalculateCorrectedSTTHit(*fairSttHit, actTrack));
675 // correctedSttHit->SetDxyz(fairSttHit->GetDx(), fairSttHit->GetDy(), fairSttHit->GetDz());
676 // correctedSttHit->SetEntryNr(FairLink(FairRootManager::Instance()->GetBranchId("STTHit"), i));
677  PndRiemannHit actHit(fairSttHit, fMapIDtoHit[fairSttHit->GetEntryNr()]);
678 
679 
680 
681  if (fVerbose > 2) {
682  std::cout << std::endl;
683  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: Act Hit: " << fairSttHit->GetEntryNr() << std::endl;
684  }
685 
686 // if (CheckRiemannHitStt(&actTrack, &actHit,(FairHit*)fSttHits->At(i)) != true) continue;
687  if (CheckRiemannHitStt(&actTrack, &actHit,fairSttHit) != true) continue;
688 
689  if (fVerbose > 2 )
690  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: STTHit add to Track!" << std::endl;
691 
692  actTrack.addHit(actHit);
693  startTrack.insert(fMapIDtoHit[actHit.hit()->GetEntryNr()]);
694  actTrack.refit(kFALSE);
695 // actTrack.szFit(kFALSE);
696  nHitsStt++;
697  if (fStrawMap.IsSectorBorderStraw(fairSttHit->GetTubeID()) != 0){
698  if (fVerbose > 2)
699  std::cout << " TubeID: " << fairSttHit->GetTubeID() << " Sector: " << fStrawMap.GetSector(fairSttHit->GetTubeID()) << std::endl;
700  if (nextSectorDirection == 0){
701  nextSectorDirection = fStrawMap.IsSectorBorderStraw(fairSttHit->GetTubeID());
702  } else if(nextSectorDirection != fStrawMap.IsSectorBorderStraw(fairSttHit->GetTubeID())){
703  continue;
704  }
705  searchSector = fStrawMap.GetSector(fairSttHit->GetTubeID()) + fStrawMap.IsSectorBorderStraw(fairSttHit->GetTubeID());
706  searchSector %= 6;
707  if (searchSector < 0)
708  searchSector = 5;
709 
710  if (fVerbose > 2 )
711  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: Next Sector: " << searchSector << " " << fStrawMap.IsSectorBorderStraw(fairSttHit->GetTubeID())<< std::endl;
712 
713  continueSearch = true;
714  }
715 
716  }
717  if (nHitsStt == 0){ //if you have not found any sttHits and you have not looked yet into the left sector from the start sector continue otherwise stop.
718  if (searchSector != fStrawMap.GetLeftSector(startSector))
719  continueSearch = true;
720  else
721  continueSearch = false;
722  }
723  }
724 }
725 
726 void PndMvdSttGemRiemannTrackFinder::AssignSkewedSttHits(PndRiemannTrack& actTrack, std::set<Int_t>& startTrack)
727 {
728  int startSector = GetStrawSector(actTrack);
729  int searchSector = startSector;
730  int nHitsStt = -1;
731 
732 // SetVerbose(3);
733 
734  if (fVerbose > 2)
735  std::cout << "StartSector: " << searchSector << std::endl;
736 
737  std::vector<FairHit*> sttHits = fSttSkewedHitsInSectors[searchSector];
738  int leftSector = fStrawMap.GetLeftSector(startSector);
739  int rightSector = fStrawMap.GetRightSector(startSector);
740  sttHits.insert(sttHits.end(), fSttSkewedHitsInSectors[leftSector].begin(), fSttSkewedHitsInSectors[leftSector].end());
741  sttHits.insert(sttHits.end(), fSttSkewedHitsInSectors[rightSector].begin(), fSttSkewedHitsInSectors[rightSector].end());
742 
743  if (fVerbose > 2)
744  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: Actual searchSector: " << searchSector << std::endl;
745 
746  for (size_t i = 0; i < sttHits.size(); i++){
747 
748  if (nHitsStt < 0) nHitsStt = 0;
749 
750  FairHit* fairSttHit = sttHits[i];
751 
752  PndRiemannHit actHit(fairSttHit, fMapIDtoHit[fairSttHit->GetEntryNr()]);
753 
754 
755 
756  if (fVerbose > 2) {
757  std::cout << std::endl;
758  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: Act Hit: " << fairSttHit->GetEntryNr() << std::endl;
759  }
760 
761  if (CheckRiemannHitSkewedStt(&actTrack, &actHit,fairSttHit) != true) continue;
762 
763  if (fVerbose > 2 )
764  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: STTHit add to Track!" << std::endl;
765 
766  actTrack.addHit(actHit);
767  startTrack.insert(fMapIDtoHit[actHit.hit()->GetEntryNr()]);
768  actTrack.refit(kFALSE);
769 // actTrack.szFit(kFALSE);
770  nHitsStt++;
771 
772  }
773 // SetVerbose(0);
774 }
775 
776 void PndMvdSttGemRiemannTrackFinder::AssignGemHits(PndRiemannTrack& actTrack, set<Int_t>& startTrack){
777 
778  if (fVerbose > 1) std::cout << "NGemHits: " << fGemHits->GetEntriesFast() << std::endl;
779  for (int i = 0; i < fGemHits->GetEntriesFast(); i++){
780  FairHit* fairGemHit = (FairHit*)fGemHits->At(i);
781  PndRiemannHit actHit(fairGemHit, fMapIDtoHit[fairGemHit->GetEntryNr()]);
782 
783  if (CheckRiemannHitGem(&actTrack, &actHit,fairGemHit) != true) continue;
784 
785  if (fVerbose > 1 )
786  std::cout << "-I- PndMvdSttGemRiemannTrackFinder SearchSector: GemHit add to Track!" << std::endl;
787 
788  actTrack.addHit(actHit);
789  startTrack.insert(fMapIDtoHit[actHit.hit()->GetEntryNr()]);
790  }
791 }
792 
793 bool PndMvdSttGemRiemannTrackFinder::CheckSZ(PndRiemannTrack aTrack)
794 {
795  aTrack.szFit(false);
796  double r = aTrack.r();
797  double dip=aTrack.dip();
798  bool sign;
799  if ((aTrack.getHit(0)->z())>0 )
800  sign=true;
801  else sign=false;
802  if (aTrack.szChi2() > GetMaxSZChi2(r,dip,sign) || aTrack.szChi2() < 0){
803  if (fVerbose > 1) std::cout << "sz-Fit does not match, Chi2: " << aTrack.szChi2() << " max: " << GetMaxSZChi2(r,dip,sign) << std::endl;
804  return false;
805  }
806  return true;
807 }
808 
809 bool PndMvdSttGemRiemannTrackFinder::CheckRiemannHit(PndRiemannTrack* track, PndRiemannHit* hit, FairHit* fairHit)
810 {
811  if (fairHit->GetEntryNr().GetType() == FairRootManager::Instance()->GetBranchId("MVDHitsStrip") ||
812  fairHit->GetEntryNr().GetType() == FairRootManager::Instance()->GetBranchId("MVDHitsPixel") ){
813  return CheckRiemannHitMvd(track, hit, fairHit);
814  } else if (fairHit->GetEntryNr().GetType() == FairRootManager::Instance()->GetBranchId("GEMHit")) {
815  return CheckRiemannHitGem(track, hit, fairHit);
816  } else if (fairHit->GetEntryNr().GetType() == FairRootManager::Instance()->GetBranchId("STTHit")) {
817  return CheckRiemannHitStt(track, hit, fairHit);
818  } else {
819  std::cout << "-E- PndMvdRIemannTrackFinder::CheckRiemannHit: Detector " << FairRootManager::Instance()->GetBranchName(fairHit->GetEntryNr().GetType()) <<
820  " not supported!" << std::endl;
821  }
822  return false;
823 }
824 
825 bool PndMvdSttGemRiemannTrackFinder::CheckRiemannHitStt(PndRiemannTrack* track, PndRiemannHit* hit, FairHit* fairHit)
826 {
827  PndRiemannHit tempHit = track->correctSttHit((PndSttHit*)fairHit);
828 
829 // std::cout << "CheckRiemannHit PndRiemannHit: " << fairHit->GetEntryNr() << " : " << hit->x().X() << "/" << hit->x().Y() << "/" << hit->x().Z() << " " << hit->z() << std::endl;
830 // std::cout << "CheckRiemannHit PndRiemannHit: " << tempHit.x().X() << "/" << tempHit.x().Y() << "/" << tempHit.x().Z() << " " << tempHit.z() << std::endl;
831 
832  hit->setXYZ(tempHit.x().X(), tempHit.x().Y(), tempHit.x().Z());
833  hit->setDXYZ(tempHit.sigmaX(), tempHit.sigmaY(), 200);
834 
835 // std::cout << "CheckRiemannHit PndRiemannHit: " << hit->x().X() << "/" << hit->x().Y() << "/" << hit->x().Z() << " " << hit->z() << std::endl;
836 
837  double dist = track->dist(hit);
838 // std::cout << "CheckRiemannHit PndRiemannHit:Dist " << dist << std::endl;
839 
840  double maxDist = 1;
841 // double szDist = track->szDist(hit);
842 // double maxSZDist = 10;
843 // double szChi2 = track->calcSZChi2(hit);
844 // double maxSZChi2 = 25;
845  //double r = track->r(); //[R.K. 01/2017] unused variable?
846  //bool sign; //[R.K.02/2017] Unused variable?
847  //if ((track->getHit(1)->z())>0 ) //[R.K.02/2017] Unused variable?
848  //sign=true; //[R.K.02/2017] Unused variable?
849  //else sign=false; //[R.K.02/2017] Unused variable?
850  if (fVerbose > 1) std::cout << "STTHit: dist " << dist << std::endl;
851 
852  if (fabs(dist) > maxDist){
853  if (fVerbose > 1)
854  std::cout << "dist larger than " << maxDist << std::endl;
855  return false;
856  } else {
857 
858 // std::cout << "PndMvdSttGemRiemannTrackFinder::CheckRiemannHitStt Hit taken! Dist: " << dist << " < " << maxDist << std::endl << std::endl;
859  }
860 //
861 // if (szChi2 > maxSZChi2){
862 // if (fVerbose > 1) std::cout << " SZ Chi2 too big! Cut at: " << maxSZChi2 << std::endl;
863 // return false;
864 // }
865 // if (fabs(szDist) > maxSZDist){
866 // if (fVerbose > 1) std::cout << "SZ Dist too big! Cut at: " << maxSZDist << std::endl;
867 // return false;
868 // }
869  return true;
870 }
871 
872 bool PndMvdSttGemRiemannTrackFinder::CheckRiemannHitSkewedStt(PndRiemannTrack* track, PndRiemannHit* hit, FairHit* ) // fairHit //[R.K.03/2017] unused variable(s)
873 {
874 
875 
876 // std::cout << "CheckRiemannHit PndRiemannHit: " << hit->x().X() << "/" << hit->x().Y() << "/" << hit->x().Z() << " " << hit->z() << std::endl;
877 
878  double dist = track->dist(hit);
879 // std::cout << "CheckRiemannHit PndRiemannHit:Dist " << dist << std::endl;
880 
881  double maxDist = 1;
882 // double szDist = track->szDist(hit);
883 // double maxSZDist = 10;
884 // double szChi2 = track->calcSZChi2(hit);
885 // double maxSZChi2 = 25;
886  //double r = track->r(); //[R.K. 01/2017] unused variable?
887  //bool sign; //[R.K.02/2017] Unused variable?
888  //if ((track->getHit(1)->z())>0 ) //[R.K.02/2017] Unused variable?
889  //sign=true; //[R.K.02/2017] Unused variable?
890  //else sign=false; //[R.K.02/2017] Unused variable?
891  if (fVerbose > 1) std::cout << "STTHit: dist " << dist << std::endl;
892 
893  if (fabs(dist) > maxDist){
894  if (fVerbose > 1)
895  std::cout << "dist larger than " << maxDist << std::endl;
896  return false;
897  } else {
898 
899 // std::cout << "PndMvdSttGemRiemannTrackFinder::CheckRiemannHitStt Hit taken! Dist: " << dist << " < " << maxDist << std::endl << std::endl;
900  }
901 //
902 // if (szChi2 > maxSZChi2){
903 // if (fVerbose > 1) std::cout << " SZ Chi2 too big! Cut at: " << maxSZChi2 << std::endl;
904 // return false;
905 // }
906 // if (fabs(szDist) > maxSZDist){
907 // if (fVerbose > 1) std::cout << "SZ Dist too big! Cut at: " << maxSZDist << std::endl;
908 // return false;
909 // }
910  return true;
911 }
912 
913 bool PndMvdSttGemRiemannTrackFinder::CheckRiemannHitGem(PndRiemannTrack* track, PndRiemannHit* hit, FairHit* fairHit)
914 {
915  double dist = track->dist(hit);
916  double maxDist = 2;
917  double szDist = track->szDist(hit);
918  double maxSZDist = 10;
919  double szChi2 = track->calcSZChi2(hit);
920  double maxSZChi2 = 25;
921  //double r = track->r(); //[R.K. 01/2017] unused variable?
922  //double dip=track->dip(); //[R.K. 01/2017] unused variable?
923  //bool sign; //[R.K.02/2017] Unused variable?
924  //if ((track->getHit(1)->z())>0 ) //[R.K.02/2017] Unused variable?
925  //sign=true; //[R.K.02/2017] Unused variable?
926  //else sign=false; //[R.K.02/2017] Unused variable?
927  if (fVerbose > 1) std::cout << "Checking GemHit: " << fairHit->GetEntryNr() << std::endl;
928  if (fVerbose > 1) std::cout << ": dist " << dist << " szDist " << szDist << " szChi2 " << szChi2 << std::endl;
929 
930  if (fabs(dist) > maxDist){
931  if (fVerbose > 1) std::cout << "dist larger than " << maxDist << std::endl;
932  return false;
933  }
934 
935  if (szChi2 > maxSZChi2){
936  if (fVerbose > 1) std::cout << " SZ Chi2 too big! Cut at: " << maxSZChi2 << std::endl;
937  return false;
938  }
939  if (fabs(szDist) > maxSZDist){
940  if (fVerbose > 1) std::cout << "SZ Dist too big! Cut at: " << maxSZDist << std::endl;
941  return false;
942  }
943  return true;
944 }
945 
946 bool PndMvdSttGemRiemannTrackFinder::CheckRiemannHitMvd(PndRiemannTrack* track, PndRiemannHit* hit, FairHit* ) // fairHit //[R.K.03/2017] unused variable(s)
947 {
948  //double oldSZChi2 = track->szChi2(); //[R.K. 01/2017] unused variable?
949  double dist = track->dist(hit);
950  double szDist = track->szDist(hit);
951  double szChi2 = track->calcSZChi2(hit);
952  double r = track->r();
953  double dip=track->dip();
954  bool sign;
955  if ((track->getHit(1)->z())>0 )
956  sign=true;
957  else sign=false;
958  if (fVerbose > 1) std::cout << ": dist " << dist << " szDist " << szDist << " szChi2 " << szChi2 << std::endl;
959 
960  if (fabs(dist) > GetMaxPlaneDist(r,dip,sign)){
961  if (fVerbose > 1) std::cout << "dist larger than " << GetMaxPlaneDist(r,dip,sign) << std::endl;
962  return false;
963  }
964 
965 // if (fVerbose > 1) std::cout << "CheckRiemannHitMvd oldSZChi2 " << oldSZChi2 << " new SZChi2 " << szChi2 << std::endl;
966 
967 // if (oldSZChi2 < szChi2){
968 // if (fVerbose > 1) std::cout << " SZ Chi2 " << szChi2 << " bigger than before! " << oldSZChi2 << std::endl;
969 // return false;
970 // }
971 
972 // if (szChi2 > GetMaxSZChi2(r,dip,sign)){
973 // if (fVerbose > 1) std::cout << " SZ Chi2 too big! Cut at: " << GetMaxSZChi2(r,dip,sign) << std::endl;
974 // return false;
975 // }
976  if (fabs(szDist) > fMaxSZDist){
977  if (fVerbose > 1) std::cout << "SZ Dist too big! Cut at: " << fMaxSZDist << std::endl;
978  return false;
979  }
980  return true;
981 }
982 
983 bool PndMvdSttGemRiemannTrackFinder::CheckBoarderHitsStt(PndTrackCand* track)
984 {
985  FairMultiLinkedData sttHitsInTrack = track->GetLinksWithType(FairRootManager::Instance()->GetBranchId("STTHit"));
986 
987  for (int i = 0; i < sttHitsInTrack.GetNLinks(); i++){
988  PndSttHit* mySttHit = (PndSttHit*)fHits[fMapIDtoHit[sttHitsInTrack.GetLink(i)]];
989  if (fVerbose > 1) std::cout << "-I- PndMvdSttGemRiemannTrackFinder::CheckBoarderHitsStt Stt Index: " << mySttHit->GetTubeID() << std::endl;
990  if (fStrawMap.GetRow(mySttHit->GetTubeID()) == 0)
991  return true;
992  }
993  return false;
994 }
995 
996 int PndMvdSttGemRiemannTrackFinder::GetStrawSector(PndRiemannTrack& track)
997 {
998  double radiusSq(.0), oldRadiusSq(.0);
999  int indexHitMaxRadius(-1);
1000 
1001 
1002  std::vector<PndRiemannHit> riemannHits = track.getHits();
1003 
1004  for (size_t i = 0; i < riemannHits.size(); i++){
1005  FairLink hitLink = fMapHitToID[riemannHits[i].hitID()];
1006  if (fVerbose > 1) std::cout << "HitID Track : " << i << " : " << riemannHits[i].hitID() << " " << hitLink << std::endl;
1007  int branchId = hitLink.GetType();
1008  FairRootManager* man = FairRootManager::Instance();
1009  if (branchId == man->GetBranchId("MVDHitsStrip") || branchId == man->GetBranchId("MVDHitsPixel")){
1010  FairHit* FairHit = (PndSdsHit*)fHits[fMapIDtoHit[hitLink]];
1011  radiusSq = FairHit->GetX()*FairHit->GetX() + FairHit->GetY()*FairHit->GetY();
1012  if (fVerbose > 1) std::cout << "-I- PndMvdSttGemRiemannTrackFinder::GetStrawSector FairHit: " << i << " : " << radiusSq << " " << *FairHit << std::endl;
1013 
1014  if (radiusSq > oldRadiusSq){
1015  oldRadiusSq = radiusSq;
1016  indexHitMaxRadius = fMapIDtoHit[hitLink];
1017  }
1018  }
1019  }
1020  if (indexHitMaxRadius > -1){
1021  if (fVerbose > 1) std::cout << "-I- PndMvdSttGemRiemannTrackFinder::GetStrawSector IndexHitMaxRadius: " << indexHitMaxRadius << " Link: " << fMapHitToID[indexHitMaxRadius] << std::endl;
1022  FairHit* FairHit = (PndSdsHit*)fHits[indexHitMaxRadius];
1023  TVector3 hitPos(FairHit->GetX(), FairHit->GetY(), FairHit->GetZ());
1024  if (fVerbose > 1) std::cout << "-I- PndMvdSttGemRiemannTrackFinder::GetStrawSector Phi: " << hitPos.Phi() << " PhiSector: " << fStrawMap.FindPhiSector(hitPos.Phi()) << std::endl;
1025  return fStrawMap.FindPhiSector(hitPos.Phi());
1026  }
1027  return -1;
1028 }
1029 
1030 std::set<int> PndMvdSttGemRiemannTrackFinder::GetTooCloseHitsInLayer(int LayerNumber , int HitNumber )
1031 {
1032  std::set<int> result;
1033  int testN;
1034  for(unsigned int i=0;i<fLayers[LayerNumber].size();i++){
1035  testN=fLayers[LayerNumber][i];
1036  if (CheckHitDistance(HitNumber, testN)!= true){
1037  if (testN != HitNumber){
1038  result.insert(testN);
1039  }
1040  }
1041  }
1042  return result;
1043 
1044 }
1046 double PndMvdSttGemRiemannTrackFinder::GetMaxPlaneDist(double radius, double dip , bool sign)
1047 {
1048  double Pt=((radius/100)*2*3*1E8)/1E9;
1049  double Theta;
1050  if (fabs(TMath::ACos(dip)) < 1E-100){
1051  Theta = 0;
1052  }
1053  else {
1054  if (sign) Theta=(TMath::ATan(TMath::Power(TMath::Tan(TMath::ACos(dip)),-1)))*180/TMath::Pi();
1055  else Theta=(TMath::Pi()-TMath::ATan(TMath::Power(TMath::Tan(TMath::ACos(dip)),-1)))*180/TMath::Pi();
1056  }
1057 
1058  if(fCutDistH!=NULL){
1059  double minPt=fCutDistH->GetXaxis()->GetXmin();
1060  double maxPt=fCutDistH->GetXaxis()->GetXmax();
1061  double minTh=fCutDistH->GetYaxis()->GetXmin();
1062  double maxTh=fCutDistH->GetYaxis()->GetXmax();
1063 
1064  int binPt=int(floor((Pt-minPt)*fCutDistH->GetXaxis()->GetNbins()/(maxPt-minPt)))+1;
1065  int binTh=int(floor((Theta-minTh)*fCutDistH->GetYaxis()->GetNbins()/(maxTh-minTh)))+1;
1066 
1067  if (binPt<1) binPt=1;
1068  if (binPt>fCutDistH->GetXaxis()->GetNbins()) binPt=fCutDistH->GetXaxis()->GetNbins();
1069  if (binTh<1) binTh=1;
1070  if (binTh>fCutDistH->GetYaxis()->GetNbins()) binTh=fCutDistH->GetYaxis()->GetNbins();
1071  return fCutDistH->GetBinContent(binPt,binTh);
1072 
1073  }
1074  else
1075  return fMaxPlaneDist;
1076 }
1077 
1078 double PndMvdSttGemRiemannTrackFinder::GetMaxSZChi2(double radius, double dip , bool sign)
1079 {
1080  double Pt=((radius/100)*2*3*1E8)/1E9;
1081  double Theta;
1082  if (fabs(TMath::ACos(dip)) < 1E-100){
1083  Theta = 0;
1084  }
1085  else {
1086  if (sign) Theta=(TMath::ATan(TMath::Power(TMath::Tan(TMath::ACos(dip)),-1)))*180/TMath::Pi();
1087  else Theta=(TMath::Pi()-TMath::ATan(TMath::Power(TMath::Tan(TMath::ACos(dip)),-1)))*180/TMath::Pi();
1088  }
1089  if(fCutChi2H!=NULL){
1090  double minPt=fCutChi2H->GetXaxis()->GetXmin();
1091  double maxPt=fCutChi2H->GetXaxis()->GetXmax();
1092  double minTh=fCutChi2H->GetYaxis()->GetXmin();
1093  double maxTh=fCutChi2H->GetYaxis()->GetXmax();
1094 
1095  int binPt=int(floor((Pt-minPt)*fCutChi2H->GetXaxis()->GetNbins()/(maxPt-minPt)))+1;
1096  int binTh=int(floor((Theta-minTh)*fCutChi2H->GetYaxis()->GetNbins()/(maxTh-minTh)))+1;
1097 
1098  if (binPt<1) binPt=1;
1099  if (binPt>fCutChi2H->GetXaxis()->GetNbins()) binPt=fCutChi2H->GetXaxis()->GetNbins();
1100  if (binTh<1) binTh=1;
1101  if (binTh>fCutChi2H->GetYaxis()->GetNbins()) binTh=fCutChi2H->GetYaxis()->GetNbins();
1102  return fCutChi2H->GetBinContent(binPt,binTh);
1103 
1104  }
1105  else
1106  return fMaxSZChi2;
1107 }
1109 
1110 
PndRiemannTrackFinder::fHitsTooClose
std::vector< std::set< int > > fHitsTooClose
matrix of TrackNr and hits which are too close to one of the three starting points ...
Definition: PndRiemannTrackFinder.h:54
PndRiemannTrackFinder::fCurvAndDipOfCand
std::vector< std::pair< double, double > > fCurvAndDipOfCand
Curvature and dip of fPndTrackCand.
Definition: PndRiemannTrackFinder.h:52
PndMvdSttGemRiemannTrackFinder::GetTooCloseHitsInLayer
std::set< int > GetTooCloseHitsInLayer(int LayerNumber, int HitNumber)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:1030
PndRiemannTrack::getNumHits
unsigned int getNumHits()
Definition: PndRiemannTrack.h:74
PndRiemannTrack::calcSZChi2
double calcSZChi2(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:522
PndMvdSttGemRiemannTrackFinder::CheckRiemannHitMvd
bool CheckRiemannHitMvd(PndRiemannTrack *track, PndRiemannHit *hit, FairHit *fairHit)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:946
r
double r
Definition: RiemannTest.C:14
PndTrackCand::Print
void Print() const
Definition: PndTrackCand.cxx:164
i
Int_t i
Definition: run_full.C:25
PndMvdSttGemRiemannTrackFinder::GetMaxPlaneDist
double GetMaxPlaneDist(double radius, double dip, bool sign)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:1046
PndMvdSttGemRiemannTrackFinder::CheckRiemannHit
bool CheckRiemannHit(PndRiemannTrack *track, PndRiemannHit *hit, FairHit *fairHit)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:809
PndMvdSttGemRiemannTrackFinder::CheckRiemannHitStt
bool CheckRiemannHitStt(PndRiemannTrack *track, PndRiemannHit *hit, FairHit *fairHit)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:825
PndRiemannTrack::correctSttHit
PndRiemannHit correctSttHit(PndSttHit *mySttHit)
Definition: PndRiemannTrack.cxx:1174
PndRiemannTrack::szChi2
double szChi2() const
Definition: PndRiemannTrack.h:120
PndRiemannTrackFinder::fMaxSZChi2
double fMaxSZChi2
Maximum allowed Chi2 in an sz fit.
Definition: PndRiemannTrackFinder.h:60
PndMvdSttGemRiemannTrackFinder::fSttHitsInSectors
std::vector< std::vector< PndSttHit * > > fSttHitsInSectors
Definition: PndMvdSttGemRiemannTrackFinder.h:65
PndRiemannTrackFinder::fTracks
std::vector< PndRiemannTrack > fTracks
Resulting Riemann Tracks.
Definition: PndRiemannTrackFinder.h:49
PndRiemannTrackFinder::TrackExists
bool TrackExists(std::set< Int_t > hitsInTrack)
Definition: PndRiemannTrackFinder.cxx:506
PndMvdSttGemRiemannTrackFinder::fLayerMap
std::map< TString, int > fLayerMap
identifier string, assigned layer id
Definition: PndMvdSttGemRiemannTrackFinder.h:71
track
PndRiemannTrack track
Definition: RiemannTest.C:33
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndRiemannTrackFinder::fMapHitToID
std::map< int, FairLink > fMapHitToID
map to convert the list of hits back into a FairLink
Definition: PndRiemannTrackFinder.h:55
s
TLorentzVector s
Definition: Pnd2DStar.C:50
PndMvdSttGemRiemannTrackFinder::CheckRiemannHitSkewedStt
bool CheckRiemannHitSkewedStt(PndRiemannTrack *track, PndRiemannHit *hit, FairHit *fairHit)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:872
PndRiemannTrackFinder::fMapIDtoHit
std::map< FairLink, int > fMapIDtoHit
map to convert the list of detID/hitID hits into the list of hits for track finding ...
Definition: PndRiemannTrackFinder.h:56
push_back
labels push_back("electron")
PndRiemannTrack::szDist
double szDist(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:900
n
int n
Definition: HypStatDecay.cxx:21
PndSttStrawMap::IsSectorBorderStraw
int IsSectorBorderStraw(int strawindex) const
Definition: PndSttStrawMap.cxx:128
PndSttStrawMap::GetSector
int GetSector(int strawindex) const
Definition: PndSttStrawMap.h:33
CAMath::Tan
float Tan(float x)
Definition: PndCAMath.h:165
PndRiemannTrack::dist
double dist(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:129
PndRiemannTrackFinder::fHits
std::vector< FairHit * > fHits
Vector of all FairHits used for track finding (fitting)
Definition: PndRiemannTrackFinder.h:48
PndGeoHandling::GetPath
TString GetPath(Int_t shortID)
for a given shortID the path is returned
Definition: PndGeoHandling.cxx:222
PndMvdSttGemRiemannTrackFinder::AddHits
void AddHits(TClonesArray *hits, Int_t branchId)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:95
PndRiemannHit::setDXYZ
void setDXYZ(double dx, double dy, double dz)
Definition: PndRiemannHit.cxx:70
PndRiemannTrackFinder::fMinNumberOfHits
int fMinNumberOfHits
Minimum number of hits in track necessary for a match.
Definition: PndRiemannTrackFinder.h:61
PndRiemannTrackFinder::fTrackCand
std::vector< PndTrackCand > fTrackCand
List of track candidates.
Definition: PndRiemannTrackFinder.h:51
p
Double_t p
Definition: anasim.C:58
PndRiemannTrack::getSZt
double getSZt() const
Definition: PndRiemannTrack.h:73
PndGemHit::GetSensorNr
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
PndTrackCand::AddHit
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Definition: PndTrackCand.cxx:36
PndSttStrawMap::GetRightSector
int GetRightSector(int sector) const
Definition: PndSttStrawMap.h:37
PndMvdSttGemRiemannTrackFinder::AssignSttHits
void AssignSttHits(PndRiemannTrack &actTrack, std::set< Int_t > &startTrack)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:636
PndSttStrawMap::IsSkewedStraw
bool IsSkewedStraw(int strawindex) const
Definition: PndSttStrawMap.h:43
PndRiemannHit::setXYZ
void setXYZ(double x, double y, double z)
Definition: PndRiemannHit.cxx:57
PndRiemannHit::sigmaX
double sigmaX() const
Definition: PndRiemannHit.h:78
PndRiemannTrack::refit
void refit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:172
PndMvdSttGemRiemannTrackFinder::fSttSkewedHitsInSectors
std::vector< std::vector< FairHit * > > fSttSkewedHitsInSectors
Definition: PndMvdSttGemRiemannTrackFinder.h:66
PndSttHit::GetTubeID
Int_t GetTubeID() const
Definition: PndSttHit.h:75
PndRiemannTrackFinder::fMaxPlaneDist
double fMaxPlaneDist
Distance cut between new point and riemann plane.
Definition: PndRiemannTrackFinder.h:57
PndMvdSttGemRiemannTrackFinder::AssignSkewedSttHits
void AssignSkewedSttHits(PndRiemannTrack &actTrack, std::set< Int_t > &startTrack)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:726
PndSttStrawMap::GetLeftSector
int GetLeftSector(int sector) const
Definition: PndSttStrawMap.h:36
PndRiemannTrackFinder::CheckHitInSameSensor
bool CheckHitInSameSensor(int hit1, int hit2)
Tests if hits in the same sensor are selected.
Definition: PndRiemannTrackFinder.cxx:305
PndRiemannTrackFinder::CheckHitDistance
bool CheckHitDistance(int hit1, int hit2)
Tests if the distance is larger than fMinPointDistance.
Definition: PndRiemannTrackFinder.cxx:238
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndRiemannTrack::getSZm
double getSZm() const
Definition: PndRiemannTrack.h:72
PndGeoHandling::Instance
static PndGeoHandling * Instance()
Definition: PndGeoHandling.cxx:36
PndGemHit::GetStationNr
Int_t GetStationNr() const
Definition: PndGemHit.h:81
PndSttStrawMap::FindPhiSector
int FindPhiSector(double phi) const
Definition: PndSttStrawMap.cxx:98
phi
double phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
PndRiemannTrack::r
double r() const
Definition: PndRiemannTrack.cxx:416
PndRiemannHit::z
double z() const
Definition: PndRiemannHit.cxx:87
PndMvdSttGemRiemannTrackFinder::CheckSZ
bool CheckSZ(PndRiemannTrack aTrack)
Tests the results of the sz fit.
Definition: PndMvdSttGemRiemannTrackFinder.cxx:793
hit
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
Definition: hit.C:1
PndRiemannTrackFinder::CreateRiemannTrack
PndRiemannTrack CreateRiemannTrack(std::set< Int_t > aHits)
Creates a PndRiemannTrack from an array of indices of Hits.
Definition: PndRiemannTrackFinder.cxx:292
PndRiemannTrackFinder::CheckZeroPassing
bool CheckZeroPassing(std::set< int > hitIds, int hit)
If the track contains (0,0) all points have to go forward or all have to go backward.
Definition: PndRiemannTrackFinder.cxx:259
PndRiemannTrack::orig
TVectorD orig() const
Definition: PndRiemannTrack.cxx:405
PndTrackCand::GetNHits
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
PndRiemannHit::x
const TVector3 & x() const
Definition: PndRiemannHit.h:71
TrackHits
int TrackHits()
Definition: TrackHits.C:8
ClassImp
ClassImp(PndMvdSttGemRiemannTrackFinder)
PndRiemannTrack::szFit
void szFit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:441
PndRiemannTrack::addHit
void addHit(PndRiemannHit &hit)
Definition: PndRiemannTrack.cxx:102
PndMvdSttGemRiemannTrackFinder::GetMaxSZChi2
double GetMaxSZChi2(double radius, double dip, bool sign)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:1078
PndMvdSttGemRiemannTrackFinder::CheckRiemannHitGem
bool CheckRiemannHitGem(PndRiemannTrack *track, PndRiemannHit *hit, FairHit *fairHit)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:913
PndSttStrawMap::GetRow
int GetRow(int strawindex) const
Definition: PndSttStrawMap.h:35
sign
int sign(T val)
Definition: PndCADef.h:48
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
PndMvdSttGemRiemannTrackFinder::AssignGemHits
void AssignGemHits(PndRiemannTrack &actTrack, std::set< Int_t > &startTrack)
Definition: PndMvdSttGemRiemannTrackFinder.cxx:776
PndMvdSttGemRiemannTrackFinder::FindTracks
void FindTracks()
Main function to start the riemann track finding.
Definition: PndMvdSttGemRiemannTrackFinder.cxx:259
PndRiemannTrack::getHits
std::vector< PndRiemannHit > getHits() const
Definition: PndRiemannTrack.h:77
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndRiemannHit::sigmaY
double sigmaY() const
Definition: PndRiemannHit.h:79
PndRiemannTrackFinder::fHitsInTracks
std::vector< std::set< Int_t > > fHitsInTracks
Vector of indizes which hits where used in which track.
Definition: PndRiemannTrackFinder.h:50
PndMvdSttGemRiemannTrackFinder::fLayers
std::vector< std::vector< int > > fLayers
contains layer information of hits
Definition: PndMvdSttGemRiemannTrackFinder.h:70
PndRiemannTrackFinder::fMaxSZDist
double fMaxSZDist
Distance cut between s-z coordinate of a new point and the sz-fit of the hits in the track...
Definition: PndRiemannTrackFinder.h:58
PndMvdSttGemRiemannTrackFinder::GetStartTracks
std::vector< std::set< Int_t > > GetStartTracks()
Definition: PndMvdSttGemRiemannTrackFinder.cxx:514
PndRiemannTrack::getHit
PndRiemannHit * getHit(unsigned int i)
Definition: PndRiemannTrack.h:75
PndRiemannTrackFinder::CheckHitInTrack
bool CheckHitInTrack(std::set< int > hitIds, int hit)
Check if this HitId is used in the track already.
Definition: PndRiemannTrackFinder.cxx:320
Generated on Sun Aug 18 2019 03:01:40 for FairRoot/PandaRoot by   doxygen 1.8.8