FairRoot/PandaRoot
PndStraightLineTrackFinderTask.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndStraightLineTrackFinderTask -----
3 // ----- Created 22/10/09 by M. Michel -----
4 // -------------------------------------------------------------------------
5 
6 #include "TClonesArray.h"
7 #include "TArrayD.h"
8 #include "TGeoManager.h"
9 
10 #include "FairRootManager.h"
11 #include "FairRun.h"
12 #include "FairRuntimeDb.h"
13 
14 #include "PndStraightLineTrackFinderTask.h"
15 
16 #include "PndSdsDigiStrip.h"
17 #include "TStopwatch.h"
18 // #include "PndSdsPixelCluster.h"
19 
20 
21 // ----- Default constructor -------------------------------------------
22 PndStraightLineTrackFinderTask::PndStraightLineTrackFinderTask(Int_t inFinderMode, TString hitBranch, TString clusterBranch, TString digiBranch, Int_t innSensPP) :
23  FairTask("LMD Track Finding Task"), nSensPP(innSensPP)
24 {
25  fFinderMode = inFinderMode;
26  fHitBranchStrip = hitBranch;
27  fClusterBranchStrip = clusterBranch;
28  fDigiBranchStrip = digiBranch;
29  dXY = 0.1;
30  flagStipSens = false;
31  flagPixelSens = true;
32  // dXY = 0.01;//TEST
33 }
34 
35 
36 // -------------------------------------------------------------------------
37 
38 
39 // ----- Destructor ----------------------------------------------------
40 PndStraightLineTrackFinderTask::~PndStraightLineTrackFinderTask()
41 {
42 }
43 // -------------------------------------------------------------------------
44 
45 // ----- Initialization of Parameter Containers -------------------------
46 void PndStraightLineTrackFinderTask::SetParContainers()
47 {
48  // Get Base Container
49 /*
50  FairRun* ana = FairRun::Instance();
51  FairRuntimeDb* rtdb=ana->GetRuntimeDb();
52  fGeoPar = (PndSdsGeoPar*)(rtdb->getContainer("PndSdsGeoPar"));
53 */
54 }
55 
56 InitStatus PndStraightLineTrackFinderTask::ReInit()
57 {
58 
59  InitStatus stat=kERROR;
60  return stat;
61 
62  /*
63  FairRun* ana = FairRun::Instance();
64  FairRuntimeDb* rtdb=ana->GetRuntimeDb();
65  fGeoPar=(PndSdsGeoPar*)(rtdb->getContainer("PndSdsGeoPar"));
66 
67  return kSUCCESS;
68  */
69 }
70 
71 // ----- Public method Init --------------------------------------------
72 InitStatus PndStraightLineTrackFinderTask::Init()
73 {
74  //lmddim = PndLmdDim::Instance();
75  // lmddim -> Read_transformation_matrices("matrices.txt", true);
76  //lmddim -> Read_transformation_matrices("matrices_perfect.txt", false);
77 
78  FairRootManager* ioman = FairRootManager::Instance();
79 
80  if ( ! ioman )
81  {
82  std::cout << "-E- PndStraightLineTrackFinderTask::Init: "
83  << "RootManager not instantiated!" << std::endl;
84  return kFATAL;
85  }
86 
87  // Get input array
88  fStripHitArray = (TClonesArray*) ioman->GetObject(fHitBranchStrip);
89  if ( !fStripHitArray){
90  std::cout << "-W- PndStraightLineTrackFinderTask::Init: " << "No fStripHitArray!" << std::endl;
91  return kERROR;
92  }
93 
94  fStripClusterArray = (TClonesArray*) ioman->GetObject(fClusterBranchStrip);
95  if ( !fStripClusterArray){
96  std::cout << "-W- PndStraightLineTrackFinderTask::Init: " << "No StripclusterArray!" << std::endl;
97  //return kERROR;
98  }
99 
100  fStripDigiArray = (TClonesArray*) ioman->GetObject(fDigiBranchStrip);
101  if ( !fStripDigiArray){
102  std::cout << "-W- PndStraightLineTrackFinderTask::Init: " << "No StripdigiArray!" << std::endl;
103  //return kERROR;
104  }
105 
106  fTrackCandArray = new TClonesArray("PndTrackCand");
107  ioman->Register("MVDTestBeamTrackCand", "Mvd", fTrackCandArray, kTRUE);
108 
109  // fTruePointArray = new TClonesArray("PndSdsMCPoint");
110 
111 
112 
113  std::cout << "-I- PndStraightLineTrackFinderTask: Initialisation successfull" << std::endl;
114  return kSUCCESS;
115 }
116 // -------------------------------------------------------------------------
117 
118 // ----- Private method SortHitsByDet --------------------------------------------
119 bool PndStraightLineTrackFinderTask::SortHitsByDet(std::vector< std::vector< std::pair<Int_t,bool> > > &hitsd, Int_t nStripHits)
120 {
121  Int_t nPlanes=0;
122 
123 //sort in plane's
124  for(Int_t iHit = 0; iHit < nStripHits; iHit++){
125  PndSdsHit* myHit = (PndSdsHit*)(fStripHitArray->At(iHit));
126 
127  Int_t sensid = myHit->GetSensorID(); // Sensors: 1..32
128  // cout<<"sensid = "<<sensid<<endl;
129  Int_t planeid = floor((sensid)/(double)nSensPP); //nSensPP sensors/plane => Planes: 0..3
130  hitsd.at(planeid).push_back( make_pair (iHit,false) );
131  }
132 
133  for(Int_t iPlane = 0; iPlane < 4; iPlane++){
134  if(hitsd.at(iPlane).size()>0) nPlanes++;
135  }
136 
137  // cout << "Hits: " << nStripHits << endl;
138  if(fVerbose>2) {
139  cout << "Hits: " << nStripHits << " in " << nPlanes << " plane(s)." << endl;
140  for(Int_t idet = 0; idet < 4; idet++)
141  cout << "Plane: "<< idet <<" DiscHits: "<< hitsd.at(idet).size() <<endl;
142  }
143 
144  if(nPlanes>2) return true;
145  return false;
146 }
147 // -------------------------------------------------------------------------
148 
149 // ----- Private method SortHitsByDet2 --------------------------------------------
150 bool PndStraightLineTrackFinderTask::SortHitsByDet2(std::vector< std::vector< std::pair<Int_t,bool> > > &hitsd, Int_t nStripHits)
151 {
152  Int_t nPlanes=0;
153 
154 //sort in plane's
155  for(Int_t iHit = 0; iHit < nStripHits; iHit++){
156  PndSdsHit* myHit = (PndSdsHit*)(fStripHitArray->At(iHit));
157 
158  Int_t sensid = myHit->GetSensorID();
159  //int ihalf,iplane,imodule,iside,idie,isensor;
160  //lmddim->Get_sensor_by_id(sensid,ihalf,iplane,imodule,iside,idie,isensor);
161  // hitsd.at(iplane).push_back(iHit);
162  hitsd.at(sensid).push_back( make_pair (iHit,false) );
163  // cout<<"sensid = "<<sensid<<endl;
164  // Int_t planeid = floor((sensid)/(double)nSensPP); //nSensPP sensors/plane => Planes: 0..3
165  // hitsd.at(planeid).push_back( make_pair (iHit,false) );
166  }
167 
168  for(Int_t iPlane = 0; iPlane < 4; iPlane++){
169  if(hitsd.at(iPlane).size()>0) nPlanes++;
170  }
171 
172  // cout << "Hits: " << nStripHits << endl;
173  if(fVerbose>2) {
174  cout << "Hits: " << nStripHits << " in " << nPlanes << " plane(s)." << endl;
175  for(Int_t idet = 0; idet < 4; idet++)
176  cout << "Plane: "<< idet <<" DiscHits: "<< hitsd.at(idet).size() <<endl;
177  }
178 
179  if(nPlanes>2) return true;
180  return false;
181 }
182 // -------------------------------------------------------------------------
183 
184 // ----- Private method SortHitsByZ --------------------------------------------
185 bool PndStraightLineTrackFinderTask::SortHitsByZ(std::vector< std::vector< std::pair<Int_t,bool> > > &hitsd, Int_t nStripHits)
186 {
187  std::vector<Double_t> detZ;
188 
189 //find plane positions
190  for (Int_t iHit = 0; iHit < nStripHits; iHit++){
191  Double_t tmp = ((PndSdsHit*) (fStripHitArray->At(iHit)))->GetZ();
192  bool newZ = true;
193  for(Int_t idet = 0; idet < detZ.size(); idet++){
194  // if(tmp == detZ.at(idet)){ //check if already found
195  // cout<<"tmp = "<<tmp<<" detZ.at(idet) = "<<detZ.at(idet)
196  // <<" fabs(tmp-detZ.at(idet))="<<fabs(tmp-detZ.at(idet))<<endl;
197  if(fabs(tmp-detZ.at(idet))<9.){ //check if already found [for using with Dipole]
198  newZ = false;
199  }
200  }
201  if(newZ){
202  detZ.push_back(tmp);
203 
204  //sort positions
205  Int_t pos=-1;
206  for(Int_t idet = detZ.size()-1; idet >= 0; idet--){
207  if(tmp < detZ.at(idet))
208  pos=idet;
209  }
210  if(pos!=-1){
211  Double_t swap = detZ.at(pos);
212  detZ.at(pos) = tmp;
213  tmp = swap;
214  for(Int_t i=pos+1; i<detZ.size(); i++){
215  swap = detZ.at(i);
216  detZ.at(i) = tmp;
217  tmp = swap;
218  }
219  }
220  }
221  }
222  // cout<<"Attention! detZ.size()="<<detZ.size()<<endl;
223 
224 //sort in plane's
225  for(Int_t iHit = 0; iHit < nStripHits; iHit++){
226  PndSdsHit* myHit = (PndSdsHit*)(fStripHitArray->At(iHit));
227 
228  Double_t z = myHit->GetZ();
229  for(Int_t idet = 0; idet < detZ.size(); idet++){ //planes
230 
231  // if( z == detZ.at(idet) ){
232  if( fabs(z-detZ.at(idet))<9. ){ //[for using with Dipole]
233  hitsd.at(idet).push_back( make_pair (iHit,false) );
234  // cout<<"detZ.at("<<idet<<")="<<detZ.at(idet)<<" z="<<z<<endl;
235  }
236  }
237  }
238 
239  // cout << "Hits: " << nStripHits << endl;
240  if(fVerbose>2) {
241  cout << "Hits: " << nStripHits << " in " << detZ.size() << " plane(s)." << endl;
242  for(Int_t idet = 0; idet < detZ.size(); idet++)
243  cout << "Plane: "<< idet <<" DiscHits: "<< hitsd.at(idet).size() <<endl;
244  }
245 
246  if(detZ.size()>3) return true;
247  return false;
248 }
249 // -------------------------------------------------------------------------
250 
251 // ----- Private method FindHitsIII --------------------------------------------
252 void PndStraightLineTrackFinderTask::FindHitsIII(std::vector<PndTrackCand> &tofill, std::vector< std::vector< std::pair<Int_t,bool> > > &hitsd, Int_t nStripHits)
253 {
254 
255  std::vector<TVector3> trackStart, trackVec;//pseudo tracks
256  std::vector<TVector3> trackStartd, trackVecd;//save errors
257  std::vector< Int_t > trackID2, trackID3; //for TrackCand
258 
259 //iterate first discs-hits with all seconds, save pseudo-tracks
260  TVector3 start, tmp, vec, dstart, dvec; //temp-vars
261 
262  if(hitsd.size()<3) return;
263  for (Int_t i=0; i<hitsd.at(1).size(); i++)
264  {
265  if(hitsd.at(1).at(i).second) continue; //if already used
266  PndSdsHit *hit1=(PndSdsHit*)fStripHitArray->At(hitsd.at(1).at(i).first);
267  start.SetXYZ(hit1->GetX(), hit1->GetY(), hit1->GetZ());
268  dstart.SetXYZ(hit1->GetDx(), hit1->GetDy(), hit1->GetDz());
269  for (Int_t k=0; k<hitsd.at(2).size(); k++)
270  {
271  if(hitsd.at(2).at(k).second) continue; //if already used
272  PndSdsHit *hit2=(PndSdsHit*)fStripHitArray->At(hitsd.at(2).at(k).first);
273  tmp.SetXYZ(hit2->GetX(), hit2->GetY(), hit2->GetZ());
274  vec = tmp - start; //calc direction vector for FINDING
275  dvec.SetXYZ(hit2->GetDx(), hit2->GetDy(), hit2->GetDz());
276  if(vec.Theta()>0.03 && vec.Theta()<0.05 && vec.Phi()>-0.25 && vec.Phi()<0.25){ //ignore vectors with theta outside 2-9 mrad
277  //if(0<1){
278  // if(vec.Theta()<0.01){ //ignore vectors with theta outside 10 mrad
279  trackStart.push_back(start);
280  trackStartd.push_back(dstart);
281  trackVec.push_back(vec); //save vector from start to second
282  trackVecd.push_back(dvec); //save error of second point for FIT, NOT error of direction vector
283  //trackID1.push_back(hitsd.at(0).at(i).first); //save Hit-Id's for TrackCand
284  //trackID2.push_back(hitsd.at(1).at(k).first);
285  trackID2.push_back(i); //save Hit-Id's for TrackCand
286  trackID3.push_back(k);
287  }
288  }//end of disc 2 (or 3 if none in disc 2) hits
289  }//end of disc 1 hits
290 
291  if(fVerbose>1) cout << "Pseudos L2L3: "<< trackStart.size() <<endl;
292 
293  std::vector<Int_t> ids;
294 
295 //check if other discs have hits in track, add points for fitting ---------------
296  for (Int_t i=0; i<trackStart.size(); i++) //pseudo-loop
297  {
298  ids.clear();
299  Int_t pntcnt=2;
300 
301  start = trackStart.at(i);
302  dstart = trackStartd.at(i);
303  vec = trackVec.at(i);
304  dvec = trackVecd.at(i);
305  ids.push_back(hitsd.at(1).at(trackID2.at(i)).first);
306  ids.push_back(hitsd.at(2).at(trackID3.at(i)).first);
307  std::vector< std::pair<Int_t,Int_t> > otherIDs;
308 
309  for (Int_t idet=3; idet < 4; idet++){ //i know this is not a loop, but in case we someday will have more planes
310  Double_t distClosest = 2*idet*dXY; //just bigger as possible
311  Bool_t firstp = true;
312  for (Int_t ihit=0; ihit<hitsd.at(idet).size(); ihit++)
313  {
314  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(idet).at(ihit).first);
315  Double_t scale = (hit->GetZ()-start.z())/vec.z();
316  tmp = start + scale*vec; //extend search-vector to hit-plane
317  Double_t distTmp = sqrt((tmp.x()-hit->GetX())*(tmp.x()-hit->GetX()) + (tmp.y()-hit->GetY())*(tmp.y()-hit->GetY()));
318  if( distTmp<(idet*dXY) ){ //if in diameter
319  if(firstp){
320  pntcnt++;
321  ids.push_back(hitsd.at(idet).at(ihit).first);
322  otherIDs.push_back( make_pair (idet,ihit) );
323  distClosest=distTmp;
324  firstp=false;
325  }else{
326  if(distTmp<distClosest){ //change if nearer
327  ids.pop_back();
328  otherIDs.pop_back();
329  ids.push_back(hitsd.at(idet).at(ihit).first);
330  otherIDs.push_back( make_pair (idet,ihit) );
331  }
332  }
333  }//endif in diameter
334  }//end of disc n hits
335  }//end of discs
336 
337  if(fVerbose>2) cout << " Track L2L3: "<< i << "#Planes: " << ids.size() <<endl;
338 
339  //create track for fitting
340  if(ids.size()>2){ //third hit found <=> !track found!
341  PndTrackCand *myTCand = new PndTrackCand();
342 
343  for (Int_t id=0; id<ids.size(); id++){
344  PndSdsHit* myHit = (PndSdsHit*)(fStripHitArray->At(ids.at(id)));
345  PndSdsClusterStrip* myCluster;
346  PndSdsDigiStrip* astripdigi;
347  if ( !fStripClusterArray || !fStripDigiArray){ //for ideal/fast Hit-Reco
348  myTCand->AddHit(0,ids.at(id),myHit->GetPosition().Mag());
349  }else{
350  myCluster = (PndSdsClusterStrip*)(fStripClusterArray->At(myHit->GetClusterIndex()));
351  astripdigi = (PndSdsDigiStrip*)fStripDigiArray->At(myCluster->GetDigiIndex(0));
352  myTCand->AddHit(astripdigi->GetDetID(),ids.at(id),myHit->GetPosition().Mag());
353  }
354  }
355 
356  //mark used hits---------
357  hitsd.at(1).at(trackID2.at(i)).second=true;
358  hitsd.at(2).at(trackID3.at(i)).second=true;
359  for(Int_t id=0; id<otherIDs.size(); id++){
360  hitsd.at(otherIDs.at(id).first).at(otherIDs.at(id).second).second=true;
361  }
362 
363  // ///Add seed information to track------------
364  // PndSdsHit* myHit0 = (PndSdsHit*)(fStripHitArray->At(ids.at(0)));
365  // PndSdsHit* myHit1 = (PndSdsHit*)(fStripHitArray->At(ids.at(1)));
366  // TVector3 hit0 = myHit0->GetPosition(); TVector3 hit1 = myHit1->GetPosition();
367  // TVector3 posSeed(hit0.X(),hit0.Y(),hit0.Z());
368  // vec*=1./vec.Mag();
369  // //shift trk out of plane [needed for correct treatment in Kalman Fillter and GEANE]
370  // double sh_z = -0.035; //350 mkm
371  // double sh_x = vec.X()*sh_z;
372  // double sh_y = vec.Y()*sh_y;
373  // TVector3 sh_point(sh_x,sh_y,sh_z);
374  // posSeed +=sh_point;
375  // myTCand->setTrackSeed(posSeed,vec,-1);
376  // ///-------------------------------------
377 
378  tofill.push_back(*(myTCand)); //save Track Candidate
379  //new((*fTrackCandArray)[trackCnt]) PndTrackCand(*(myTCand));
380  delete myTCand;
381  }//Track Cand build
382  }//end of pseudo-tracks
383 
384 }
385 // -------------------------------------------------------------------------
386 
387 // ----- Private method FindHitsII --------------------------------------------
388 void PndStraightLineTrackFinderTask::FindHitsII(std::vector<PndTrackCand> &tofill, std::vector< std::vector< std::pair<Int_t,bool> > > &hitsd, Int_t nStripHits)
389 {
390 
391  std::vector<TVector3> trackStart, trackVec;//pseudo tracks
392  std::vector<TVector3> trackStartd, trackVecd;//save errors
393  std::vector< Int_t > trackID1, trackID3; //for TrackCand
394 
395 //iterate first discs-hits with all seconds, save pseudo-tracks
396  TVector3 start, tmp, vec, dstart, dvec; //temp-vars
397 
398  if(hitsd.size()<3) return;
399  for (Int_t i=0; i<hitsd.at(0).size(); i++)
400  {
401  if(hitsd.at(0).at(i).second) continue; //if already used
402  PndSdsHit *hit1=(PndSdsHit*)fStripHitArray->At(hitsd.at(0).at(i).first);
403  start.SetXYZ(hit1->GetX(), hit1->GetY(), hit1->GetZ());
404  dstart.SetXYZ(hit1->GetDx(), hit1->GetDy(), hit1->GetDz());
405  for (Int_t k=0; k<hitsd.at(2).size(); k++)
406  {
407  if(hitsd.at(2).at(k).second) continue; //if already used
408  PndSdsHit *hit2=(PndSdsHit*)fStripHitArray->At(hitsd.at(2).at(k).first);
409  tmp.SetXYZ(hit2->GetX(), hit2->GetY(), hit2->GetZ());
410  vec = tmp - start; //calc direction vector for FINDING
411  dvec.SetXYZ(hit2->GetDx(), hit2->GetDy(), hit2->GetDz());
412  if(vec.Theta()>0.03 && vec.Theta()<0.05 && vec.Phi()>-0.25 && vec.Phi()<0.25){ //ignore vectors with theta outside 2-9 mrad
413  // if(0<1){
414  // if(vec.Theta()<0.01){ //ignore vectors with theta outside 10 mrad
415  trackStart.push_back(start);
416  trackStartd.push_back(dstart);
417  trackVec.push_back(vec); //save vector from start to second
418  trackVecd.push_back(dvec); //save error of second point for FIT, NOT error of direction vector
419  //trackID1.push_back(hitsd.at(0).at(i).first); //save Hit-Id's for TrackCand
420  //trackID2.push_back(hitsd.at(1).at(k).first);
421  trackID1.push_back(i); //save Hit-Id's for TrackCand
422  trackID3.push_back(k);
423  }
424  }//end of disc 2 (or 3 if none in disc 2) hits
425  }//end of disc 1 hits
426 
427  if(fVerbose>1) cout << "Pseudos L1L3: "<< trackStart.size() <<endl;
428 
429  std::vector<Int_t> ids;
430 
431 //check if other discs have hits in track, add points for fitting ---------------
432  for (Int_t i=0; i<trackStart.size(); i++) //pseudo-loop
433  {
434  ids.clear();
435  Int_t pntcnt=2;
436 
437  start = trackStart.at(i);
438  dstart = trackStartd.at(i);
439  vec = trackVec.at(i);
440  dvec = trackVecd.at(i);
441  ids.push_back(hitsd.at(0).at(trackID1.at(i)).first);
442  ids.push_back(hitsd.at(2).at(trackID3.at(i)).first);
443  std::vector< std::pair<Int_t,Int_t> > otherIDs;
444 
445  for (Int_t idet=3; idet < 4; idet++){ //i know this is not a loop, but in case we someday will have more planes
446  Double_t distClosest = 2*idet*dXY; //just bigger as possible
447  Bool_t firstp = true;
448  for (Int_t ihit=0; ihit<hitsd.at(idet).size(); ihit++)
449  {
450  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(idet).at(ihit).first);
451  Double_t scale = (hit->GetZ()-start.z())/vec.z();
452  tmp = start + scale*vec; //extend search-vector to hit-plane
453  Double_t distTmp = sqrt((tmp.x()-hit->GetX())*(tmp.x()-hit->GetX()) + (tmp.y()-hit->GetY())*(tmp.y()-hit->GetY()));
454  if( distTmp<(idet*dXY) ){ //if in diameter
455  if(firstp){
456  pntcnt++;
457  ids.push_back(hitsd.at(idet).at(ihit).first);
458  otherIDs.push_back( make_pair (idet,ihit) );
459  distClosest=distTmp;
460  firstp=false;
461  }else{
462  if(distTmp<distClosest){ //change if nearer
463  ids.pop_back();
464  otherIDs.pop_back();
465  ids.push_back(hitsd.at(idet).at(ihit).first);
466  otherIDs.push_back( make_pair (idet,ihit) );
467  }
468  }
469  }//endif in diameter
470  }//end of disc n hits
471  }//end of discs
472 
473  if(fVerbose>2) cout << " Track L1L3: "<< i << "#Planes: " << ids.size() <<endl;
474 
475  //create track für fitting
476  if(ids.size()>2){ //third hit found <=> !track found!
477  PndTrackCand *myTCand = new PndTrackCand();
478 
479  for (Int_t id=0; id<ids.size(); id++){
480  PndSdsHit* myHit = (PndSdsHit*)(fStripHitArray->At(ids.at(id)));
481  PndSdsClusterStrip* myCluster;
482  PndSdsDigiStrip* astripdigi;
483  if ( !fStripClusterArray || !fStripDigiArray){ //for ideal/fast Hit-Reco
484  myTCand->AddHit(0,ids.at(id),myHit->GetPosition().Mag());
485  }else{
486  myCluster = (PndSdsClusterStrip*)(fStripClusterArray->At(myHit->GetClusterIndex()));
487  astripdigi = (PndSdsDigiStrip*)fStripDigiArray->At(myCluster->GetDigiIndex(0));
488  myTCand->AddHit(astripdigi->GetDetID(),ids.at(id),myHit->GetPosition().Mag());
489  }
490  }
491 
492  //mark used hits---------
493  hitsd.at(0).at(trackID1.at(i)).second=true;
494  hitsd.at(2).at(trackID3.at(i)).second=true;
495  for(Int_t id=0; id<otherIDs.size(); id++){
496  hitsd.at(otherIDs.at(id).first).at(otherIDs.at(id).second).second=true;
497  }
498 
499  // ///Add seed information to track------------
500  // PndSdsHit* myHit0 = (PndSdsHit*)(fStripHitArray->At(ids.at(0)));
501  // PndSdsHit* myHit1 = (PndSdsHit*)(fStripHitArray->At(ids.at(1)));
502  // TVector3 hit0 = myHit0->GetPosition(); TVector3 hit1 = myHit1->GetPosition();
503  // TVector3 posSeed(hit0.X(),hit0.Y(),hit0.Z());
504  // vec*=1./vec.Mag();
505  // //shift trk out of plane [needed for correct treatment in Kalman Fillter and GEANE]
506  // double sh_z = -0.035; //350 mkm
507  // double sh_x = vec.X()*sh_z;
508  // double sh_y = vec.Y()*sh_y;
509  // TVector3 sh_point(sh_x,sh_y,sh_z);
510  // posSeed +=sh_point;
511  // myTCand->setTrackSeed(posSeed,vec,-1);
512  // ///-------------------------------------
513 
514  tofill.push_back(*(myTCand)); //save Track Candidate
515  //new((*fTrackCandArray)[trackCnt]) PndTrackCand(*(myTCand));
516  delete myTCand;
517  }//Track Cand build
518  }//end of pseudo-tracks
519 
520 }
521 // -------------------------------------------------------------------------
522 
523 // ----- Private method FindHitsI --------------------------------------------
524 void PndStraightLineTrackFinderTask::FindHitsI(std::vector<PndTrackCand> &tofill, std::vector< std::vector< std::pair<Int_t,bool> > > &hitsd, Int_t nStripHits)
525 {
526 
527  std::vector<TVector3> trackStart, trackVec;//pseudo tracks
528  std::vector<TVector3> trackStartd, trackVecd;//save errors
529  std::vector< Int_t > trackID1, trackID2; //for TrackCand
530 
531 //iterate first discs-hits with all seconds, save pseudo-tracks
532  TVector3 start, tmp, vec, dstart, dvec; //temp-vars
533 
534  if(hitsd.size()<2) return;
535  for (Int_t i=0; i<hitsd.at(0).size(); i++)
536  {
537  PndSdsHit *hit1=(PndSdsHit*)fStripHitArray->At(hitsd.at(0).at(i).first);
538  start.SetXYZ(hit1->GetX(), hit1->GetY(), hit1->GetZ());
539  dstart.SetXYZ(hit1->GetDx(), hit1->GetDy(), hit1->GetDz());
540  for (Int_t k=0; k<hitsd.at(1).size(); k++)
541  {
542  PndSdsHit *hit2=(PndSdsHit*)fStripHitArray->At(hitsd.at(1).at(k).first);
543  tmp.SetXYZ(hit2->GetX(), hit2->GetY(), hit2->GetZ());
544  vec = tmp - start; //calc direction vector for FINDING
545  dvec.SetXYZ(hit2->GetDx(), hit2->GetDy(), hit2->GetDz());
546  // if(vec.Theta()>0.03 && vec.Theta()<0.05 && vec.Phi()>-0.25 && vec.Phi()<0.25){ //ignore vectors with theta outside 2-9 mrad
547  // if(vec.Theta()<0.01){ //ignore vectors with theta outside 10 mrad
548  trackStart.push_back(start);
549  trackStartd.push_back(dstart);
550  trackVec.push_back(vec); //save vector from start to second
551  trackVecd.push_back(dvec); //save error of second point for FIT, NOT error of direction vector
552  //trackID1.push_back(hitsd.at(0).at(i).first); //save Hit-Id's for TrackCand
553  //trackID2.push_back(hitsd.at(1).at(k).first);
554  trackID1.push_back(i); //save Hit-Id's for TrackCand
555  trackID2.push_back(k);
556  // }
557  }//end of disc 2 (or 3 if none in disc 2) hits
558  }//end of disc 1 hits
559 
560  if(fVerbose>1) cout << "Pseudos L1L2: "<< trackStart.size() <<endl;
561 
562  std::vector<Int_t> ids;
563 
564 //check if other discs have hits in track, add points for fitting ---------------
565  for (Int_t i=0; i<trackStart.size(); i++) //pseudo-loop
566  {
567  ids.clear();
568  Int_t pntcnt=2;
569 
570  start = trackStart.at(i);
571  dstart = trackStartd.at(i);
572  vec = trackVec.at(i);
573  dvec = trackVecd.at(i);
574  ids.push_back(hitsd.at(0).at(trackID1.at(i)).first);
575  ids.push_back(hitsd.at(1).at(trackID2.at(i)).first);
576  std::vector< std::pair<Int_t,Int_t> > otherIDs;
577 
578  for (Int_t idet=2; idet < 4; idet++){
579  Double_t distClosest = 2*idet*dXY; //just bigger as possible
580  Bool_t firstp = true;
581  for (Int_t ihit=0; ihit<hitsd.at(idet).size(); ihit++)
582  {
583  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(idet).at(ihit).first);
584  Double_t scale = (hit->GetZ()-start.z())/vec.z();
585  tmp = start + scale*vec; //extend search-vector to hit-plane
586  Double_t distTmp = sqrt((tmp.x()-hit->GetX())*(tmp.x()-hit->GetX()) + (tmp.y()-hit->GetY())*(tmp.y()-hit->GetY()));
587  if( distTmp<(idet*dXY) ){ //if in diameter
588  if(firstp){
589  pntcnt++;
590  ids.push_back(hitsd.at(idet).at(ihit).first);
591  otherIDs.push_back( make_pair (idet,ihit) );
592  distClosest=distTmp;
593  firstp=false;
594  }else{
595  if(distTmp<distClosest){ //change if nearer
596  ids.pop_back();
597  otherIDs.pop_back();
598  ids.push_back(hitsd.at(idet).at(ihit).first);
599  otherIDs.push_back( make_pair (idet,ihit) );
600  }
601  }
602  }//endif in diameter
603  }//end of disc n hits
604  }//end of discs
605 
606  if(fVerbose>2) cout << " Track L1L2: "<< i << "#Planes: " << ids.size() <<endl;
607 
608  //create track für fitting
609  if(ids.size()>2){ //third hit found <=> !track found!
610  // if(ids.size()>3){ //4 hits are needed because 6 parameters are used in trk-fit!
611  // // //third hit found <=> !track found!
612  // // //check direction for reduction of ghost tracks
613  // if(ids.size()>2 &&
614  // fabs(vec.Phi())<0.26 && vec.Theta()>3e-2 && vec.Theta()<5e-2){
615  if(fVerbose>2) cout << " Track: "<< i << "#Planes: " << ids.size() <<endl;
616  PndTrackCand *myTCand = new PndTrackCand();
617 
618  for (Int_t id=0; id<ids.size(); id++){
619  PndSdsHit* myHit = (PndSdsHit*)(fStripHitArray->At(ids.at(id)));
620  PndSdsClusterStrip* myCluster;
621  PndSdsDigiStrip* astripdigi;
622 
623  myTCand->AddHit(FairRootManager::Instance()->GetBranchId(fHitBranchStrip),ids.at(id),myHit->GetPosition().Mag());
624 
625  }
626 
627  //mark used hits---------
628  hitsd.at(0).at(trackID1.at(i)).second=true;
629  hitsd.at(1).at(trackID2.at(i)).second=true;
630  for(Int_t id=0; id<otherIDs.size(); id++){
631  hitsd.at(otherIDs.at(id).first).at(otherIDs.at(id).second).second=true;
632  }
633 
635  // Double_t Z0 = 1099.;
636  // // Double_t Z0 = 0;
637  // PndSdsHit* myHit0 = (PndSdsHit*)(fStripHitArray->At(ids.at(0)));
638  // PndSdsHit* myHit1 = (PndSdsHit*)(fStripHitArray->At(ids.at(1)));
639  // TVector3 hit0 = myHit0->GetPosition(); TVector3 hit1 = myHit1->GetPosition();
640  // double p1seed = (hit0.X()-hit1.X())/(hit0.Z()-hit1.Z());
641  // double p0seed = hit0.X() - p1seed*(hit0.Z()-Z0);
642  // // double p0seed = 0.5*(hit0.X()+hit1.X()-p1seed*(hit0.Z()+hit1.Z()-2*1099.)); //TO DO: don't use const
643  // double p3seed = (hit0.Y()-hit1.Y())/(hit0.Z()-hit1.Z());
644  // double p2seed = hit0.Y() - p3seed*(hit0.Z()-Z0);
645  // //double p2seed = 0.5*(hit0.Y()+hit1.Y()-p1seed*(hit0.Z()+hit1.Z()-2*1099.)); //TO DO: don't use const
646  // TVector3 posSeed(p0seed,p2seed,Z0);
647  // TVector3 dirSeed(p1seed,p3seed,1.);
648  // cout<<"posSeed:"<<endl;
649  // posSeed.Print();
650  // cout<<"dirSeed:"<<endl;
651  // dirSeed.Print();
652 
653  // // Double_t Z0 = 1099.;
654  // // TVector3 posSeed(start.X(),start.Y(),Z0);
655  // TVector3 posSeed(start.X(),start.Y(),start.Z());
656  // // cout<<"Trk cand: pos"<<endl;
657  // // posSeed.Print();
658  // vec*=1./vec.Mag();
659  // //shift trk out of plane [needed for correct treatment in Kalman Fillter and GEANE]
660  // double sh_z = -0.035; //350 mkm
661  // double sh_x = vec.X()*sh_z;
662  // double sh_y = vec.Y()*sh_y;
663  // TVector3 sh_point(sh_x,sh_y,sh_z);
664  // posSeed +=sh_point;
665 
666  // // cout<<"Trk cand: vec"<<endl;
667  // // vec.Print();
668  // myTCand->setTrackSeed(posSeed,vec,-1);
669  // ///-------------------------------------
670 
671 
672  tofill.push_back(*(myTCand)); //save Track Candidate
673  //new((*fTrackCandArray)[trackCnt]) PndTrackCand(*(myTCand));
674  delete myTCand;
675  }//Track Cand build
676  }//end of pseudo-tracks
677 
678 }
679 // -------------------------------------------------------------------------
680 
681 // ----- Public method Exec --------------------------------------------
682 void PndStraightLineTrackFinderTask::Exec(Option_t*)
683 {
684  TStopwatch *timer_exec = new TStopwatch();
685  if(fVerbose>2) timer_exec->Start();
686  if(fVerbose>2) cout << "Evt " << FairRootManager::Instance()->GetEntryNr() << " started--------------"<<endl<<endl;
687 
688  //which combinations of planes to build pseudo-vectos
689  //Int_t planeFlag=0; //flag=... 0: planes 1-2, 1: 1-2 & 1-3, 2: 1-2 & 1-3 & 2-3
690 
691  // Reset output array
692  if ( ! fTrackCandArray )
693  Fatal("Exec", "No trackCandArray");
694  fTrackCandArray->Clear();
695 
696  Int_t nStripHits = fStripHitArray->GetEntriesFast();
697  if(fVerbose>2) cout << "# Hits: \t"<< nStripHits <<endl<<endl;
698  bool usedFlag[nStripHits]; //for pseudo-vector building
699  for(Int_t seti=0; seti<nStripHits; seti++)
700  usedFlag[seti]=false;
701 
702  if(nStripHits<2){
703  if(fVerbose>2) cout << "Evt finsihed: too less hits-----"<<endl<<endl;
704  return;
705  }
706 
707  //-----do some sorting first----------
708  std::vector< std::vector< std::pair<int,bool> > > hitsd(4);
709  bool resSortHits;
710  if(flagStipSens) resSortHits = SortHitsByDet(hitsd, nStripHits);
711  else{
712  if(flagPixelSens) resSortHits = SortHitsByDet2(hitsd, nStripHits);
713  else{
714  std::cout<<"Algorithm is needed sensor type! Please, set it via SetSensStripFlag(bool fS) or SetSensPixelFlag(bool fS)"<<std::endl;
715  return;
716  }
717  }
718  if(!resSortHits){
719  if(fVerbose>2) cout << "Evt finsihed: too less planes-----"<<endl<<endl;
720  return;
721  }
722 
723 
724  if(fVerbose>2){
725  cout<<"HitMap size: "<< hitsd.size() <<endl;
726  if(hitsd.size()>0){
727  for (Int_t i=0; i<hitsd.at(0).size(); i++){
728  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(0).at(i).first);
729  cout<<"Plane0 Hit=("<<hit->GetX()<<", "<<hit->GetY()<<", "<<hit->GetZ()<<")"<<endl;
730  }
731  }
732  if(hitsd.size()>1){
733  for (Int_t i=0; i<hitsd.at(1).size(); i++){
734  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(1).at(i).first);
735  cout<<"Plane1 Hit=("<<hit->GetX()<<", "<<hit->GetY()<<", "<<hit->GetZ()<<")"<<endl;
736  }
737  }
738  if(hitsd.size()>2){
739  for (Int_t i=0; i<hitsd.at(2).size(); i++){
740  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(2).at(i).first);
741  cout<<"Plane2 Hit=("<<hit->GetX()<<", "<<hit->GetY()<<", "<<hit->GetZ()<<")"<<endl;
742  }
743  }
744  if(hitsd.size()>3){
745  for (Int_t i=0; i<hitsd.at(3).size(); i++){
746  PndSdsHit *hit=(PndSdsHit*)fStripHitArray->At(hitsd.at(3).at(i).first);
747  cout<<"Plane3 Hit=("<<hit->GetX()<<", "<<hit->GetY()<<", "<<hit->GetZ()<<")"<<endl;
748  }
749  }
750  }
751 
752  std::vector<PndTrackCand> theCands; //Track Candidates
753 
754  //--------find some tracks--------------
755  FindHitsI(theCands, hitsd, nStripHits);
756  if(fFinderMode){
757  FindHitsII(theCands, hitsd, nStripHits);
758  FindHitsIII(theCands, hitsd, nStripHits);
759  }
760 
761  //fill tracklist für fitting
762  for(int t=0; t<theCands.size(); t++)
763  new((*fTrackCandArray)[t]) PndTrackCand(theCands.at(t));
764 
765  if(fVerbose>2) cout << "Evt finsihed--------------"<<endl<<endl;
766  if(fVerbose>2){
767  timer_exec->Stop();
768  Double_t rtime_exec = timer_exec->RealTime();
769  Double_t ctime_exec = timer_exec->CpuTime();
770  cout << "Real time for Exec:" << rtime_exec << " s, CPU time " << ctime_exec << " s" << endl;
771  cout << endl;
772  }
773 }
774 // -------------------------------------------------------------------------
775 
776 Double_t PndStraightLineTrackFinderTask::GetTrackCurvature(PndMCTrack* myTrack)
777 {
778  TVector3 p = myTrack->GetMomentum();
779  return (2/TMath::Sqrt(p.Px()*p.Px() + p.Py()*p.Py()));
780 }
781 // -------------------------------------------------------------------------
782 
783 Double_t PndStraightLineTrackFinderTask::GetTrackDip(PndMCTrack* myTrack)
784 {
785  TVector3 p= myTrack->GetMomentum();
786  return (p.Mag()/TMath::Sqrt(p.Px()*p.Px() + p.Py()*p.Py()));
787 }
788 
789 
790 // -------------------------------------------------------------------------
791 ClassImp(PndStraightLineTrackFinderTask);
792 
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
p
Double_t p
Definition: anasim.C:58
fVerbose
int fVerbose
Definition: poormantracks.C:24
i
Int_t i
Definition: run_full.C:25
PndSdsHit::GetPosition
TVector3 GetPosition() const
Definition: PndSdsHit.h:93
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndSdsDigiStrip
Class for digitised strip hits.
Definition: PndSdsDigiStrip.h:29
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndStraightLineTrackFinderTask::FindHitsII
void FindHitsII(std::vector< PndTrackCand > &tofill, std::vector< std::vector< std::pair< Int_t, bool > > > &hitsd, Int_t nStripHits)
Definition: PndStraightLineTrackFinderTask.cxx:388
PndStraightLineTrackFinderTask::SortHitsByDet2
bool SortHitsByDet2(std::vector< std::vector< std::pair< Int_t, bool > > > &hitsd, Int_t nStripHits)
Definition: PndStraightLineTrackFinderTask.cxx:150
PndSdsCluster::GetDigiIndex
Int_t GetDigiIndex(Int_t i) const
Definition: PndSdsCluster.h:40
PndStraightLineTrackFinderTask::FindHitsIII
void FindHitsIII(std::vector< PndTrackCand > &tofill, std::vector< std::vector< std::pair< Int_t, bool > > > &hitsd, Int_t nStripHits)
Definition: PndStraightLineTrackFinderTask.cxx:252
PndTrackCand::AddHit
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Definition: PndTrackCand.cxx:36
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
z
Double_t z
Definition: createRootFscGeometryFile.C:146
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndStraightLineTrackFinderTask::SortHitsByDet
bool SortHitsByDet(std::vector< std::vector< std::pair< Int_t, bool > > > &hitsd, Int_t nStripHits)
Definition: PndStraightLineTrackFinderTask.cxx:119
PndStraightLineTrackFinderTask::GetTrackCurvature
Double_t GetTrackCurvature(PndMCTrack *myTrack)
Definition: PndStraightLineTrackFinderTask.cxx:776
PndSdsDigi::GetDetID
Int_t GetDetID() const
Definition: PndSdsDigi.h:61
ClassImp
ClassImp(PndAnaContFact)
PndStraightLineTrackFinderTask::FindHitsI
void FindHitsI(std::vector< PndTrackCand > &tofill, std::vector< std::vector< std::pair< Int_t, bool > > > &hitsd, Int_t nStripHits)
Definition: PndStraightLineTrackFinderTask.cxx:524
t
TTree * t
Definition: bump_analys.C:13
PndSdsHit::GetClusterIndex
Int_t GetClusterIndex() const
Definition: PndSdsHit.h:94
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
PndStraightLineTrackFinderTask::SortHitsByZ
bool SortHitsByZ(std::vector< std::vector< std::pair< Int_t, bool > > > &hitsd, Int_t nStripHits)
Definition: PndStraightLineTrackFinderTask.cxx:185
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
vec
dble_vec_t vec[12]
Definition: ranlxd.cxx:380
PndStraightLineTrackFinderTask::PndStraightLineTrackFinderTask
PndStraightLineTrackFinderTask(Int_t inFinderMode=0, TString hitBranch="LMDHitsStrip", TString clusterBranch="LMDStripClusterCand", TString digiBranch="LMDStripDigis", Int_t innSensPP=8)
Definition: PndStraightLineTrackFinderTask.cxx:22
Generated on Wed Apr 3 2019 08:44:01 for FairRoot/PandaRoot by   doxygen 1.8.5