FairRoot/PandaRoot
Public Member Functions | Private Member Functions | Private Attributes | List of all members
PndSttMvdGemTracking Class Reference

#include <PndSttMvdGemTracking.h>

Inheritance diagram for PndSttMvdGemTracking:
PndPersistencyTask

Public Member Functions

 PndSttMvdGemTracking ()
 
 PndSttMvdGemTracking (Int_t verbose)
 
 ~PndSttMvdGemTracking ()
 
virtual InitStatus Init ()
 
virtual void Exec (Option_t *opt)
 
void SetPersistenc (Bool_t persistence)
 
virtual void SetParContainers ()
 
void Copy (PndTrackCand *completeCand, PndTrack *completeTrack, PndTrackCand *sttmvdCand, PndTrack *sttmvd)
 
void SetTimes (Int_t times)
 
void SetMaxDistance (Double_t maxdistance)
 
void SwitchOnDisplay ()
 
void SetupGEMPlanes ()
 
void Reset (Int_t nhits, Int_t ntracks)
 
void OrderGemHits (Int_t nhits)
 
Int_t GetTrackIndex (Int_t i)
 
Int_t GetHitIndex (Int_t i)
 
Int_t GetPosIndex (PndGemHit *hit)
 
Int_t CountTracks ()
 
Int_t CountHitsInTrack (Int_t itrk)
 
void AddHitToTrack (Int_t ihit, Int_t itrk)
 
void DeleteHitFromTrack (Int_t ihit, Int_t itrk)
 
std::vector< int > GetTracksAssociatedToHit (Int_t ihit)
 
std::vector< int > GetHitsAssociatedToTrack (Int_t itrk)
 
std::vector< int > GetHitsAssociatedToTrackOnPlane (Int_t itrk, Int_t ipos)
 
void AddRemainingHits (Int_t ntracks)
 
void CheckCombinatorial (Int_t nhits, Int_t ntracks)
 
void ForbidMultiAssignedHits (Int_t nhits, Int_t ntracks)
 
void OnlyOneHitToEachTrack (Int_t nhits, Int_t ntracks)
 
void Retrack ()
 
Bool_t PropagateToGemPlane (FairTrackParP *tmppar, FairTrackParP *gempar, Int_t ipos)
 
Bool_t PropagateToGemPlaneAsHelix (PndTrack *sttmvd, FairTrackParP *gempar, Int_t ipos)
 
Double_t IsAssignable (FairTrackParP *gempar, PndGemHit *gemhit)
 
std::vector< int > AssignHits (Int_t itrk, FairTrackParP *gempar, Int_t ipos)
 
void EvaluatePerformances (Int_t nhits, Int_t ntracks)
 
void FillTrueDistances ()
 
Int_t SelectPdgCode (PndTrack *sttmvd)
 
void Kalman (TMatrixT< double > extrap, TMatrixT< double > measurement, TMatrixT< double > H, TMatrixT< double > extrap_cov, TMatrixT< double > measurement_cov, TMatrixT< double > &kalman, TMatrixT< double > &kalman_cov)
 
FairTrackParP SetStartParameters (PndTrack *sttmvd, PndTrackCand *sttmvdCand)
 
Int_t GetClosestOnFirst (FairTrackParP *gempar, Int_t ipos, Double_t &closestdistance)
 
Bool_t Prefit (PndTrack *sttmvdTrack, PndTrackCand *sttmvdCand, TVector3 &lastpos, TVector3 &lastmom)
 
Bool_t IntersectionFinder (Double_t xc, Double_t yc, Double_t radius, PndSttHit *stthit, TVector3 &xyz, TVector3 &dxyz)
 
Bool_t Fit (TMatrixT< double > points, Double_t &outxc, Double_t &outyc, Double_t &outradius)
 
Bool_t ZFit (TMatrixT< double > points, Int_t charge, Double_t xc, Double_t yc, Double_t radius, Double_t &fitm, Double_t &fip)
 
Bool_t GetInitialParams (PndTrack *sttmvd, Double_t &xc, Double_t &yc, Double_t &radius, Double_t &fitm, Double_t &fitp)
 
Double_t CalculatePhi (TVector2 v, TVector2 p, double alpha, double Phi0, int charge)
 
Double_t CompareToPreviousPhi (Double_t Fi, Double_t Fi_pre, int charge)
 
Bool_t ZFind (Int_t nhits, TMatrixT< double > points, Double_t xc, Double_t yc, Double_t radius)
 
void ConsiderCombinatorialEffect (Int_t nhits)
 
void SetCombinatorialDistance (Double_t combidistance)
 
void UpdateMCTrackId (PndTrackCand *completeCand)
 
void UseMonteCarlo ()
 
void WriteHistograms ()
 
void SetEvaluateFlag (Bool_t flag)
 
Int_t GetPdgFromMC (int trackid)
 
Int_t GetChargeCorrectedPdgFromMC (int trackid, int charge)
 
void SetPdgFromMC ()
 
void SetPdgFromMC (TString trackid)
 
void SetDefaultPdg (int pdg)
 
void SetBranchNames (TString mvdpixel, TString mvdstrip, TString stt, TString gem)
 
void SetTrackBranchNames (TString startingtrack, TString startingtrackcand)
 
void SetPersistency (Bool_t val=kTRUE)
 
Bool_t GetPersistency ()
 

Private Member Functions

 ClassDef (PndSttMvdGemTracking, 1)
 

Private Attributes

TClonesArray * fGemHitArray
 
TClonesArray * fGemPointArray
 
TClonesArray * fMvdPixelHitArray
 
TClonesArray * fMvdStripHitArray
 
TClonesArray * fSttHitArray
 
TClonesArray * fMvdPointArray
 
TClonesArray * fSttPointArray
 
TClonesArray * fMCTrackArray
 
TClonesArray * fTrackArray
 
TClonesArray * fTrackIDArray
 
TClonesArray * fTrackCandArray
 
TClonesArray * fCompleteTrackCandArray
 
TClonesArray * fCompleteTrackArray
 
TClonesArray * fTubeArray
 
Bool_t fPdgFromMC
 
PndGeoSttParfSttParameters
 
PndGemDigiParfGemParameters
 
FairGeanePro * fPro
 
Int_t countgood
 
Int_t countbad
 
Int_t countdoubt
 
Int_t countreconstructablehit
 
Int_t countprinotassigned
 
Int_t countsecnotassigned
 
Int_t countreconstructablepoint
 
Int_t countplane [8]
 
Int_t countprop [2][8]
 
Int_t countsttmvd
 
Int_t countsttmvdusable
 
Int_t countnohitonplane [8]
 
Int_t evt
 
std::vector< int > usabletracks
 
std::vector< int > fOrdering
 
std::vector< int >::iterator fOrderingIterator
 
TClonesArray * fSensPositions
 
Int_t fNPositions
 
TVector3 fDj
 
TVector3 fDk
 
Int_t fPdgCode
 
Int_t fDefaultPdgCode
 
std::vector< std::pair< int,
int > > 
trackvector
 
std::vector< int > trackindexes
 
std::vector< int > notassignedhits
 
std::vector< int > notassignedtracks
 
std::map< int, bool > towhichplane
 
std::map< int, bool > fProTracks
 
TMatrixT< float > hitmap
 
TMatrixT< float > hitcounter
 
TMatrixT< double > distancemap
 
Double_t fMaxDistance
 
Int_t fTimes
 
Int_t fTurn
 
Bool_t fDisplayOn
 
TCanvas * display
 
TH2F * h [8]
 
TH2F * hnotskewed
 
TH2F * hskewed
 
TH1F * hdist [8]
 
TH1F * hdist2 [8]
 
TH1F * hsigma [8]
 
TH1F * hsigma2 [8]
 
TH1F * hchosen [8]
 
TH1F * hchosen2 [8]
 
TH1F * hmcdist [8]
 
TH1F * hmcx [8]
 
TH1F * hmcy [8]
 
Bool_t fUseMC
 
TString fMvdPixelBranchName
 
TString fMvdStripBranchName
 
TString fSttBranchName
 
TString fGemBranchName
 
TString fStartTrackBranchName
 
TString fStartTrackCandBranchName
 
TString fStartTrackIDBranchName
 
std::map< int, int > fCombiMap
 
Double_t fCombiDistance
 
Bool_t fEvaluate
 

Detailed Description

Definition at line 18 of file PndSttMvdGemTracking.h.

Constructor & Destructor Documentation

PndSttMvdGemTracking::PndSttMvdGemTracking ( )

Default constructor

Definition at line 60 of file PndSttMvdGemTracking.cxx.

References fCombiDistance, fDefaultPdgCode, fDisplayOn, fEvaluate, fGemBranchName, fMaxDistance, fMvdPixelBranchName, fMvdStripBranchName, fPdgFromMC, fStartTrackBranchName, fStartTrackCandBranchName, fStartTrackIDBranchName, fSttBranchName, fTimes, fTurn, fUseMC, fVerbose, and PndPersistencyTask::SetPersistency().

60  :
61  PndPersistencyTask("MVD-STT-GEM tracking") {
62  fVerbose = 1;
63  fEvaluate = kTRUE;
64  fDisplayOn = false;
65  fTimes = 5;
66  fTurn = 1;
67  fMaxDistance = -1;
68  fUseMC = kFALSE;
69  fCombiDistance = 1.;
70  fDefaultPdgCode = -13;
71  fPdgFromMC = kFALSE;
72 
73  fMvdPixelBranchName = "MVDHitsPixel";
74  fMvdStripBranchName = "MVDHitsStrip";
75  fSttBranchName = "STTHit";
76  fGemBranchName = "GEMHit";
77 
78  fStartTrackBranchName = "SttMvdTrack";
79  fStartTrackCandBranchName = "SttMvdTrackCand";
80  fStartTrackIDBranchName = "SttMvdTrackID";
81 
82  SetPersistency(kTRUE);
83 
84 }
int fVerbose
Definition: poormantracks.C:24
void SetPersistency(Bool_t val=kTRUE)
PndSttMvdGemTracking::PndSttMvdGemTracking ( Int_t  verbose)

Definition at line 87 of file PndSttMvdGemTracking.cxx.

References fCombiDistance, fDefaultPdgCode, fDisplayOn, fEvaluate, fGemBranchName, fMaxDistance, fMvdPixelBranchName, fMvdStripBranchName, fPdgFromMC, fStartTrackBranchName, fStartTrackCandBranchName, fStartTrackIDBranchName, fSttBranchName, fTimes, fTurn, fUseMC, fVerbose, PndPersistencyTask::SetPersistency(), and verbose.

87  :
88  PndPersistencyTask("MVD-STT-GEM tracking") {
89  SetPersistency(kTRUE);
90  fVerbose = verbose;
91  if(verbose > 0) fEvaluate = kTRUE;
92  else fEvaluate = kFALSE;
93  fDisplayOn = false;
94  fTimes = 5;
95  fTurn = 1;
96  fMaxDistance = -1;
97  fUseMC = kFALSE;
98  fCombiDistance = 1.;
99  fDefaultPdgCode = -13;
100  fPdgFromMC = kFALSE;
101 
102  fMvdPixelBranchName = "MVDHitsPixel";
103  fMvdStripBranchName = "MVDHitsStrip";
104  fSttBranchName = "STTHit";
105  fGemBranchName = "GEMHit";
106 
107  fStartTrackBranchName = "SttMvdTrack";
108  fStartTrackCandBranchName = "SttMvdTrackCand";
109  fStartTrackIDBranchName = "SttMvdTrackID";
110 
111 }
int fVerbose
Definition: poormantracks.C:24
#define verbose
void SetPersistency(Bool_t val=kTRUE)
PndSttMvdGemTracking::~PndSttMvdGemTracking ( )

Destructor

Definition at line 116 of file PndSttMvdGemTracking.cxx.

116  {
117 }

Member Function Documentation

void PndSttMvdGemTracking::AddHitToTrack ( Int_t  ihit,
Int_t  itrk 
)

Definition at line 1292 of file PndSttMvdGemTracking.cxx.

References fVerbose, notassignedhits, and trackvector.

Referenced by AddRemainingHits(), Exec(), and Retrack().

1292  {
1293  std::pair<int, int> thispair(itrk, ihit);
1294  trackvector.push_back(thispair);
1295  std::vector<int>::iterator iter;
1296  iter = std::find(notassignedhits.begin(), notassignedhits.end(), ihit);
1297  int where = iter - notassignedhits.begin();
1298  if(where != (int)notassignedhits.size()) notassignedhits.erase(iter); // remove from not assigned list
1299  if(fVerbose > 0) cout << "ADD HIT " << ihit << " TO TRK " << itrk << endl;
1300 }
int fVerbose
Definition: poormantracks.C:24
std::vector< std::pair< int, int > > trackvector
std::vector< int > notassignedhits
void PndSttMvdGemTracking::AddRemainingHits ( Int_t  ntracks)

Definition at line 1400 of file PndSttMvdGemTracking.cxx.

References AddHitToTrack(), CHECKCOMBI, CountTracks(), distancemap, Double_t, fCombiMap, fGemHitArray, fNPositions, fOrdering, fOrderingIterator, fProTracks, fVerbose, GetHitsAssociatedToTrack(), GetPosIndex(), GetTrackIndex(), GetTracksAssociatedToHit(), hitcounter, hitmap, i, and towhichplane.

Referenced by Exec().

1400  { // ntracks //[R.K.03/2017] unused variable(s)
1401 
1402  if(fVerbose > 0) cout << "ADD REMAINING HITS to " << CountTracks() << " tracks " << endl;
1403 
1404  for(Int_t k = 0; k < CountTracks(); k++) {
1405  int itrk = GetTrackIndex(k);
1406  if(fProTracks[itrk] == false) continue;
1407  if(fVerbose > 0) cout << k << " ITRK " << itrk << endl;
1408 
1409  std::vector<int> hitvector = GetHitsAssociatedToTrack(itrk);
1410  Int_t nhitinthistrack = hitvector.size();
1411 
1412  // if all the positions are filled then continue
1413  if(nhitinthistrack == fNPositions) { if(fVerbose > 0) cout << "all pos filled" << endl; continue; }
1414 
1415  // create a vector of all the positions and
1416  // initialize it as all full
1417  std::vector<int> missingpositions;
1418  for(int ipos = 0; ipos < fNPositions; ipos++) missingpositions.push_back(ipos);
1419 
1420  // run over the track hits and delete from
1421  // the list the positions already filled
1422  std::vector<int>::iterator iter;
1423  for(size_t i = 0; i < hitvector.size(); i++) {
1424  PndGemHit *gemhit =(PndGemHit*) fGemHitArray->At(hitvector[i]);
1425  if(!gemhit) continue;
1426 
1427  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), GetPosIndex(gemhit));
1428  int ipos = fOrderingIterator - fOrdering.begin();
1429  iter = std::find(missingpositions.begin(), missingpositions.end(), ipos);
1430  Int_t where = iter - missingpositions.begin();
1431  if(where != (int)missingpositions.size()) {
1432  if(fVerbose > 0) cout << "here" << endl; // CHECK delete it
1433  missingpositions.erase(iter);
1434  }
1435  }
1436 
1437  // for each missing pos check
1438  // whether there are hits there
1439  // or not
1440  for(size_t i = 0; i < missingpositions.size(); i++) {
1441  Int_t ipos = missingpositions[i];
1442  if(fVerbose > 0) cout << "missing position " << ipos << endl;
1443  // if there are not hits continue
1444  if(towhichplane[ipos] == false) continue;
1445 
1446  // if there are, run over them and pick the closest
1447  // if it is below limit distance
1448  Double_t tmpdist = 1000;
1449  Int_t tmphitindex = -1;
1450  for(Int_t ihit = 0; ihit < hitcounter(ipos, 0); ihit++) {
1451  int hitindex = (int) hitmap(ipos, ihit);
1452  if(GetTracksAssociatedToHit(hitindex).size() != 0) continue;
1453  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(hitindex);
1454  if(!gemhit) continue;
1455  if(CHECKCOMBI > 0) if(fCombiMap[hitindex] != 0) continue;
1456  if(fVerbose > 0) cout << "distance " << distancemap[itrk][hitindex] << endl;
1457  if(distancemap[itrk][hitindex] == -1) continue;
1458  if(distancemap[itrk][hitindex] < tmpdist) {
1459  tmpdist = distancemap[itrk][hitindex];
1460  tmphitindex = hitindex;
1461  }
1462  }
1463  if(tmphitindex != -1) {
1464  AddHitToTrack(tmphitindex, itrk);
1465  if(fVerbose > 0) cout << "ADDED " << tmphitindex << " TO TRK " << itrk << " at distance " << tmpdist << endl;
1466  }
1467  }
1468  }
1469 
1470 }
std::vector< int >::iterator fOrderingIterator
int fVerbose
Definition: poormantracks.C:24
TMatrixT< float > hitmap
Int_t GetPosIndex(PndGemHit *hit)
Int_t i
Definition: run_full.C:25
std::vector< int > GetHitsAssociatedToTrack(Int_t itrk)
std::vector< int > fOrdering
Double_t
TMatrixT< double > distancemap
TMatrixT< float > hitcounter
std::map< int, bool > fProTracks
#define CHECKCOMBI
void AddHitToTrack(Int_t ihit, Int_t itrk)
std::vector< int > GetTracksAssociatedToHit(Int_t ihit)
std::map< int, bool > towhichplane
std::map< int, int > fCombiMap
std::vector< int > PndSttMvdGemTracking::AssignHits ( Int_t  itrk,
FairTrackParP *  gempar,
Int_t  ipos 
)

Definition at line 1777 of file PndSttMvdGemTracking.cxx.

References CHECKCOMBI, display, distancemap, Double_t, fCombiMap, fDisplayOn, fGemHitArray, fVerbose, hitcounter, hitmap, and IsAssignable().

Referenced by Exec(), and Retrack().

1777  {
1778  std::vector<int> assignedhits; // CHECK need init?
1779  // compare the position with each *******************************
1780  // GEM hit on THIS sensor
1781  int hitonsensor = (int) hitcounter(ipos, 0);
1782  for(int ihit = 0; ihit < hitonsensor; ihit++) {
1783  int hitindex = (int) hitmap(ipos, ihit);
1784  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(hitindex);
1785  if(!gemhit) continue;
1786  if(CHECKCOMBI == 2 && fCombiMap[hitindex] != 0) continue;
1787  Double_t distance = IsAssignable(gempar, gemhit);
1788  // if the track to hit distance is below threshold
1789  // "distance" is filled (-1 otherwise)
1790  if(distance >= 0.) {
1791  // the separate coordinates are not working at all ****
1792  // if((distVec.X()/ TMath::Sqrt(hitErr.X() * hitErr.X() + gemErr.X() * gemErr.X())) <= 5 && (distVec.Y()/ TMath::Sqrt(hitErr.Y() * hitErr.Y() + gemErr.Y() * gemErr.Y())) <= 5) { // CHECK delete this and do not use it!
1793  // COUNTERS 167 0 0 210 43 43 0
1794  // COUNTERS 740 627 16 1101 -266 234 2
1795  // COUNTERS 241 385 8 359 -267 64 0
1796 
1797  // COUNTERS 189 0 0 215 26 26 4
1798  // COUNTERS 729 318 32 770 -277 87 3
1799  // COUNTERS 87 66 0 91 -62 1 0
1800 
1801  // COUNTERS 205 1 0 211 5 6 7
1802  // COUNTERS 806 300 13 838 -268 53 14
1803  // COUNTERS 1209 1087 0 1281 -1015 31 2
1804  // ****************************************************
1805 
1806  if(fVerbose > 0) cout << "this hit " << hitindex << " BELONGS to this track " << itrk << "!!!!!!!!!!!!!" << endl;
1807  // fill the distancemap: for each track compute its
1808  // distance from each hit and put it in the map
1809  distancemap[itrk][hitindex] = distance;
1810 
1811  // assign it to all the tracks it might belong to (for now)
1812  // completeCand->AddHit(FairRootManager::Instance()->GetBranchId(fGemBranchName), hitindex, gemhit->GetPosition().Mag()); // CHECK rho and kGemHit
1813  if(fVerbose != 0) cout << "assign " << hitindex << " to track " << itrk << endl;
1814  // AddHitToTrack(hitindex, itrk);
1815  assignedhits.push_back(hitindex);
1816 
1817  if(fDisplayOn) {
1818  TMarker *g1mrk = new TMarker(gemhit->GetX(), gemhit->GetY(), 7);
1819  g1mrk->SetMarkerColor(4); // blue
1820  display->cd(1 + ipos);
1821  g1mrk->Draw("SAME");
1822  display->Update();
1823  display->Modified();
1824  }
1825  }
1826  else if(fVerbose > 0) cout << "this hit " << hitindex << " DOES NOT BELONG to this track " << itrk << "!!!!!!!!!!!!!" << endl;
1827  }
1828 
1829  return assignedhits;
1830 }
int fVerbose
Definition: poormantracks.C:24
TMatrixT< float > hitmap
Double_t IsAssignable(FairTrackParP *gempar, PndGemHit *gemhit)
Double_t
TMatrixT< double > distancemap
TMatrixT< float > hitcounter
#define CHECKCOMBI
std::map< int, int > fCombiMap
Double_t PndSttMvdGemTracking::CalculatePhi ( TVector2  v,
TVector2  p,
double  alpha,
double  Phi0,
int  charge 
)

Definition at line 2914 of file PndSttMvdGemTracking.cxx.

References Double_t, Pi, and pi.

Referenced by GetInitialParams(), Prefit(), PropagateToGemPlaneAsHelix(), and ZFit().

2915 {
2916  Double_t Fi = - charge * TMath::ACos(v * p / (v.Mod() * p.Mod()));
2917  double pi = TMath::Pi();
2918  double pi2 = 2 * pi;
2919 
2920  // Fi = h * (pi2 - h * Fi) // should be correct
2921  if((charge > 0 && ((Phi0 > 0 && ((alpha > 0 && alpha > Phi0) ||
2922  (alpha < 0 && alpha < Phi0 - pi)))
2923  ||
2924  ((Phi0 < 0 && ((alpha > 0 && alpha < pi + Phi0) ||
2925  (alpha < 0 && alpha > Phi0)))) ))) Fi = - (pi2 + Fi) ;
2926  else if((charge < 0 && ((Phi0 > 0 && ((alpha > 0 && alpha < Phi0) ||
2927  (alpha < 0 && alpha > Phi0 - pi)))
2928  ||
2929  ((Phi0 < 0 && ((alpha > 0 && alpha > pi + Phi0) ||
2930  (alpha < 0 && alpha < Phi0)))) ))) Fi = pi2 - Fi ;
2931 
2932  return Fi;
2933 }
Double_t p
Definition: anasim.C:58
#define pi
Definition: createSTT.C:60
__m128 v
Definition: P4_F32vec4.h:4
Double_t
double alpha
Definition: f_Init.h:9
Double_t Pi
void PndSttMvdGemTracking::CheckCombinatorial ( Int_t  nhits,
Int_t  ntracks 
)

Definition at line 3207 of file PndSttMvdGemTracking.cxx.

References combi, count, DeleteHitFromTrack(), fCombiMap, fGemHitArray, fNPositions, fOrdering, fOrderingIterator, fProTracks, fVerbose, GetHitsAssociatedToTrack(), GetPosIndex(), GetTrackIndex(), and i.

Referenced by Exec().

3208 {
3209 
3210  if(nhits == 0) return;
3211 
3212  if(fVerbose > 0) cout << "CHECK COMBINATORIAL: DELETE FAKE HITS" << endl;
3213 
3214  // loop over the tracks
3215  for(Int_t it = 0; it < ntracks; it++) {
3216  int itrk = GetTrackIndex(it);
3217  if(fProTracks[itrk] == false) continue;
3218  std::vector<int> thistrackhits = GetHitsAssociatedToTrack(itrk);
3219  if(thistrackhits.size() == 0) continue;
3220  // fill table:
3221  // each row is a sensor plane and contains
3222  // the hits associtaed to this track on it
3223  TMatrixT<double> combi(fNPositions, thistrackhits.size());
3224  TMatrixT<double> addhit(fNPositions, 1);
3225 
3226  // init
3227  for(Int_t ipos = 0; ipos < fNPositions; ipos++) {
3228  addhit[ipos][0] = 0;
3229  for(size_t j = 0; j < thistrackhits.size(); j++) combi[ipos][j] = -1;
3230  }
3231 
3232  for(size_t j = 0; j < thistrackhits.size(); j++) {
3233  int ihit = thistrackhits[j];
3234 
3235  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(ihit);
3236  if(!gemhit) continue;
3237 
3238  int posindex = GetPosIndex(gemhit);
3239  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
3240  int index = fOrderingIterator - fOrdering.begin();
3241 
3242  int hitpos = (int) addhit[index][0];
3243  combi[index][hitpos] = ihit;
3244  addhit[index][0]++;
3245  }
3246 
3247  if(fVerbose > 0) {
3248  cout << "itrk " << itrk << " " << nhits << " " << thistrackhits.size() << endl;
3249  combi.Print();
3250  }
3251 
3252 
3253  // loop on sensor planes for this track
3254  for(Int_t ipos = 0; ipos < fNPositions; ipos++) {
3255 
3256  int count = 0;
3257  // count hits associated on this sensor plane
3258  for(size_t i = 0; i < thistrackhits.size(); i++) {
3259  if(combi[ipos][i] != -1) count++;
3260  }
3261 
3262  if(fVerbose > 0) cout << count << " hits on pos " << ipos << endl;
3263 
3264  // if there is no hit or only one, continue ...
3265  if(count <= 1) continue;
3266  // ... else
3267  int count2 = 0;
3268  // count how many true (non combi) hits are there
3269  for(size_t i = 0; i < thistrackhits.size(); i++) {
3270  int ihit = (int) combi[ipos][i];
3271  if(ihit == -1) continue;
3272  if(fCombiMap[ihit] == 0) count2++;
3273  }
3274 
3275  if(fVerbose > 0) cout << count2 << " true hits on pos " << ipos << endl;
3276 
3277  // if there is no true hit, continue ...
3278  if(count2 == 0) continue;
3279  // ... else, clean up from combinatorial hits
3280  for(size_t i = 0; i < thistrackhits.size(); i++) {
3281  int ihit = (int) combi[ipos][i];
3282  if(fCombiMap[ihit] != 0) {
3283  if(fVerbose > 0) cout << "delete " << ihit << " from track " << itrk << endl;
3284  DeleteHitFromTrack(ihit, itrk);
3285  }
3286  }
3287  }
3288  }
3289 }
std::vector< int >::iterator fOrderingIterator
int fVerbose
Definition: poormantracks.C:24
Int_t GetPosIndex(PndGemHit *hit)
Int_t i
Definition: run_full.C:25
std::vector< int > GetHitsAssociatedToTrack(Int_t itrk)
std::vector< int > fOrdering
TGeoCombiTrans * combi
std::map< int, bool > fProTracks
void DeleteHitFromTrack(Int_t ihit, Int_t itrk)
int count
std::map< int, int > fCombiMap
PndSttMvdGemTracking::ClassDef ( PndSttMvdGemTracking  ,
 
)
private
Double_t PndSttMvdGemTracking::CompareToPreviousPhi ( Double_t  Fi,
Double_t  Fi_pre,
int  charge 
)

Definition at line 2935 of file PndSttMvdGemTracking.cxx.

References Pi, and pi.

Referenced by GetInitialParams(), Prefit(), and ZFit().

2936 {
2937  // if(fabs(Fi) < fabs(Fi_pre)) Fi += h * pi2 // CHECK should be ok
2938  double pi = TMath::Pi();
2939  double pi2 = 2 * pi;
2940 
2941  if(charge < 0 && Fi < Fi_pre) Fi += pi2;
2942  else if(charge > 0 && Fi > Fi_pre) Fi -= pi2;
2943  Fi_pre = Fi;
2944  return Fi;
2945 }
#define pi
Definition: createSTT.C:60
Double_t Pi
void PndSttMvdGemTracking::ConsiderCombinatorialEffect ( Int_t  nhits)

sensor[ihit][3] = posindex; switch(posindex) { case 11: sensor[ihit][3] = 0; break; case 12: sensor[ihit][3] = 1; break; case 21: sensor[ihit][3] = 2; break; case 22: sensor[ihit][3] = 3; break; case 31: sensor[ihit][3] = 4; break; case 32: sensor[ihit][3] = 5; break; }

Definition at line 2948 of file PndSttMvdGemTracking.cxx.

References counter1, counter2, display, fCombiDistance, fCombiMap, fDisplayOn, fGemHitArray, fNPositions, fOrdering, fOrderingIterator, PndGemHit::GetDigiNr(), GetPosIndex(), PndGemHit::GetSensorNr(), PndGemHit::GetStationNr(), hit, nhits, sensor, and CAMath::Sqrt().

Referenced by Exec().

2948  {
2949 
2950  if(nhits == 0) return;
2951 
2952  // matrix hitid x y posindex = istat * 10 + isens iflag
2953  TMatrixT<double> sensor(nhits, 5);
2954  TMatrixT<double> nhitsonsensor(fNPositions, 1);
2955  TMatrixT<double> nhitsonsensor2(fNPositions, nhits);
2956 
2957 // cout << "gem n hits " << nhits << endl;
2958 // std::vector<int> mcpoints[fNPositions]; // CHECK
2959  std::vector< std::vector<int> > mcpoints; // CHECK
2960 
2961  for (Int_t ihit = 0; ihit < nhits; ihit++) {
2962  PndGemHit *hit = (PndGemHit*) fGemHitArray->At(ihit);
2963  if(!hit) continue;
2964  // cout << "ihit " << ihit << endl;
2965  sensor[ihit][0] = ihit;
2966  sensor[ihit][1] = hit->GetX();
2967  sensor[ihit][2] = hit->GetY();
2968  //int istat = hit->GetStationNr(); //[R.K. 01/2017] unused variable?
2969  //int isens = hit->GetSensorNr(); //[R.K. 01/2017] unused variable?
2970  int posindex = GetPosIndex(hit);
2971  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
2972  int ipos = fOrderingIterator - fOrdering.begin();
2973  sensor[ihit][3] = ipos;
2974  sensor[ihit][4] = true;
2993  nhitsonsensor2[ipos][(int) nhitsonsensor[ipos][0]] = ihit;
2994  nhitsonsensor[ipos][0]++;
2995  }
2996 
2997 
2998  // std::vector<int> accepted[fNPositions];
2999  std::vector< std::vector<int> > accepted;
3000  for(int ipos = 0; ipos < fNPositions; ipos++) {
3001  std::vector<int> acc;
3002  accepted.push_back(acc);
3003  }
3004 
3005  // nhitsonsensor2.Print();
3006 
3007 // cout << "SENS I, STAT I " << nhitsonsensor[0][0] << endl;
3008 // cout << "SENS 2, STAT I " << nhitsonsensor[1][0] << endl;
3009 // cout << "SENS I, STAT II " << nhitsonsensor[2][0] << endl;
3010 // cout << "SENS 2, STAT II " << nhitsonsensor[3][0] << endl;
3011 // cout << "SENS I, STAT III " << nhitsonsensor[4][0] << endl;
3012 // cout << "SENS 2, STAT III " << nhitsonsensor[5][0] << endl;
3013 
3014 
3015  int counter1 = 0;
3016 // cout << "nhits " << nhits << endl;
3017  for(int ihit = 0; ihit < nhits; ihit++) {
3018 
3019  if(sensor[ihit][3] != 0 &&
3020  sensor[ihit][3] != 2 &&
3021  sensor[ihit][3] != 4) continue;
3022  int first = (int) sensor[ihit][3];
3023  if(sensor[ihit][4] == false) continue;
3024  // cout << "FIRST " << first << endl;
3025  double x1 = sensor[ihit][1];
3026  double y1 = sensor[ihit][2];
3027  counter1++;
3028 
3029  int counter2 = 0;
3030  for(int jhit = 0; jhit < nhits; jhit++) {
3031 
3032  if(sensor[jhit][3] != first + 1) continue;
3033  int second = (int) sensor[jhit][3];
3034  if(sensor[jhit][4] == false) continue;
3035  // cout << "SECOND " << second << endl;
3036 
3037  double x2 = sensor[jhit][1];
3038  double y2 = sensor[jhit][2];
3039  counter2++;
3040  double distance = TMath::Sqrt((x1 - x2) * (x1 - x2) +
3041  (y1 - y2) * (y1 - y2));
3042 
3043  if(distance < fCombiDistance) {
3044  std::vector< int > acc1 = accepted[first];
3045  std::vector< int > acc2 = accepted[second];
3046  bool alreadythere1 = false, alreadythere2 = false;
3047  for(size_t j = 0; j < acc1.size(); j++) {
3048  if(sensor[ihit][0] == acc1[j]) {
3049  alreadythere1 = true;
3050  break;
3051  }
3052  }
3053  for(size_t j = 0; j < acc2.size(); j++) {
3054  if(sensor[jhit][0] == acc2[j]) {
3055  alreadythere2 = true;
3056  break;
3057  }
3058  }
3059 
3060  if(alreadythere1 == false) {
3061  // accepted[first].push_back((int) sensor[ihit][0]);
3062  acc1.push_back((int) sensor[ihit][0]);
3063  accepted[first] = acc1;
3064  fCombiMap[ihit] = 0;
3065 
3066  // ^^^^^^^
3067  PndGemHit *hit = (PndGemHit*) fGemHitArray->At(ihit);
3068  if(!hit) continue;
3069  double digis[2];
3070  digis[0] = hit->GetDigiNr(0);
3071  digis[1] = hit->GetDigiNr(1);
3072  int stat = hit->GetStationNr();
3073  int sens = hit->GetSensorNr();
3074  int posindex = stat * 10 + sens;
3075  switch(posindex) {
3076  case 11:
3077  posindex = 0; break;
3078  case 12:
3079  posindex = 1; break;
3080  case 21:
3081  posindex = 2; break;
3082  case 22:
3083  posindex = 3; break;
3084  case 31:
3085  posindex = 4; break;
3086  case 32:
3087  posindex = 5; break;
3088  }
3089 
3090 
3091  // cout << "SETTING " << ihit << " to false with ";
3092  for(int khit = 0; khit < nhits; khit++) {
3093 
3094  if(ihit == khit) continue;
3095  if(sensor[khit][3] != posindex) continue;
3096  // if(posindex == 1) cout << "? " << khit << endl;
3097  PndGemHit *ghit = (PndGemHit*) fGemHitArray->At(khit);
3098  if(!ghit) continue;
3099 
3100  if(ghit->GetDigiNr(0) == digis[0] || ghit->GetDigiNr(0) == digis[1]) {
3101  // cout << khit << " ";
3102  sensor[khit][4] = false;
3103  }
3104  if(ghit->GetDigiNr(1) == digis[0] || ghit->GetDigiNr(1) == digis[1]) {
3105  // cout << khit << " ";
3106  sensor[khit][4] = false;
3107  }
3108  }
3109  sensor[ihit][4] = false;
3110 
3111  // cout << endl;
3112  // ^^^^^^^
3113 
3114  if(fDisplayOn == kTRUE) {
3115  TMarker *amrk = new TMarker(x1, y1, 21);
3116  amrk->SetMarkerColor(5);
3117  display->cd(first + 1);
3118  amrk->Draw("SAME");
3119  display->Update();
3120  display->Modified();
3121  }
3122  }
3123 
3124  if(alreadythere2 == false) {
3125 // accepted[second].push_back((int) sensor[jhit][0]);
3126  acc2.push_back((int) sensor[jhit][0]);
3127  accepted[second] = acc2;
3128 
3129  fCombiMap[jhit] = 0;
3130 
3131  // ^^^^^^^
3132  PndGemHit *hit = (PndGemHit*) fGemHitArray->At(jhit);
3133  if(!hit) continue;
3134  double digis[2];
3135  digis[0] = hit->GetDigiNr(0);
3136  digis[1] = hit->GetDigiNr(1);
3137  int stat = hit->GetStationNr();
3138  int sens = hit->GetSensorNr();
3139  int posindex = stat * 10 + sens;
3140  switch(posindex) {
3141  case 11:
3142  posindex = 0; break;
3143  case 12:
3144  posindex = 1; break;
3145  case 21:
3146  posindex = 2; break;
3147  case 22:
3148  posindex = 3; break;
3149  case 31:
3150  posindex = 4; break;
3151  case 32:
3152  posindex = 5; break;
3153  }
3154 
3155  // cout << "SETTING " << jhit << " to false with ";
3156  for(int khit = 0; khit < nhits; khit++) {
3157  if(jhit == khit) continue;
3158  if(sensor[khit][3] != posindex) continue;
3159 
3160  PndGemHit *ghit = (PndGemHit*) fGemHitArray->At(khit);
3161  if(!ghit) continue;
3162  if(ghit->GetDigiNr(0) == digis[0] || ghit->GetDigiNr(0) == digis[1]) {
3163  // cout << khit << " ";
3164  sensor[khit][4] = false;
3165  }
3166  if(ghit->GetDigiNr(1) == digis[0] || ghit->GetDigiNr(1) == digis[1]) {
3167  // cout << khit << " ";
3168  sensor[khit][4] = false;
3169  }
3170  }
3171  sensor[jhit][4] = false;
3172  // cout << endl;
3173  // ^^^^^^^
3174 
3175  if(fDisplayOn == kTRUE) {
3176  TMarker *amrk2 = new TMarker(x2, y2, 21);
3177  amrk2->SetMarkerColor(5);
3178  display->cd(second + 1);
3179  amrk2->Draw("SAME");
3180  display->Update();
3181  display->Modified();
3182  }
3183  }
3184  }
3185  if(counter2 == nhitsonsensor[first + 1][0]) break;
3186  }
3187  if(counter1 == (nhitsonsensor[0][0] + nhitsonsensor[2][0] + nhitsonsensor[4][0])) break;
3188  }
3189 
3190  // cout << "MCPOINTS: " << endl;
3191 // for(int istat = 0; istat < fNPositions; istat++) {
3192 // cout << "istat " << istat << " " ;
3193 // for(int j = 0; j < mcpoints[istat].size(); j++) cout << " " << mcpoints[istat][j];
3194 // cout << endl;
3195 // }
3196 // cout << "ACCEPTED " << endl;
3197 // for(int istat = 0; istat < fNPositions; istat++) {
3198 // cout << "istat " << istat << " ";
3199 // for(int i = 0; i < accepted[istat].size(); i++) cout << " " << accepted[istat][i];
3200 // cout << endl;
3201 // }
3202 
3203 }
std::vector< int >::iterator fOrderingIterator
Int_t GetPosIndex(PndGemHit *hit)
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
TGeoVolume * sensor
static int counter2
Definition: createSTT.C:28
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
std::vector< int > fOrdering
Int_t GetStationNr() const
Definition: PndGemHit.h:81
Int_t GetDigiNr(Int_t iside) const
Definition: PndGemHit.h:77
PndSdsMCPoint * hit
Definition: anasim.C:70
static int counter1
Definition: createSTT.C:27
std::map< int, int > fCombiMap
void PndSttMvdGemTracking::Copy ( PndTrackCand completeCand,
PndTrack completeTrack,
PndTrackCand sttmvdCand,
PndTrack sttmvd 
)

Definition at line 759 of file PndSttMvdGemTracking.cxx.

References PndTrackCand::AddHit(), fCompleteTrackArray, fCompleteTrackCandArray, fVerbose, PndTrackCandHit::GetDetId(), PndTrack::GetFlag(), PndTrackCandHit::GetHitId(), PndTrackCand::getMcTrackId(), PndTrackCand::GetNHits(), PndTrack::GetParamFirst(), PndTrack::GetParamLast(), PndTrackCandHit::GetRho(), PndTrackCand::GetSortedHit(), PndTrackCand::GetSortedHits(), PndTrack::SetFlag(), PndTrackCand::setMcTrackId(), and PndTrack::SetRefIndex().

Referenced by Exec().

760  {
761 
762  // ------------------------------------------------- CHECK this has to be changed
763  // completeCand = sttmvdCand; // CHECK!
764  // CHECK ref index, to assign completeCand to sttmvdCand (itrk) URGENT
765 
766  // copy to the complete track (CHECK will be probably changed)
767  TClonesArray& clref = *fCompleteTrackCandArray;
768  Int_t size = clref.GetEntriesFast();
769  completeCand = new(clref[size]) PndTrackCand();
770 
771 
772 
773  completeCand->setMcTrackId(sttmvdCand->getMcTrackId());
774 
775  std::vector<PndTrackCandHit> sttmvdhits = sttmvdCand->GetSortedHits();
776  for(size_t ihit = 0; ihit < sttmvdCand->GetNHits(); ihit++) {
777  completeCand->AddHit(sttmvdCand->GetSortedHit(ihit).GetDetId(),
778  sttmvdCand->GetSortedHit(ihit).GetHitId(),
779  sttmvdCand->GetSortedHit(ihit).GetRho());
780 
781 
782  if(fVerbose > 0) cout << "STT + MVD iHit " << sttmvdCand->GetSortedHit(ihit).GetHitId() << " detId " << sttmvdCand->GetSortedHit(ihit).GetDetId() << endl;
783 
784  }
785 
786  TClonesArray& clref2 = *fCompleteTrackArray;
787  Int_t size2 = clref2.GetEntriesFast();
788  completeTrack = new(clref2[size2]) PndTrack(sttmvd->GetParamFirst(),
789  sttmvd->GetParamLast(),
790  *completeCand);
791  completeTrack->SetRefIndex("SttMvdGemTrackCand", size); // CHECK size or size2?
792  completeTrack->SetFlag(sttmvd->GetFlag());
793 }
int fVerbose
Definition: poormantracks.C:24
void SetRefIndex(TString branch, Int_t i)
Definition: PndTrack.h:41
int getMcTrackId() const
Definition: PndTrackCand.h:60
std::vector< PndTrackCandHit > GetSortedHits()
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
void setMcTrackId(int i)
Definition: PndTrackCand.h:72
Int_t GetFlag() const
Definition: PndTrack.h:33
TClonesArray * fCompleteTrackCandArray
TClonesArray * fCompleteTrackArray
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Double_t GetRho() const
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
void SetFlag(Int_t i)
Definition: PndTrack.h:38
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
Int_t GetHitId() const
Int_t GetDetId() const
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
Int_t PndSttMvdGemTracking::CountHitsInTrack ( Int_t  itrk)

Definition at line 1280 of file PndSttMvdGemTracking.cxx.

References counter, fVerbose, and trackvector.

Referenced by Exec().

1280  {
1281 
1282  int counter = 0;
1283  std::vector< std::pair<int, int> >::iterator iter;
1284  for(iter = trackvector.begin(); iter != trackvector.end(); iter++) {
1285  std::pair<int, int> thispair = *iter;
1286  if(thispair.first == itrk) counter++;
1287  }
1288  if(fVerbose > 0) std::cout << "track " << itrk << " has " << counter << " hits" << std::endl;
1289  return counter;
1290 }
int fVerbose
Definition: poormantracks.C:24
std::vector< std::pair< int, int > > trackvector
int counter
Definition: ZeeAnalysis.C:59
Int_t PndSttMvdGemTracking::CountTracks ( )

Definition at line 1260 of file PndSttMvdGemTracking.cxx.

References trackindexes.

Referenced by AddRemainingHits(), Exec(), OnlyOneHitToEachTrack(), and Retrack().

1260  {
1261  // probably this was wrong
1262  // int counter = 0;
1263 // std::vector< std::pair<int, int> >::iterator iter;
1264 // std::vector<int> counttracks;
1265 // std::vector<int>::iterator iter2;
1266 // for(iter = trackvector.begin(); iter != trackvector.end(); iter++) {
1267 // std::pair<int, int> thispair = *iter;
1268 // iter2 = find(counttracks.begin(), counttracks.end(), thispair.first);
1269 // Int_t where = iter2 - counttracks.begin();
1270 // if(where == counttracks.size()) {
1271 // counter++;
1272 // counttracks.push_back(thispair.first);
1273 // }
1274 // }
1275 // return counter;
1276 
1277  return trackindexes.size();
1278 }
std::vector< int > trackindexes
void PndSttMvdGemTracking::DeleteHitFromTrack ( Int_t  ihit,
Int_t  itrk 
)

Definition at line 1303 of file PndSttMvdGemTracking.cxx.

References fVerbose, notassignedhits, and trackvector.

Referenced by CheckCombinatorial(), ForbidMultiAssignedHits(), OnlyOneHitToEachTrack(), and Retrack().

1303  {
1304  std::pair<int, int> thispair(itrk, ihit);
1305  std::vector< std::pair<int, int> >::iterator iter;
1306  iter = std::find(trackvector.begin(), trackvector.end(), thispair);
1307  int where = iter - trackvector.begin();
1308  if(fVerbose != 0) std::cout << "where " << where << std::endl;
1309  trackvector.erase(iter);
1310  std::vector<int>::iterator iter2;
1311  iter2 = std::find(notassignedhits.begin(), notassignedhits.end(), ihit);
1312  int where2 = iter2 - notassignedhits.begin();
1313  if(where2 != (int)notassignedhits.size()) notassignedhits.push_back(ihit); // put it in not assigned list
1314  if(fVerbose > 0) cout << "DELETE HIT " << ihit << " FROM TRK " << itrk << endl;
1315 
1316 }
int fVerbose
Definition: poormantracks.C:24
std::vector< std::pair< int, int > > trackvector
std::vector< int > notassignedhits
void PndSttMvdGemTracking::EvaluatePerformances ( Int_t  nhits,
Int_t  ntracks 
)

Definition at line 796 of file PndSttMvdGemTracking.cxx.

References countbad, countdoubt, countgood, countnohitonplane, countplane, countprinotassigned, countprop, countreconstructablehit, countreconstructablepoint, countsecnotassigned, countsttmvd, countsttmvdusable, display, fabs(), fCompleteTrackCandArray, fDisplayOn, fGemBranchName, fGemHitArray, fGemPointArray, FillTrueDistances(), fMCTrackArray, fOrdering, fOrderingIterator, fTrackArray, fVerbose, PndTrackCandHit::GetDetId(), PndTrackCandHit::GetHitId(), PndTrackCand::getMcTrackId(), PndMCTrack::GetMotherID(), PndTrackCand::GetNHits(), GetPosIndex(), PndTrackCand::GetSortedHit(), GetTracksAssociatedToHit(), i, nhits, notassignedhits, pnt, refIndex, usabletracks, x0, y0, and z0.

Referenced by Exec().

796  {
797  // CHECK
798 
799  for(int itrk = 0; itrk < fCompleteTrackCandArray->GetEntriesFast(); itrk++) {
800  PndTrackCand *completeCand = (PndTrackCand* ) fCompleteTrackCandArray->At(itrk);
801  if(!completeCand) continue;
802  if(fVerbose > 0) cout << "complete cand " << itrk << " has nhits " << completeCand->GetNHits() << endl;
803 
804  // UpdateMCTrackId(completeCand);
805 
806  for (size_t ihit = 0; ihit < completeCand->GetNHits(); ihit++) {
807  PndTrackCandHit candhit = completeCand->GetSortedHit(ihit);
808  Int_t iHit = candhit.GetHitId();
809  Int_t detId = candhit.GetDetId();
810  // if(fVerbose != 0) // CHECK
811  if(fVerbose > 0) cout << "iHit " << iHit << " detId " << detId << "(" << FairRootManager::Instance()->GetBranchId(fGemBranchName) << ")" << endl;
812  if(detId != FairRootManager::Instance()->GetBranchId(fGemBranchName)) continue;
813 
814  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(iHit);
815  if(!gemhit) continue;
816 
817  // only for test.. will be deleted for real events ........CHECK
818  Int_t refIndex = gemhit->GetRefIndex();
819  if(refIndex == -1) {
820  countbad++; // background hit
821  if(fDisplayOn) {
822  TMarker *g2mrkb = new TMarker(gemhit->GetX(), gemhit->GetY(), 30);
823  int posindex = GetPosIndex(gemhit);
824  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
825  int index = fOrderingIterator - fOrdering.begin();
826  display->cd(index + 1);
827  g2mrkb->SetMarkerSize(2);
828  g2mrkb->SetMarkerColor(1);
829  g2mrkb->Draw("SAME");
830  display->Update();
831  display->Modified();
832  }
833  continue;
834  }
835 
836 
837  PndGemMCPoint *gempnt = (PndGemMCPoint*) fGemPointArray->At(refIndex);
838 
839  Int_t mcIndex = gempnt->GetTrackID();
840 
841  if(mcIndex == -1 || completeCand->getMcTrackId() == -1) countdoubt++;
842  else if(mcIndex == completeCand->getMcTrackId()) countgood++;
843  else if (mcIndex != completeCand->getMcTrackId()) countbad++;
844 
845  if(fDisplayOn) {
846  TMarker *g2mrk = new TMarker(gemhit->GetX(), gemhit->GetY(), 30);
847 
848  int posindex = GetPosIndex(gemhit);
849 
850  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
851  int index = fOrderingIterator - fOrdering.begin();
852  display->cd(index + 1);
853  g2mrk->SetMarkerSize(2);
854  if(mcIndex == completeCand->getMcTrackId()) g2mrk->SetMarkerColor(2); // red
855  else if (mcIndex != completeCand->getMcTrackId()) g2mrk->SetMarkerColor(1); // black
856  g2mrk->Draw("SAME");
857  display->Update();
858  display->Modified();
859  }
860 
861  }
862  }
863 
864  std::vector<int> recoTrack;
865  // the left hits have not been assigned
866  for(int itrk = 0; itrk < ntracks; itrk++) {
867  PndTrackCand *completeCand = (PndTrackCand* ) fCompleteTrackCandArray->At(itrk);
868  if(!completeCand) continue;
869  recoTrack.push_back(completeCand->getMcTrackId());
870  }
871 
872  for(int ihit = 0; ihit < nhits; ihit++) {
873  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(ihit);
874  if(!gemhit) continue;
875  Int_t refIndex = gemhit->GetRefIndex();
876  if(refIndex == -1) continue;
877  PndGemMCPoint *gempnt = (PndGemMCPoint*) fGemPointArray->At(refIndex);
878  if(!gempnt) continue;
879 
880  Int_t mcIndex = gempnt->GetTrackID();
881  if(mcIndex == -1) continue; // CHECK THIS
882  std::vector<int>::iterator recoTrackIter;
883  recoTrackIter = find(recoTrack.begin(), recoTrack.end(), mcIndex);
884  Int_t where= recoTrackIter - recoTrack.begin();
885 
886  // take only real hits: no combinatories CHECK
887  std::vector<int> truepoints;
888  if(fabs(gempnt->GetX() - gemhit->GetX()) < 1 && // CHECK if 1 is ok as distance mcpoint to hit
889  fabs(gempnt->GetY() - gemhit->GetY()) < 1) {
890  std::vector<int>::iterator iter;
891  iter = find(truepoints.begin(), truepoints.end(), refIndex);
892  Int_t where2 = iter - truepoints.begin();
893 
894 
895  // CHECK 4 PERFORMACE delete these .....
896  std::vector<int>::iterator iter3;
897  iter3 = find(usabletracks.begin(), usabletracks.end(), mcIndex);
898  Int_t where3 = iter3 - usabletracks.begin();
899  if(where != (int)recoTrack.size() && where2 == (int)truepoints.size() && where3 != (int)usabletracks.size()) {
900  // .............. and uncomment this
901  // if(where != recoTrack.size() && where2 == truepoints.size()) {
902 
904  truepoints.push_back(refIndex);
905 
906  if(fDisplayOn) {
907  int posindex = GetPosIndex(gemhit);
908  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
909  int ipos = fOrderingIterator - fOrdering.begin();
910 
911  TMarker *gpntmrk = new TMarker(gempnt->GetX(), gempnt->GetY(), 7);
912  gpntmrk->SetMarkerColor(5); //
913  display->cd(1 + ipos);
914  gpntmrk->Draw("SAME");
915  display->Update();
916  display->Modified();
917  }
918  }
919  }
920 
921  if(where != (int)recoTrack.size()) countreconstructablehit++; // if the hit is found
922 
923  if(GetTracksAssociatedToHit(ihit).size() != 0) continue;
924  PndMCTrack *mctrk = (PndMCTrack*) fMCTrackArray->At(mcIndex);
925  if(!mctrk) continue;
926  if(mctrk->GetMotherID() == -1) countprinotassigned++;
927  else countsecnotassigned++;
928  }
929 
930  for(size_t i = 0; i < notassignedhits.size(); i++) {
931  if(fVerbose != 0) cout << "hit " << notassignedhits[i] << "NOT ASSIGNED" << endl;
932  }
933 
934 
936 
937 
938  // draw ALL the MC points
939  for(int ipnt = 0; ipnt < fGemPointArray->GetEntriesFast(); ipnt++) {
941  if(!pnt) continue;
942  double x0 = pnt->GetX();
943  double y0 = pnt->GetY();
944  double z0 = pnt->GetZ();
945  int ipos = -1;
946  if(fDisplayOn) { // CHECK
947  if(z0 < 116.5) ipos = 0;
948  else if(z0 < 118.) ipos = 1;
949  else if(z0 < 152.4) ipos = 2;
950  else if(z0 < 154) ipos = 3;
951  else if(z0 < 188.4) ipos = 4;
952  else if(z0 < 190) ipos = 5;
953 
954  cout << "MC POINT " << x0 << " " << y0 << " " << z0 << " " << pnt->GetTrackID() << endl;
955 
956  TMarker *gpntmrk = new TMarker(x0, y0, 3);
957  gpntmrk->SetMarkerColor(7); //
958  display->cd(1 + ipos);
959  // cout << "DRAWING" << endl;
960  gpntmrk->Draw("SAME");
961  display->Update();
962  display->Modified();
963  }
964  }
965 
966  if(fVerbose > 0) {
967  cout << "SUMMARY OF EVENT" << endl;
968  cout << "mvd + stt tracks " << fTrackArray->GetEntriesFast() << endl;
969  cout << "gem hits " << fGemHitArray->GetEntriesFast() << endl;
970  cout << "complete track cands " << fCompleteTrackCandArray->GetEntriesFast() << endl;
971  }
972  cout << "COUNTERS " << countgood << " " << countbad << " " << countdoubt << " "
975  cout << "STTMVD " << countsttmvd << " of which usable " << countsttmvdusable << endl;
976 
977  cout << "planes: ";
978  for(int iplane = 0; iplane < 8; iplane++) cout << countplane[iplane] << " ";
979  cout << endl;
980  cout << "nohit on plane: ";
981  for(int iplane = 0; iplane < 8; iplane++) cout << countnohitonplane[iplane] << " ";
982  cout << endl;
983  cout << "prop 1 turn: ";
984  for(int iplane = 0; iplane < 8; iplane++) cout << countprop[0][iplane] << " ";
985  cout << endl;
986  cout << "prop 2 turn: ";
987  for(int iplane = 0; iplane < 8; iplane++) cout << countprop[1][iplane] << " ";
988  cout << endl;
989 
990 }
std::vector< int >::iterator fOrderingIterator
Double_t z0
Definition: checkhelixhit.C:62
Double_t x0
Definition: checkhelixhit.C:70
int fVerbose
Definition: poormantracks.C:24
Int_t GetPosIndex(PndGemHit *hit)
int getMcTrackId() const
Definition: PndTrackCand.h:60
Int_t i
Definition: run_full.C:25
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
TClonesArray * pnt
TClonesArray * fCompleteTrackCandArray
std::vector< int > fOrdering
Int_t refIndex
Definition: checkmomentum.C:30
Double_t y0
Definition: checkhelixhit.C:71
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
std::vector< int > GetTracksAssociatedToHit(Int_t ihit)
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
Int_t GetHitId() const
std::vector< int > notassignedhits
Int_t GetDetId() const
std::vector< int > usabletracks
void PndSttMvdGemTracking::Exec ( Option_t *  opt)
virtual

Virtual method Exec

Definition at line 487 of file PndSttMvdGemTracking.cxx.

References PndTrackCand::AddHit(), AddHitToTrack(), AddRemainingHits(), AssignHits(), Bool_t, CHECKCOMBI, CheckCombinatorial(), ConsiderCombinatorialEffect(), Copy(), CountHitsInTrack(), countnohitonplane, countplane, countprop, countsttmvd, countsttmvdusable, CountTracks(), distancemap, Double_t, EvaluatePerformances(), evt, fabs(), fCompleteTrackArray, fCompleteTrackCandArray, fDefaultPdgCode, fDisplayOn, fEvaluate, fGemBranchName, fGemHitArray, fGemParameters, fMCTrackArray, fMvdPixelBranchName, fMvdStripBranchName, fNPositions, ForbidMultiAssignedHits(), fPdgCode, fPdgFromMC, fProTracks, fSttBranchName, fTrackArray, fTurn, fVerbose, GetChargeCorrectedPdgFromMC(), GetClosestOnFirst(), GetHitsAssociatedToTrack(), PndTrackCand::getMcTrackId(), PndMCTrack::GetMotherID(), PndTrack::GetParamLast(), PndMCTrack::GetPdgCode(), PndGemHit::GetPosition(), PndTrack::GetRefIndex(), PndGemStation::GetSensor(), PndGemDigiPar::GetStation(), PndTrack::GetTrackCandPtr(), GetTrackIndex(), PndGemSensor::GetZ0(), i, nhits, OnlyOneHitToEachTrack(), OrderGemHits(), PndGemSensor::Print(), PropagateToGemPlane(), PropagateToGemPlaneAsHelix(), Reset(), Retrack(), sensor, PndTrack::SetFlag(), SetStartParameters(), PndTrack::SetTrackCand(), sttmvd, THROWAWAY, towhichplane, trackindexes, and usabletracks.

487  {
488  if(fVerbose > 0) {
489  cout << "==================== EVENT " << evt << endl;
490  cout << "detId " << FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName)
491  << " " << FairRootManager::Instance()->GetBranchId(fMvdStripBranchName)
492  << " " << FairRootManager::Instance()->GetBranchId(fGemBranchName)
493  << " " << FairRootManager::Instance()->GetBranchId(fSttBranchName) << endl;
494  }
495 
496 
497  evt++;
498  fTurn = 1;
499  fCompleteTrackCandArray->Delete();
500  fCompleteTrackArray->Delete();
501 
502  Int_t nhits = fGemHitArray->GetEntriesFast();
503  Int_t ntracks = fTrackArray->GetEntriesFast();
504 
505  if(fVerbose > 0) {
506  cout << "stt + mvd track array " << ntracks << endl;
507  cout << "gem hits " << nhits << endl;
508  }
509 
510  // if(nhits == 0) return; // if there is no GEM hit, there is no need to waste time! CHECK this has to be changed
511  if(ntracks == 0) return; // if there is no track, there is no need to waste time!
512  countsttmvd += ntracks;
513 
514  PndTrack *sttmvd = NULL;
515  PndTrackCand *sttmvdCand = NULL;
516  PndTrackCand *completeCand = NULL;
517  PndTrack *completeTrack = NULL;
518 
519  // RESET all the maps and counters for this event
520  Reset(nhits, ntracks);
521 
522  // Order GEM hits
523  OrderGemHits(nhits);
524  Int_t flag[ntracks];
525 
527 
528  // list of tracks which can/cannot be propagated
529  fProTracks.clear();
530 
531 
532  // loop on the tracks found in mvd + stt ************
533  for (Int_t itrk = 0; itrk < ntracks; itrk++)
534  {
535  fProTracks[itrk] = false;
536 
537  if(fVerbose > 0) cout << "----------- TRACK " << itrk << "------------" << endl;
538 
539  flag[itrk] = 0;
540  sttmvd = (PndTrack*) fTrackArray->At(itrk);
541  if (!sttmvd)
542  {
543  cout << "-W- PndSttMvdGemTracking::Exec: Empty PndTrack at " << itrk << endl;
544  continue;
545  }
546  sttmvdCand = sttmvd->GetTrackCandPtr();
547  if (!sttmvdCand)
548  {
549  cout << "-W- PndSttMvdGemTracking::Exec: Empty PndTrackCand at " << itrk << endl;
550  continue;
551  }
552 
553  if(fDisplayOn) {
554  char goOnChar;
555  cout << "press any key" << endl;
556  cin >> goOnChar;
557  cout << "GOING ON, nex track" << endl;
558  }
559 
560  // ------------------------------------------------- CHECK this has to be changed
561  Copy(completeCand, completeTrack, sttmvdCand, sttmvd);
562  trackindexes.push_back(itrk);
563  // -------------------------------------------------
564 
565  if(nhits == 0) { flag[itrk] = -1; continue; } // CHECK
566 
567  // cout << "copied completeCand from sttmvdCand @ " << itrk << " has hits " << completeCand->GetNHits() << endl;
568  //
569  FairTrackParP lastpar = sttmvd->GetParamLast();
570  if(lastpar.GetMomentum().Z() > 999990.) { flag[itrk] = -2; continue; } // CHECK
571  FairTrackParP *gempar = new FairTrackParP();
572 
573  if(fabs(lastpar.GetPosition().X()) > 42. || fabs(lastpar.GetPosition().Y()) > 42.) { flag[itrk] = -3; continue; } // CHECK 4 PERFORMANCE delete this!!!!
574  if(lastpar.GetMomentum().Mag() < 0.15) {
575  if(fVerbose > 0) {
576  cout << "TOO LOW MOMENTUM " << itrk << endl;
577  lastpar.GetPosition().Print();
578  lastpar.GetMomentum().Print();
579  }
580  flag[itrk] = 1;
581 
582  // continue; // CHECK 4 PERFORMANCE delete this!!!!
583  }
584 
585  // last z position
586  // get 1st sensor in 1st station
587  PndGemStation *station = fGemParameters->GetStation(0);
588  PndGemSensor *sensor = station->GetSensor(0);
589  if(lastpar.GetPosition().Z() > sensor->GetZ0()) {
590  if(fVerbose > 0) {
591  cout << "Z OUT OF BOUNDS" << endl;
592  lastpar.GetPosition().Print();
593  lastpar.GetMomentum().Print();
594  }
595  flag[itrk] = 2;
596 
597  // continue; // CHECK 4 PERFORMANCE delete this!!!!
598 
599  }
600 
601  if(fVerbose > 0) cout << "SetParameters for track " << itrk << endl;
602  FairTrackParP tmppar = SetStartParameters(sttmvd, sttmvdCand);
603 
604  if (fabs(tmppar.GetMomentum().Z()) < 1.e-5) {
605  if(fVerbose > 0) cout << " CANNOT PROPAGATE because z mom == 0" << endl;
606  flag[itrk] = -7;
607  continue;
608  }
609 
610  // =========== test of prop on 1st plane
611  FairTrackParP *gempartest = new FairTrackParP();
612  if(PropagateToGemPlaneAsHelix(sttmvd, gempartest, 0) == kFALSE) {
613  if(fVerbose > 0) cout << " CANNOT PROPAGATE " << endl;
614  delete gempartest;
615  flag[itrk] = -6;
616  continue;
617  }
618 
619  // ===========
620 
621  int charge = tmppar.GetQ();
622  // fPdgCode = -13 * charge;
623  if(fPdgFromMC) fPdgCode = GetChargeCorrectedPdgFromMC(itrk, charge);
624  else fPdgCode = fDefaultPdgCode * charge;
625 
626 // cout << "START MOM " << lastpar.GetMomentum().Mag() << " " << tmppar.GetMomentum().Mag() << endl;
627 
628  // **************
629  if(fEvaluate) {
630  int mcIndex = sttmvdCand->getMcTrackId();
631  if(mcIndex == -1) {
632  flag[itrk] = -4;
633  if(THROWAWAY) continue;
634  } // CHECK 4 PERFORMANCE delete this!!!!
635  else {
636  // cout << "FROM MC " << mcIndex << endl;
637  PndMCTrack *mctrk = (PndMCTrack*) fMCTrackArray->At(mcIndex);
638  if(mctrk) {
639  // cout << "PDG " << mctrk->GetPdgCode() << " MOTHERID " << mctrk->GetMotherID() << endl;
640  int mccharge = (int) TMath::Sign(1., TDatabasePDG::Instance()->GetParticle(mctrk->GetPdgCode())->Charge()/3.);
641  if(mccharge != charge && fVerbose > 0) cout << "WRONG CHARGE " << charge << " " << mccharge << " " << mctrk->GetPdgCode() << endl;
642 
643 
644  if(mctrk->GetMotherID() != -1) {
645  flag[itrk] = -5;
646  if(THROWAWAY) continue;
647  } // CHECK 4 PERFORMANCE delete this!!!!
648  }
649  }
650  // **************
652  usabletracks.push_back(mcIndex); // CHECK 4 PERFORMANCE delete this!!!!
653  }
654 
655  fProTracks[itrk] = true;
656 
657  Int_t closestonfirst = -1;
658  Double_t closestdistance = -1;
659  // loop over the GEM hits from the closest to 0, 0, 0 to the most external
660  // and
661  // propagate to track from mvd + stt to each of them
662  for(int ipos = 0; ipos < fNPositions; ipos++) {
663  // if no hit then continue
664  if(towhichplane[ipos] == false) { countnohitonplane[ipos]++; continue; }
665 
666  // ...else propagate here
667  Bool_t prop = PropagateToGemPlane(&tmppar, gempar, ipos);
668  if(prop == kFALSE) prop = PropagateToGemPlaneAsHelix(sttmvd, gempar, ipos);
669 
670  if (prop == kFALSE) {
671  countprop[fTurn-1][ipos]++;
672  continue; // CHECK (continue or break?)
673  }
674  countplane[ipos]++;
675 
676  // if the propagation was successful ...
677  // assign the hits to track itrk, extrapolated to gempar on ipos
678  std::vector<int> assignedhits = AssignHits(itrk, gempar, ipos);
679  for(size_t ihit = 0; ihit < assignedhits.size(); ihit++) AddHitToTrack(assignedhits[ihit], itrk);
680 
681  if(closestonfirst == -1) closestonfirst = GetClosestOnFirst(gempar, ipos, closestdistance);
682 
683  Double_t covMat[15]; gempar->GetCov(covMat); // CHECK
684  tmppar.SetTrackPar(gempar->GetX(), gempar->GetY(), gempar->GetZ(),
685  gempar->GetPx(), gempar->GetPy(), gempar->GetPz(), gempar->GetQ(),
686  covMat,
687  gempar->GetOrigin(), gempar->GetIVer(), gempar->GetJVer(), gempar->GetKVer());
688  }
689 
690  // use closest on first CHECK with more tracks
691  if(GetHitsAssociatedToTrack(itrk).size() == 0 && closestonfirst != -1) {
692  // cout << "CLOSEST HIT " << closestonfirst << " " << closestdistance << endl;
693  distancemap[itrk][closestonfirst] = closestdistance;
694  AddHitToTrack(closestonfirst, itrk);
695  }
696 
697 
698  }
699 
700  if(nhits == 0) return;
701 
702  if(fVerbose > 0 && CountTracks() != fCompleteTrackCandArray->GetEntriesFast()) cout << "ERROR!!! " << CountTracks() << " " << fCompleteTrackCandArray->GetEntriesFast() << endl;
703 
704  // CLEANUP =================================
705  if(CHECKCOMBI == 1) CheckCombinatorial(nhits, ntracks);
706  ForbidMultiAssignedHits(nhits, ntracks);
707  OnlyOneHitToEachTrack(nhits, ntracks);
708  AddRemainingHits(ntracks);
709  fTurn = 2;
710  Retrack();
711 
712 
713  // CHECK to test
714  if(fVerbose > 0) {
715  cout << "N OF TRACKS " << CountTracks() << endl;
716  for(int i = 0; i < CountTracks(); i++) {
717  int itrk = GetTrackIndex(i);
718  Int_t nofhits = CountHitsInTrack(itrk);
719  cout << "track " << itrk << " has " << nofhits << endl;
720  std::vector<int> hitsoftrack;
721  hitsoftrack = GetHitsAssociatedToTrack(itrk);
722  for(int ihit = 0; ihit < nofhits; ihit++) cout << " " << hitsoftrack[ihit];
723  cout << endl;
724  }
725  }
726  // ----------------------------------
727 
728  // ACTUALLY ASSIGN HITS TO TRACK
729  PndGemHit *gemhit = NULL;
730  for(Int_t i = 0; i < CountTracks(); i++) {
731  int itrk = GetTrackIndex(i);
732  if(fVerbose > 0) cout << i << " ACTUALLY ASSIGNING TRACK " << itrk << endl;
733  completeCand = (PndTrackCand* ) fCompleteTrackCandArray->At(itrk);
734  std::vector<int> thistrackhits = GetHitsAssociatedToTrack(itrk);
735 
736  for(size_t j = 0; j < thistrackhits.size(); j++) {
737  int ihit = thistrackhits[j];
738  gemhit = (PndGemHit*) fGemHitArray->At(ihit);
739  if(!gemhit) continue;
740  if(fVerbose > 0) cout << "ACTUALLY add hit " << ihit << " to track " << itrk << endl;
741  completeCand->AddHit(FairRootManager::Instance()->GetBranchId(fGemBranchName), ihit, gemhit->GetPosition().Mag()); // CHECK rho and kGemHit
742 
743  }
744 
745  completeTrack = (PndTrack*) fCompleteTrackArray->At(itrk);
746  if(fVerbose > 0 && completeTrack->GetRefIndex() != itrk) cout << "************** ERROR ****************" << endl;
747  completeTrack->SetTrackCand(*completeCand);
748  completeTrack->SetFlag(flag[itrk]);
749  // cout << "track " << itrk << " has flag " << flag[itrk] << endl;
750  }
751 
752 
753  // performance
754  if(fEvaluate) EvaluatePerformances(nhits, ntracks);
755 
756 }
int fVerbose
Definition: poormantracks.C:24
PndGemDigiPar * fGemParameters
#define THROWAWAY
void OrderGemHits(Int_t nhits)
int getMcTrackId() const
Definition: PndTrackCand.h:60
void Copy(PndTrackCand *completeCand, PndTrack *completeTrack, PndTrackCand *sttmvdCand, PndTrack *sttmvd)
Int_t i
Definition: run_full.C:25
void CheckCombinatorial(Int_t nhits, Int_t ntracks)
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
std::vector< int > trackindexes
TGeoVolume * sensor
Int_t CountHitsInTrack(Int_t itrk)
TClonesArray * fCompleteTrackCandArray
Int_t GetChargeCorrectedPdgFromMC(int trackid, int charge)
TClonesArray * fCompleteTrackArray
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
Int_t GetClosestOnFirst(FairTrackParP *gempar, Int_t ipos, Double_t &closestdistance)
std::vector< int > GetHitsAssociatedToTrack(Int_t itrk)
void EvaluatePerformances(Int_t nhits, Int_t ntracks)
void AddRemainingHits(Int_t ntracks)
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Bool_t PropagateToGemPlane(FairTrackParP *tmppar, FairTrackParP *gempar, Int_t ipos)
Double_t
TMatrixT< double > distancemap
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
std::map< int, bool > fProTracks
#define CHECKCOMBI
void ForbidMultiAssignedHits(Int_t nhits, Int_t ntracks)
void AddHitToTrack(Int_t ihit, Int_t itrk)
Bool_t PropagateToGemPlaneAsHelix(PndTrack *sttmvd, FairTrackParP *gempar, Int_t ipos)
TVector3 GetPosition() const
Definition: PndGemHit.h:71
void SetFlag(Int_t i)
Definition: PndTrack.h:38
FairTrackParP SetStartParameters(PndTrack *sttmvd, PndTrackCand *sttmvdCand)
PndSttMvdTracking * sttmvd
Definition: runrecoMix.C:114
void OnlyOneHitToEachTrack(Int_t nhits, Int_t ntracks)
Double_t GetZ0() const
Definition: PndGemSensor.h:100
std::vector< int > AssignHits(Int_t itrk, FairTrackParP *gempar, Int_t ipos)
Int_t GetRefIndex() const
Definition: PndTrack.h:36
PndTrackCand * GetTrackCandPtr()
Definition: PndTrack.h:48
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
void SetTrackCand(const PndTrackCand &cand)
Definition: PndTrack.h:43
void ConsiderCombinatorialEffect(Int_t nhits)
std::map< int, bool > towhichplane
PndGemStation * GetStation(Int_t iStation)
void Reset(Int_t nhits, Int_t ntracks)
std::vector< int > usabletracks
void PndSttMvdGemTracking::FillTrueDistances ( )

Definition at line 2133 of file PndSttMvdGemTracking.cxx.

References Bool_t, fabs(), fDefaultPdgCode, fGemPointArray, fMCTrackArray, fMvdPixelBranchName, fMvdPixelHitArray, fMvdPointArray, fMvdStripBranchName, fMvdStripHitArray, fNPositions, fPdgCode, fPdgFromMC, fProTracks, fSttBranchName, fSttHitArray, fSttPointArray, fTrackArray, GetChargeCorrectedPdgFromMC(), PndTrackCandHit::GetDetId(), PndTrackCandHit::GetHitId(), PndTrackCand::getMcTrackId(), PndMCTrack::GetMomentum(), PndTrackCand::GetNHits(), PndTrackCand::GetSortedHit(), PndMCTrack::GetStartVertex(), PndTrack::GetTrackCandPtr(), hmcdist, hmcx, hmcy, PropagateToGemPlane(), PropagateToGemPlaneAsHelix(), and SetStartParameters().

Referenced by EvaluatePerformances().

2133  {
2134 
2135  // loop on tracks
2136  // cout << "FOUND TRACKS " << fTrackArray->GetEntriesFast() << endl;
2137  for(Int_t itrk = 0; itrk < fTrackArray->GetEntriesFast(); itrk++) {
2138  PndTrack* sttmvdTrack = (PndTrack*) fTrackArray->At(itrk);
2139  if(!sttmvdTrack) continue;
2140  PndTrackCand *sttmvdCand = sttmvdTrack->GetTrackCandPtr();
2141  if(!sttmvdCand) continue;
2142  if(fProTracks[itrk] == kFALSE) continue;
2143  Int_t mcIndex = sttmvdCand->getMcTrackId();
2144  // cout << "TRACK " << itrk << " HAS MC " << mcIndex << endl;
2145  if(mcIndex == -1) continue;
2146  PndMCTrack *mctrk = (PndMCTrack*) fMCTrackArray->At(mcIndex);
2147  if(mctrk) {
2148  //int pdgcode = mctrk->GetPdgCode(); //[R.K. 01/2017] unused variable?
2149  //int motherid = mctrk->GetMotherID(); //[R.K. 01/2017] unused variable?
2150  TVector3 vertex = mctrk->GetStartVertex();
2151  TVector3 vtxmomentum = mctrk->GetMomentum();
2152  // cout << "PDG " << pdgcode << " MOTHID " << motherid << " VTX POS " << vertex.Mag() << " MOM " << vtxmomentum.Mag() << endl;
2153 
2154  // start from MC last point on STT
2155  Int_t nsttmvdhits = sttmvdCand->GetNHits();
2156  PndTrackCandHit candhit = sttmvdCand->GetSortedHit(nsttmvdhits - 1);
2157  Int_t iHit = candhit.GetHitId();
2158  Int_t detId = candhit.GetDetId();
2159  FairHit *fhit = NULL;
2160  FairMCPoint *fpnt = NULL;
2161  // cout << "HIT " << iHit << " FOR " << detId << endl;
2162  if(detId == FairRootManager::Instance()->GetBranchId(fMvdStripBranchName)) { // CHECK
2163  // if(detId == 26) {
2164  fhit = (FairHit*) fMvdStripHitArray->At(iHit);
2165  if(fhit->GetRefIndex() != -1 && fhit) fpnt = (FairMCPoint*) fMvdPointArray->At(fhit->GetRefIndex());
2166  }
2167  else if(detId == FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName)) { // CHECK
2168  // else if(detId == 24) {
2169  fhit = (FairHit*) fMvdPixelHitArray->At(iHit);
2170  if(fhit->GetRefIndex() != -1 && fhit) fpnt = (FairMCPoint*) fMvdPointArray->At(fhit->GetRefIndex());
2171  }
2172  else if(detId == FairRootManager::Instance()->GetBranchId(fSttBranchName)) {
2173  fhit = (FairHit*) fSttHitArray->At(iHit);
2174  if(fhit->GetRefIndex() != -1 && fhit) fpnt = (FairMCPoint*) fSttPointArray->At(fhit->GetRefIndex());
2175  }
2176  if(fpnt) {
2177  TVector3 position(fpnt->GetX(), fpnt->GetY(), fpnt->GetZ());
2178  TVector3 momentum(fpnt->GetPx(), fpnt->GetPy(), fpnt->GetPz());
2179 
2180 // cout << "MC/PR POS" << endl;
2181 // position.Print();
2182 // sttmvdTrack->GetParamLast().GetPosition().Print();
2183 // cout << "MC/PR MOM" << endl;
2184 // momentum.Print();
2185 // sttmvdTrack->GetParamLast().GetMomentum().Print();
2186  }
2187  }
2188  FairTrackParP *gempar = new FairTrackParP();
2189  FairTrackParP tmppar = SetStartParameters(sttmvdTrack, sttmvdCand);
2190  int charge = tmppar.GetQ();
2191  // fPdgCode = -13 * charge;
2192  if(fPdgFromMC) fPdgCode = GetChargeCorrectedPdgFromMC(itrk, charge);
2193  else fPdgCode = fDefaultPdgCode * charge;
2194 
2195  // extrapolate each track on each plane
2196  for(int ipos = 0; ipos < fNPositions; ipos++) {
2197 
2198  Bool_t prop = PropagateToGemPlane(&tmppar, gempar, ipos);
2199  if(prop == kFALSE) prop = PropagateToGemPlaneAsHelix(sttmvdTrack, gempar, ipos);
2200  if (prop == kFALSE) continue; // CHECK (continue or break?)
2201  TVector3 extrappos = gempar->GetPosition();
2202 
2203  // loop on MC points
2204  for(int ipnt = 0; ipnt < fGemPointArray->GetEntriesFast(); ipnt++) {
2205  PndGemMCPoint *gempnt = (PndGemMCPoint*) fGemPointArray->At(ipnt);
2206  if(!gempnt) continue;
2207  Int_t pntMcIndex = gempnt->GetTrackID();
2208  // if the point belongs to the mctrack
2209  if(pntMcIndex != mcIndex) continue;
2210  TVector3 gempos(gempnt->GetX(), gempnt->GetY(), gempnt->GetZ());
2211  // if they are on the same plane
2212  if(fabs(gempos.Z() - extrappos.Z()) < 0.1) {
2213  double distance = (gempos - extrappos).Perp();
2214  hmcdist[ipos]->Fill(distance);
2215  hmcx[ipos]->Fill(gempos.X() - extrappos.X());
2216  hmcy[ipos]->Fill(gempos.Y() - extrappos.Y());
2217  // cout << "POINT " << ipnt << " TO TRACK " << itrk << "(MC " << mcIndex << ") distance " << distance << endl;
2218  }
2219  }
2220  }
2221  }
2222 }
int getMcTrackId() const
Definition: PndTrackCand.h:60
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
Int_t GetChargeCorrectedPdgFromMC(int trackid, int charge)
TClonesArray * fMvdPixelHitArray
Bool_t PropagateToGemPlane(FairTrackParP *tmppar, FairTrackParP *gempar, Int_t ipos)
TClonesArray * fMvdStripHitArray
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
std::map< int, bool > fProTracks
Bool_t PropagateToGemPlaneAsHelix(PndTrack *sttmvd, FairTrackParP *gempar, Int_t ipos)
FairTrackParP SetStartParameters(PndTrack *sttmvd, PndTrackCand *sttmvdCand)
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndTrackCand * GetTrackCandPtr()
Definition: PndTrack.h:48
Int_t GetHitId() const
Int_t GetDetId() const
Bool_t PndSttMvdGemTracking::Fit ( TMatrixT< double >  points,
Double_t outxc,
Double_t outyc,
Double_t outradius 
)

Definition at line 2562 of file PndSttMvdGemTracking.cxx.

References a, alpha, CAMath::ATan2(), b, c, CAMath::Cos(), Double_t, fabs(), fSttBranchName, nhits, R, s, CAMath::Sin(), sqrt(), and v.

Referenced by Prefit().

2563 {
2564 
2565  Int_t nhits = points.GetRowUpb() + 1;
2566  int lasthitid = -1, firsthitid = -1;
2567  for(int ihit = 0; ihit < nhits; ihit++) {
2568  if(points[ihit][0] != -1 && points[ihit][10] != -1 && points[ihit][10] != 1) {
2569  if(firsthitid == -1) firsthitid = ihit;
2570  lasthitid = ihit;
2571  }
2572  }
2573 
2574  if(firsthitid == -1 || lasthitid == -1) return kFALSE; // CHECK
2575 
2576  Double_t trasl[2] = {points[firsthitid][2], points[firsthitid][3]};
2577 
2578  // cout << "first/last " << firsthitid << " " << lasthitid << endl;
2579  if(firsthitid >= lasthitid) return false;
2580  Double_t alpha = TMath::ATan2(points[lasthitid][3] - points[firsthitid][3],
2581  points[lasthitid][2] - points[firsthitid][2]);
2582  Double_t Suu, Su, Sv, Suv, S1, Suuu, Suuv, Suuuu;
2583 
2584  Su = 0.;
2585  Sv = 0.;
2586  Suu = 0.;
2587  Suv = 0.;
2588  Suuu = 0.;
2589  S1 = 0.;
2590  Suuv = 0.;
2591  Suuuu = 0.;
2592  Double_t s = 0.001; // CHECK
2593  // ..............................................
2594 
2595  TVector3 fitpoint;
2596  for(int ihit = 0; ihit < nhits; ihit++)
2597  {
2598  Int_t hitId = (Int_t) points[ihit][0];
2599  Int_t detId = (Int_t) points[ihit][1];
2600  if(hitId == -1) continue;
2601  Int_t fitflag = (Int_t) points[ihit][10];
2602  if(fitflag == 1 || fitflag == -1) continue;
2603  if(detId == FairRootManager::Instance()->GetBranchId(fSttBranchName) && points[ihit][8] < 0.1) continue;
2604  fitpoint.SetXYZ(points[ihit][2], points[ihit][3], points[ihit][4]);
2605  Double_t sigx = points[ihit][5];
2606  Double_t sigy = points[ihit][6];
2607 // cout << "fitpoint" << endl;
2608 // fitpoint.Print();
2609  // to the fit ================================
2610  Double_t xtrasl, ytrasl;
2611  // traslation
2612  xtrasl = fitpoint.X() - trasl[0];
2613  ytrasl = fitpoint.Y() - trasl[1];
2614 
2615  Double_t xrot, yrot;
2616  // rotation
2617  xrot = TMath::Cos(alpha)*xtrasl + TMath::Sin(alpha)*ytrasl;
2618  yrot = -TMath::Sin(alpha)*xtrasl + TMath::Cos(alpha)*ytrasl;
2619 
2620  // re-traslation
2621  xtrasl = xrot + s;
2622  ytrasl = yrot;
2623 
2624  // change coordinate
2625  Double_t u, v, sigv2; //, sigu2; //[R.K.02/2017] Unused variable?
2626  u = xtrasl / (xtrasl*xtrasl + ytrasl*ytrasl);
2627  v = ytrasl / (xtrasl*xtrasl + ytrasl*ytrasl);
2628 
2629  Double_t dvdx = (-2 * xtrasl * ytrasl)/pow((xtrasl*xtrasl + ytrasl*ytrasl),2);
2630  Double_t dvdy = (xtrasl*xtrasl - ytrasl*ytrasl) / pow((xtrasl*xtrasl + ytrasl*ytrasl),2);
2631  //Double_t dudx = (ytrasl*ytrasl - xtrasl*xtrasl) / pow((xtrasl*xtrasl + ytrasl*ytrasl),2); //[R.K.02/2017] Unused variable?
2632  //Double_t dudy = (-2 * xtrasl * ytrasl)/pow((xtrasl*xtrasl + ytrasl*ytrasl),2); //[R.K.02/2017] Unused variable?
2633 
2634  //sigu2 = dudx * dudx * sigx * sigx + dudy * dudy * sigy * sigy + 2 * dudx * dudy * sigx * sigy; //[R.K.02/2017] Unused variable?
2635  sigv2 = dvdx * dvdx * sigx * sigx + dvdy * dvdy * sigy * sigy + 2 * dvdx * dvdy * sigx * sigy;
2636 
2637  if(sigv2 == 0) sigv2 = 1e-5; // CHECK MVD covariance
2638 
2639  Su = Su + (u/sigv2);
2640  Sv = Sv + (v/sigv2);
2641 
2642  Suv = Suv + ((u*v)/sigv2);
2643  Suu = Suu + ((u*u)/sigv2);
2644 
2645  Suuu = Suuu + ((u*u*u)/sigv2);
2646  Suuv = Suuv + ((u*u*v)/sigv2);
2647 
2648  Suuuu = Suuuu + ((u*u*u*u)/sigv2);
2649 
2650  S1 = S1 + 1/sigv2;
2651  }
2652 
2653 
2654  TMatrixT<double> matrix(3,3);
2655  matrix[0][0] = S1;
2656  matrix[0][1] = Su;
2657  matrix[0][2] = Suu;
2658 
2659  matrix[1][0] = Su;
2660  matrix[1][1] = Suu;
2661  matrix[1][2] = Suuu;
2662 
2663  matrix[2][0] = Suu;
2664  matrix[2][1] = Suuu;
2665  matrix[2][2] = Suuuu;
2666 
2667  Double_t determ;
2668 
2669  determ = matrix.Determinant();
2670 
2671  if (determ != 0) {
2672  matrix.Invert();
2673  }
2674  else {
2675  // cout << "DET 0" << endl; // CHECK what to do
2676  return false;
2677  }
2678 
2679  TMatrixT<double> column(3,1);
2680  column[0][0] = Sv;
2681  column[1][0] = Suv;
2682  column[2][0] = Suuv;
2683 
2684  TMatrixT<double> column2(3,1);
2685  column2.Mult(matrix, column);
2686 
2687  Double_t a, b, c;
2688  a = column2[0][0];
2689  b = column2[1][0];
2690  c = column2[2][0];
2691 
2692  if(fabs(a)<0.000001) {
2693  // cout << "A < 1e-**" << endl;
2694  return kFALSE;
2695  }
2696 
2697  // center and radius
2698  Double_t xcrot, ycrot, xc, yc, epsilon, R;
2699  ycrot = 1/(2*a);
2700  xcrot = -b/(2*a);
2701  epsilon = -c*pow((1+(b*b)), -3/2);
2702  R = epsilon + sqrt((xcrot*xcrot)+(ycrot*ycrot));
2703 
2704  // re-rotation and re-traslation of xc and yc
2705  // translation
2706  xcrot = xcrot - s;
2707  // rotation
2708  xc = TMath::Cos(alpha)*xcrot - TMath::Sin(alpha)*ycrot;
2709  yc = TMath::Sin(alpha)*xcrot + TMath::Cos(alpha)*ycrot;
2710  // traslation
2711  xc = xc + trasl[0];
2712  yc = yc + trasl[1];
2713  //Double_t phi = TMath::ATan2(yc, xc); //[R.K.02/2017] Unused variable?
2714  //Double_t d; //[R.K.02/2017] Unused variable?
2715  //d = ((xc + yc) - R*(TMath::Cos(phi) + TMath::Sin(phi)))/(TMath::Cos(phi) + TMath::Sin(phi)); //[R.K.02/2017] Unused variable?
2716 
2717  // cout << "REFITTED FIT: " << xc << " " << yc << endl;
2718  // cout << "RAGGIO: " << R << endl;
2719 
2720  outxc = xc;
2721  outyc = yc;
2722  outradius = R;
2723  return true;
2724 }
TTree * b
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
TLorentzVector s
Definition: Pnd2DStar.C:50
static T Sin(const T &x)
Definition: PndCAMath.h:42
static T Cos(const T &x)
Definition: PndCAMath.h:43
__m128 v
Definition: P4_F32vec4.h:4
Int_t a
Definition: anaLmdDigi.C:126
Double_t
static T ATan2(const T &y, const T &x)
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
double alpha
Definition: f_Init.h:9
Double_t R
Definition: checkhelixhit.C:61
void PndSttMvdGemTracking::ForbidMultiAssignedHits ( Int_t  nhits,
Int_t  ntracks 
)

Definition at line 1122 of file PndSttMvdGemTracking.cxx.

References DeleteHitFromTrack(), distancemap, Double_t, fVerbose, GetTracksAssociatedToHit(), and nhits.

Referenced by Exec().

1123 {
1124  if(fVerbose > 0) cout << "FORBID MULTI ASSIGNED HITS : DELETING HITS" << endl;
1125  // ... use distancemap[itrk][ihit]
1126  for(int ihit = 0; ihit < nhits; ihit++) {
1127  if(fVerbose > 0) cout << "hit " << ihit << " associated to " << GetTracksAssociatedToHit(ihit).size() << " tracks" << endl;
1128  // if it is assigned to no track (0) or one track (1) we are done
1129  if(GetTracksAssociatedToHit(ihit).size() <= 1) continue;
1130 
1131  // else pick the right track
1132  Double_t tmpdistance = 1000; // CHECK set this to maxdistance
1133  Int_t trkcounter = 0;
1134  Int_t tmptrack = -1;
1135 
1136  int tracksassociated = GetTracksAssociatedToHit(ihit).size();
1137 
1138  // loop over the tracks
1139  for(int jtrk = 0; jtrk < tracksassociated; jtrk++) {
1140  Int_t itrk = GetTracksAssociatedToHit(ihit)[trkcounter];
1141  // cout << "track " << itrk << " distant " << distancemap[itrk][ihit] << " is " << trkcounter << endl;
1142  trkcounter++;
1143  if(distancemap[itrk][ihit] < tmpdistance) { // ... substitute this trk to tmp trk
1144  if(tmptrack != -1) {
1145  DeleteHitFromTrack(ihit, tmptrack);
1146  if(fVerbose > 0) cout << "hit " << ihit << " deleted from track " << tmptrack << endl;
1147  trkcounter--;
1148  }
1149  tmpdistance = distancemap[itrk][ihit];
1150  tmptrack = itrk;
1151  }
1152  else { // ... remove this hit definitively
1153  if(fVerbose > 0) cout << "hit " << ihit << " deleted from track " << itrk << endl;
1154  DeleteHitFromTrack(ihit, itrk);
1155  trkcounter--;
1156  }
1157  }
1158  }
1159 }
int fVerbose
Definition: poormantracks.C:24
Double_t
TMatrixT< double > distancemap
void DeleteHitFromTrack(Int_t ihit, Int_t itrk)
std::vector< int > GetTracksAssociatedToHit(Int_t ihit)
Int_t PndSttMvdGemTracking::GetChargeCorrectedPdgFromMC ( int  trackid,
int  charge 
)

Definition at line 3505 of file PndSttMvdGemTracking.cxx.

References GetPdgFromMC().

Referenced by Exec(), and FillTrueDistances().

3505  {
3506  int pdg = GetPdgFromMC(trackid);
3507  // is the reco track in accordance with the mc one?
3508  double mccharge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/3.;
3509  if((charge * mccharge) > 0.) return pdg;
3510  else return -pdg;
3511 }
Int_t GetPdgFromMC(int trackid)
Int_t PndSttMvdGemTracking::GetClosestOnFirst ( FairTrackParP *  gempar,
Int_t  ipos,
Double_t closestdistance 
)

Definition at line 2224 of file PndSttMvdGemTracking.cxx.

References CHECKCOMBI, Double_t, fCombiMap, fGemHitArray, PndGemHit::GetPosition(), hitcounter, and hitmap.

Referenced by Exec().

2224  {
2225  int hitonsensor = (int) hitcounter(ipos, 0);
2226  Int_t tmphit = -1;
2227  Double_t tmpdistance = 10000;
2228  for(int ihit = 0; ihit < hitonsensor; ihit++) {
2229  int hitindex = (int) hitmap(ipos, ihit);
2230  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(hitindex);
2231  if(!gemhit) continue;
2232  if(CHECKCOMBI == 2 && fCombiMap[hitindex] != 0) continue;
2233  TVector3 gemhitpos = gemhit->GetPosition();
2234  TVector3 extrapos = gempar->GetPosition();
2235 
2236  Double_t distance = (gemhitpos - extrapos).Perp();
2237  if(distance < tmpdistance) {
2238  tmpdistance = distance;
2239  tmphit = hitindex;
2240  }
2241  }
2242  closestdistance = tmpdistance;
2243 
2244  if(closestdistance < 0) return -1;
2245  return tmphit;
2246 }
TMatrixT< float > hitmap
Double_t
TMatrixT< float > hitcounter
#define CHECKCOMBI
TVector3 GetPosition() const
Definition: PndGemHit.h:71
std::map< int, int > fCombiMap
Int_t PndSttMvdGemTracking::GetHitIndex ( Int_t  i)

Definition at line 1243 of file PndSttMvdGemTracking.cxx.

References counter, and trackvector.

1243  { // CHECK THIS HAS TO BE TESTED!!
1244  int counter = -1;
1245  int tmphit = -1;
1246  std::vector< std::pair<int, int> >::iterator iter;
1247  for(iter = trackvector.begin(); iter != trackvector.end(); iter++) {
1248  std::pair<int, int> thispair = *iter;
1249  if(thispair.second != tmphit) {
1250  tmphit = thispair.second;
1251  counter++;
1252  // cout << i << " " << counter << " " << tmphit << " " << trackvector[i].second << endl;
1253  if(counter == i) return tmphit;
1254  }
1255  }
1256  cout << "PndSttMvdGemTracking::GetTrackIndex " << i << " Out Of Bounds" << endl;
1257  return -1;
1258 }
Int_t i
Definition: run_full.C:25
std::vector< std::pair< int, int > > trackvector
int counter
Definition: ZeeAnalysis.C:59
std::vector< int > PndSttMvdGemTracking::GetHitsAssociatedToTrack ( Int_t  itrk)

Definition at line 1334 of file PndSttMvdGemTracking.cxx.

References fVerbose, and trackvector.

Referenced by AddRemainingHits(), CheckCombinatorial(), Exec(), GetHitsAssociatedToTrackOnPlane(), OnlyOneHitToEachTrack(), and Retrack().

1334  {
1335  std::vector<int> thistrackhits;
1336  std::vector< std::pair<int, int> >::iterator iter;
1337  for(iter = trackvector.begin(); iter != trackvector.end(); iter++) {
1338  std::pair<int, int> thispair = *iter;
1339  if(thispair.first == itrk) thistrackhits.push_back(thispair.second);
1340  }
1341  if(fVerbose != 0) {
1342  std::cout << "track " << itrk << " has " << thistrackhits.size() << " hits: " ;
1343  for(size_t j = 0; j < thistrackhits.size(); j++) std::cout << thistrackhits[j] << " ";
1344  std::cout << std::endl;
1345  }
1346  return thistrackhits;
1347 }
int fVerbose
Definition: poormantracks.C:24
std::vector< std::pair< int, int > > trackvector
std::vector< int > PndSttMvdGemTracking::GetHitsAssociatedToTrackOnPlane ( Int_t  itrk,
Int_t  ipos 
)

Definition at line 1349 of file PndSttMvdGemTracking.cxx.

References Bool_t, GetHitsAssociatedToTrack(), hitcounter, and hitmap.

Referenced by Retrack().

1349  {
1350  // cout << "GetHitsAssociatedToTrackOnPlane " << itrk << " " << ipos << endl;
1351  int hitonsensor = (int) hitcounter(ipos, 0);
1352  std::vector<int> thistrackhitsonplane = GetHitsAssociatedToTrack(itrk);
1353  // cout << "total associated hits " << thistrackhitsonplane.size() << endl;
1354  std::vector<int>::iterator iter;
1355  for(iter = thistrackhitsonplane.begin(); iter != thistrackhitsonplane.end(); iter++) {
1356  Int_t ihit = *iter;
1357  // cout << "loop on " << hitonsensor << " hits" << endl;
1358  Bool_t hitfound = kFALSE;
1359  for(int jhit = 0; jhit < hitonsensor; jhit++) {
1360  Int_t hitindex = (int) hitmap(ipos, jhit);
1361  // cout << "ihit: ipos, jhit, hitmap " << ihit << " " << ipos << " " << jhit << " " << hitindex << endl;
1362  if(hitindex == ihit) hitfound = kTRUE;
1363  }
1364  if(hitfound == kFALSE) {
1365  // cout << "erase " << ihit << endl;
1366  thistrackhitsonplane.erase(iter);
1367  iter--;
1368  }
1369  }
1370  return thistrackhitsonplane;
1371 }
TMatrixT< float > hitmap
std::vector< int > GetHitsAssociatedToTrack(Int_t itrk)
TMatrixT< float > hitcounter
Bool_t PndSttMvdGemTracking::GetInitialParams ( PndTrack sttmvd,
Double_t xc,
Double_t yc,
Double_t radius,
Double_t fitm,
Double_t fitp 
)

Definition at line 2805 of file PndSttMvdGemTracking.cxx.

References CAMath::ATan2(), CalculatePhi(), CompareToPreviousPhi(), CAMath::Cos(), d, Double_t, fabs(), fVerbose, PndTrack::GetParamFirst(), PndTrack::GetParamLast(), p1, p2, phi, Pi, CAMath::Sin(), CAMath::Sqrt(), v, x0, and y0.

Referenced by Prefit(), and PropagateToGemPlaneAsHelix().

2806 {
2807  if(fVerbose > 0) cout << "GET INITIAL PARAMS" << endl;
2808  FairTrackParP recopar = sttmvd->GetParamFirst();
2809  TVector3 recomom = recopar.GetMomentum();
2810  TVector3 recopos = recopar.GetPosition();
2811  Int_t charge = recopar.GetQ();
2812 
2813  radius = recomom.Perp()/0.006;
2814  Double_t beta;
2815 
2816  if(fabs(recomom.X()) > 1e-10) {
2817  // track from tangent ---------------------
2818  //double reco_m1 = recomom.Y() / recomom.X(); //[R.K. 01/2017] unused variable?
2819  //double reco_q1 = recopos.Y() - recopos.X() * reco_m1; //[R.K. 01/2017] unused variable?
2820  //double reco_m2 = -1./reco_m1; //[R.K. 01/2017] unused variable?
2821  //double reco_q2 = recopos.Y() - recopos.X() * reco_m2; //[R.K. 01/2017] unused variable?
2822  beta = TMath::ATan2(recomom.X(), recomom.Y());
2823  }
2824  else beta = TMath::Sign(1., recomom.Y()) * TMath::Pi();
2825  //double recoX0, recoY0; //[R.K. 01/2017] unused variable?
2826  if(charge > 0) {
2827  xc = recopos.X() + radius * TMath::Cos(beta);
2828  yc = recopos.Y() - radius * TMath::Sin(beta);
2829  }
2830  else {
2831  xc = recopos.X() - radius * TMath::Cos(beta);
2832  yc = recopos.Y() + radius * TMath::Sin(beta);
2833  }
2834 
2835  // vector calculation (alternative): tested, it works!
2836  // TVector2 direction(recomom.X(), recomom.Y());
2837  // direction = direction.Unit();
2838  // TVector2 rad(charge * direction.Y() * R, - charge * direction.X() * R);
2839 
2840  // TVector2 center = recopos.XYvector() + rad;
2841  // xc = center.X();
2842  // yc = center.Y();
2843 
2844 
2845  // ---------------------------------------------------
2846  FairTrackParP recoparlast = sttmvd->GetParamLast();
2847  TVector3 recoposlast = recoparlast.GetPosition();
2848 
2849 
2850  if(fVerbose > 0) {
2851  cout << "charge/xc/yc/radius " << charge << " " << xc << " " << yc << " " << radius << endl;
2852  recomom.Print();
2853  recopos.Print();
2854  recoparlast.GetMomentum().Print();
2855  recoposlast.Print();
2856  }
2857 
2858 
2859  fitm = recomom.Z() / recomom.Perp(); // CHECK fitm = pz / pt :-)GOOD!
2860 
2861  // x0 y0
2862  Double_t d = TMath::Sqrt(xc * xc + yc * yc) - radius;
2863  Double_t phi = TMath::ATan2(yc, xc);
2864 
2865  Double_t x0 = d * TMath::Cos(phi);
2866  Double_t y0 = d * TMath::Sin(phi);
2867 
2868  Double_t Phi0 = TMath::ATan2((y0 - yc),(x0 - xc));
2869  Double_t scosfirst = 0, scoslast = 0.;
2870 
2871  if(fVerbose > 0) cout << "Phi0 " << Phi0 * TMath::RadToDeg() << endl;
2872  // CHECK :-)GOOD! ...
2873  TVector2 v(x0 - xc, y0 - yc);
2874  double alpha1 = TMath::ATan2(recopos.Y() - y0 + radius * TMath::Sin(Phi0), recopos.X() - x0 + radius * TMath::Cos(Phi0));
2875  TVector2 p1(recopos.X() - xc, recopos.Y() - yc);
2876  Double_t Fi1 = CalculatePhi(v, p1, alpha1, Phi0, charge);
2877 
2878  if(fVerbose > 0) {
2879  cout << "alpha1, Fi1 " << alpha1 * TMath::RadToDeg() << " " << Fi1 * TMath::RadToDeg() << endl;
2880  p1.Print();
2881  }
2882 
2883  double alpha2 = TMath::ATan2(recoposlast.Y() - y0 + radius * TMath::Sin(Phi0), recoposlast.X() - x0 + radius * TMath::Cos(Phi0));
2884  TVector2 p2(recoposlast.X() - xc, recoposlast.Y() - yc);
2885  Double_t Fi2 = CalculatePhi(v, p2, alpha2, Phi0, charge);
2886  Fi2 = CompareToPreviousPhi(Fi2, Fi1, charge); // CHECK this!
2887 
2888  if(fVerbose > 0) {
2889  cout << "alpha2, Fi2 " << alpha2 * TMath::RadToDeg() << " " << Fi2 * TMath::RadToDeg() << endl;
2890  p2.Print();
2891  }
2892 
2893  scosfirst = - charge * radius * Fi1; // scos = -q * R * phi CHECK :-)GOOD!
2894  scoslast = - charge * radius * Fi2; // CHECK :-)GOOD!
2895  // ............. :-)GOOD!
2896 
2897  // z = z0 + scos * fitm
2898  fitp = (recopos.Z() + recoposlast.Z() - fitm * (scosfirst + scoslast)) / 2.; // CHECK :-)GOOD!
2899 
2900 
2901  if(fVerbose > 0) {
2902  cout << "positions first/last" << endl;
2903  recopos.Print();
2904  recoposlast.Print();
2905 
2906  cout << "scosfirst/scoslast " << scosfirst << " " << scoslast << endl;
2907  cout << "fitm/fitp " << fitm << " " << fitp << endl;
2908  cout << "z1/z2 " << fitp + fitm * scosfirst << " " << fitp + fitm * scoslast << endl;
2909  }
2910  return true;
2911 }
Double_t CompareToPreviousPhi(Double_t Fi, Double_t Fi_pre, int charge)
Double_t x0
Definition: checkhelixhit.C:70
int fVerbose
Definition: poormantracks.C:24
TObjArray * d
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
static T Sin(const T &x)
Definition: PndCAMath.h:42
Double_t CalculatePhi(TVector2 v, TVector2 p, double alpha, double Phi0, int charge)
static T Cos(const T &x)
Definition: PndCAMath.h:43
__m128 v
Definition: P4_F32vec4.h:4
Double_t
Double_t y0
Definition: checkhelixhit.C:71
static T ATan2(const T &y, const T &x)
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TPad * p2
Definition: hist-t7.C:117
TPad * p1
Definition: hist-t7.C:116
Double_t Pi
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
Int_t PndSttMvdGemTracking::GetPdgFromMC ( int  trackid)

Definition at line 3466 of file PndSttMvdGemTracking.cxx.

References fDefaultPdgCode, fMCTrackArray, fTrackIDArray, PndTrackID::GetCorrTrackID(), PndTrackID::GetNCorrTrackId(), PndMCTrack::GetPdgCode(), and mctrack.

Referenced by GetChargeCorrectedPdgFromMC().

3466  {
3467 
3468  PndTrackID *trackID = (PndTrackID*) fTrackIDArray->At(trackid);
3469  if (trackID->GetNCorrTrackId()>0)
3470  {
3471  Int_t mctrackid = trackID->GetCorrTrackID();
3472  if (mctrackid == -1) {
3473  std::cout << "-W- PndSttMvdGemTracking: no mctrackid - use default hypo " << fDefaultPdgCode << std::endl;
3474  return fDefaultPdgCode;
3475  }
3476 
3477  PndMCTrack *mctrack = (PndMCTrack*) fMCTrackArray->At(mctrackid);
3478  if (!mctrack)
3479  {
3480  std::cout << "-W- PndSttMvdGemTracking: no mctrack " << mctrackid << " - use default hypo " << fDefaultPdgCode << std::endl;
3481  return fDefaultPdgCode;
3482  }
3483 
3484  Int_t pdg = mctrack->GetPdgCode();
3485 
3486  if(pdg >= 100000000) // ion
3487  {
3488  std::cout << "-W- PndSttMvdGemTracking: we have an ion here " << pdg << " - use default hypo " << fDefaultPdgCode << std::endl;
3489  return fDefaultPdgCode;
3490  }
3491  else if (TDatabasePDG::Instance()->GetParticle(mctrack->GetPdgCode())->Charge() == 0)
3492  {
3493  std::cout << "-W- PndSttMvdGemTracking: we have an neutral here " << pdg << " - use default hypo " << fDefaultPdgCode << std::endl;
3494  return fDefaultPdgCode;
3495  }
3496  return pdg;
3497  } // end of "at least one correlated mc index"
3498  else
3499  {
3500  std::cout << "-W- PndSttMvdGemTracking: no correlated mctrackid at all - use default hypo " << fDefaultPdgCode << std::endl;
3501  return fDefaultPdgCode;
3502  }
3503 }
PndMCTrack * mctrack
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
Int_t GetCorrTrackID(Int_t i=0) const
Definition: PndTrackID.h:35
Short_t GetNCorrTrackId(void) const
Definition: PndTrackID.h:34
Bool_t PndPersistencyTask::GetPersistency ( )
inlineinherited

Definition at line 32 of file PndPersistencyTask.h.

References PndPersistencyTask::fPersistency.

Referenced by PndLmdPixelHitProducerFast::GetPersistance(), PndMdtDigitization::Init(), PndMdtHitProducerIdeal::Init(), PndMdtClusterTask::Init(), PndFtsHitProducerRealFast::Init(), PndDiscTaskReconstruction::Init(), PndRichHitProducer::Init(), PndSttHitProducerRealFast::Init(), PndSttHelixHitProducer::Init(), PndDiscTaskPID::Init(), PndIdealTrackFinder::Init(), Init(), PndMdtTrkProducer::Init(), PndFtsHitProducerRealFull::Init(), PndLmdPixelClusterTask::Init(), PndSttHitProducerRealFull::Init(), PndLmdStripClusterTask::Init(), PndEmcApdHitProducer::Init(), PndMissingPzCleanerTask::Init(), PndEmcMakeRecoHit::Init(), PndEmcMakeClusterOnline::Init(), PndTrackSmearTask::Init(), PndEmcFWEndcapTimebasedWaveforms::Init(), PndSttHitProducerIdeal::Init(), PndEmcFWEndcapDigi::Init(), PndFtsHitProducerIdeal::Init(), PndEmcMakeCluster::Init(), PndMdtPointsToWaveform::Init(), PndDiscTaskDigitization::Init(), PndEmcMakeDigi::Init(), PndSdsTimeWalkCorrTask::Init(), PndLmdPixelHitProducerFast::Init(), PndDrcHitFinder::Init(), PndRichHitFinder::Init(), PndEmcMakeCorr::Init(), PndFtofHitProducerIdeal::Init(), PndEmcHitsToWaveform::Init(), PndSciTDigiTask::Init(), PndDrcHitProducerIdeal::Init(), PndSdsHitProducerIdeal::Init(), PndSciTHitProducerIdeal::Init(), PndRecoMultiKalmanTask2::Init(), PndEmcHitProducer::Init(), PndDrcHitProducerReal::Init(), PndDskFLGHitProducerIdeal::Init(), PndEmcTmpWaveformToDigi::Init(), PndDrcDigiTask::Init(), PndEmcWaveformToDigi::Init(), PndSttMatchTracks::Init(), PndEmcWaveformToCalibratedDigi::Init(), PndTrkTracking2::Init(), PndSttFindTracks::Init(), PndEmcMultiWaveformToCalibratedDigi::Init(), PndRecoKalmanTask2::Init(), PndDrcTimeDigiTask::Init(), PndEmcExpClusterSplitter::Init(), PndFtsHoughTrackerTask::Init(), PndSdsNoiseProducer::Init(), PndEmcPhiBumpSplitter::Init(), PndSdsIdealRecoTask::Init(), PndSdsHybridHitProducer::Init(), PndRecoMultiKalmanTask::Init(), PndSdsIdealClusterTask::Init(), PndRecoKalmanTask::Init(), PndSdsStripHitProducerDif::Init(), PndGemDigitize::Init(), PndSdsStripHitProducer::Init(), PndGemFindHits::Init(), PndSdsPixelClusterTask::Init(), PndSdsStripClusterTask::Init(), PndMvdGemTrackFinderOnHits::Init(), PndBarrelTrackFinder::Init(), PndEmcFullDigiTask::PndEmcFullDigiTask(), PndEmcMakeBump::PndEmcMakeBump(), PndUnassignedHitsTask::RegisterBranches(), PndMvdClusterTask::SetPersistance(), PndMvdDigiTask::SetPersistance(), PndEmcMakeBump::SetStorageOfData(), and PndEmcFullDigiTask::StoreDigi().

32 { return fPersistency; }
Int_t PndSttMvdGemTracking::GetPosIndex ( PndGemHit hit)

Definition at line 1374 of file PndSttMvdGemTracking.cxx.

References PndGemHit::GetSensorNr(), PndGemHit::GetStationNr(), and IDEAL.

Referenced by AddRemainingHits(), CheckCombinatorial(), ConsiderCombinatorialEffect(), EvaluatePerformances(), IsAssignable(), OnlyOneHitToEachTrack(), and OrderGemHits().

1374  {
1375 
1376  int istat = hit->GetStationNr();
1377  int isens = hit->GetSensorNr();
1378 
1379  // translate posindex to index in ordering
1380  // 11 ==> 0
1381  // 12 ==> 1
1382  // ...
1383  int posindex = istat * 10 + isens;
1384 
1385  // CHECK delete this, works only for ideal GEM hit prod *********
1386  if(IDEAL == true) {
1387  if(hit->GetZ() < 117) posindex = 11;
1388  else if(hit->GetZ() < 118) posindex = 12;
1389  else if(hit->GetZ() < 153) posindex = 21;
1390  else if(hit->GetZ() < 154) posindex = 22;
1391  else if(hit->GetZ() < 189) posindex = 31;
1392  else if(hit->GetZ() < 190) posindex = 32;
1393  }
1394  // *****************************************************************
1395 
1396  return posindex;
1397 }
#define IDEAL
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
Int_t GetStationNr() const
Definition: PndGemHit.h:81
Int_t PndSttMvdGemTracking::GetTrackIndex ( Int_t  i)

Definition at line 1221 of file PndSttMvdGemTracking.cxx.

References trackindexes.

Referenced by AddRemainingHits(), CheckCombinatorial(), Exec(), OnlyOneHitToEachTrack(), and Retrack().

1221  {
1222 
1223  // probably this was wrong
1224 // int counter = -1;
1225 // int tmptrack = -1;
1226 // std::vector< std::pair<int, int> >::iterator iter;
1227 // for(iter = trackvector.begin(); iter != trackvector.end(); iter++) {
1228 // std::pair<int, int> thispair = *iter;
1229 // if(thispair.first != tmptrack) {
1230 // tmptrack = thispair.first;
1231 // counter++;
1232 // cout << i << " " << counter << " " << tmptrack << " " << trackvector[i].first << endl;
1233 // if(counter == i) { cout << "TRACK INDEX " << tmptrack << endl; return tmptrack; }
1234 // }
1235 // }
1236 
1237  if(i < (int)trackindexes.size()) return trackindexes[i];
1238  cout << "PndSttMvdGemTracking::GetTrackIndex " << i << " Out Of Bounds" << endl;
1239  return -1;
1240 }
Int_t i
Definition: run_full.C:25
std::vector< int > trackindexes
std::vector< int > PndSttMvdGemTracking::GetTracksAssociatedToHit ( Int_t  ihit)

Definition at line 1318 of file PndSttMvdGemTracking.cxx.

References fVerbose, and trackvector.

Referenced by AddRemainingHits(), EvaluatePerformances(), and ForbidMultiAssignedHits().

1318  {
1319  // which track/s does hit # belong to
1320  std::vector<int> thishittracks;
1321  std::vector< std::pair<int, int> >::iterator iter;
1322  for(iter = trackvector.begin(); iter != trackvector.end(); iter++) {
1323  std::pair<int, int> thispair = *iter;
1324  if(thispair.second == ihit) thishittracks.push_back(thispair.first);
1325  }
1326  if(fVerbose != 0) {
1327  std::cout << "hit " << ihit << " belongs to " << thishittracks.size() << " tracks: " ;
1328  for(size_t j = 0; j < thishittracks.size(); j++) std::cout << thishittracks[j] << " ";
1329  std::cout << std::endl;
1330  }
1331  return thishittracks;
1332 }
int fVerbose
Definition: poormantracks.C:24
std::vector< std::pair< int, int > > trackvector
InitStatus PndSttMvdGemTracking::Init ( )
virtual

Virtual method Init

Definition at line 121 of file PndSttMvdGemTracking.cxx.

References countbad, countdoubt, countgood, countnohitonplane, countplane, countprinotassigned, countprop, countreconstructablehit, countsecnotassigned, countsttmvd, countsttmvdusable, evt, fCompleteTrackArray, fCompleteTrackCandArray, fDefaultPdgCode, fGemBranchName, fGemHitArray, fGemPointArray, PndSttMapCreator::FillTubeArray(), fMCTrackArray, fMvdPixelBranchName, fMvdPixelHitArray, fMvdPointArray, fMvdStripBranchName, fMvdStripHitArray, fPdgFromMC, fPro, fStartTrackBranchName, fStartTrackCandBranchName, fStartTrackIDBranchName, fSttBranchName, fSttHitArray, fSttParameters, fSttPointArray, fTrackArray, fTrackCandArray, fTrackIDArray, fTubeArray, fVerbose, PndPersistencyTask::GetPersistency(), and SetupGEMPlanes().

121  {
122 
123  evt = 0;
124  countgood = 0; // correct association
125  countbad = 0; // wrong association
126  countdoubt = 0; // doubt in ref index CHECK (temporary)
127  countprinotassigned = 0; // not assigned hits from primary tracks
128  countsecnotassigned = 0; // not assigned hits from secondary tracks
129  countreconstructablehit = 0; // n of hits which belong to a mctrack which has mvd + stt candidate
130  for(int iplane = 0; iplane < 8; iplane++) {
131  countplane[iplane] = 0;
132  countprop[0][iplane] = 0;
133  countprop[1][iplane] = 0;
134  countnohitonplane[iplane] = 0;
135 
136  }
137  countsttmvd = 0;
138  countsttmvdusable = 0;
139 
140  FairRootManager *ioman = FairRootManager::Instance();
141 
142  cout << "-I- -------------------" << endl;
143  cout << "-I- PndSttMvdGemTracking: using branches "
144  << fMvdPixelBranchName << " "
145  << fMvdStripBranchName << " "
146  << fSttBranchName << " "
147  << fGemBranchName << endl;
148  cout << "-I- to change one or more of these use PndSttMvdGemTracking:SetBranchName( TStrings ); the order of TStrings is mvd pixel name, mvd strip name, stt name, gem name" << endl;
149  cout << "starting track for extrapolation " << fStartTrackBranchName << " " << fStartTrackCandBranchName << endl;
150  if(fPdgFromMC) cout << "starting trackid for extrapolation " << fStartTrackIDBranchName << endl;
151  cout << "-I- -------------------" << endl;
152 
153 
154  if (!ioman)
155  {
156  cout << "-E- PndSttMvdGemTracking: "
157  << "RootManager not instantised!" << endl;
158  return kFATAL;
159  }
160 
161  // open SttMvdTrackCand array
162  fTrackCandArray=(TClonesArray*) ioman->GetObject(fStartTrackCandBranchName);
163  if(fTrackCandArray==0){
164  cout << "fStartTrackCandBranchName " << fStartTrackCandBranchName << " not found" << endl;
165  Error("PndSttMvdGemTracking:Init","stt + mvd trackcand - array not found!");
166  return kERROR;
167  }
168 
169  // open SttMvdTrack array
170  fTrackArray = (TClonesArray*) ioman->GetObject(fStartTrackBranchName);
171  if(!fTrackArray) {
172  Error("PndSttMvdGemTracking:Init","stt + mvd track - array not found!");
173  return kERROR;
174  }
175 
176  // open SttMvdTrack array
177  if(fPdgFromMC) {
178  fTrackIDArray = (TClonesArray*) ioman->GetObject(fStartTrackIDBranchName);
179  if(!fTrackIDArray) {
180  Error("PndSttMvdGemTracking:Init","stt + mvd trackID - array not found!");
181  return kERROR;
182  }
183  }
184 
185  // open GEM hit array
186  fGemHitArray = (TClonesArray*) ioman->GetObject(fGemBranchName);
187  if(!fGemHitArray) {
188  Error("PndSttMvdGemTracking:Init","gem hit - array not found!");
189  return kERROR;
190  }
191 
192  // open GEM point array
193  fGemPointArray = (TClonesArray*) ioman->GetObject("GEMPoint");
194  if(!fGemPointArray) {
195  Error("PndSttMvdGemTracking:Init","gem point - array not found!");
196  return kERROR;
197  }
198 
199  // open MC track array CHECK to be deleted for real data
200  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
201  if(!fMCTrackArray) {
202  Error("PndSttMvdGemTracking:Init","mctrack - array not found!");
203  return kERROR;
204  }
205 
206  // open MVD pixel hit array
207  fMvdPixelHitArray = (TClonesArray*) ioman->GetObject(fMvdPixelBranchName);
208  if(!fMvdPixelHitArray) {
209  Error("PndSttMvdGemTracking:Init","mvd pixel hit - array not found!");
210  return kERROR;
211  }
212 
213  // open MVD strip hit array
214  fMvdStripHitArray = (TClonesArray*) ioman->GetObject(fMvdStripBranchName);
215  if(!fMvdStripHitArray) {
216  Error("PndSttMvdGemTracking:Init","mvd strip hit - array not found!");
217  return kERROR;
218  }
219 
220  // open STT hit array
221  fSttHitArray = (TClonesArray*) ioman->GetObject(fSttBranchName);
222  if(!fSttHitArray) {
223  Error("PndSttMvdGemTracking:Init","stt hit - array not found!");
224  return kERROR;
225  }
226 
227  // open MVD point array
228  fMvdPointArray = (TClonesArray*) ioman->GetObject("MVDPoint");
229  if(!fMvdPointArray) {
230  Error("PndSttMvdGemTracking:Init","mvd point - array not found!");
231  return kERROR;
232  }
233 
234  // open STT point array
235  fSttPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
236  if(!fSttPointArray) {
237  Error("PndSttMvdGemTracking:Init","stt point - array not found!");
238  return kERROR;
239  }
240 
241  // Create and register PndTrack array
242  fCompleteTrackCandArray = new TClonesArray("PndTrackCand", 100);
243  ioman->Register("SttMvdGemTrackCand", "SttMvdGem", fCompleteTrackCandArray, GetPersistency());
244 
245  fCompleteTrackArray = new TClonesArray("PndTrack", 100);
246  ioman->Register("SttMvdGemTrack", "SttMvdGem", fCompleteTrackArray, GetPersistency());
247 
248  // GEANE propagation to volume
249  fPro = new FairGeanePro();
250  if (fVerbose==0) fPro->SetPrintErrors(kFALSE);
251 
252 
253  // STT mapper
255  fTubeArray = mapper->FillTubeArray();
256 
257  // set up geometry of GEMs;
258  SetupGEMPlanes();
259 
260  if(fPdgFromMC) cout << "-I- PndSttMvdGemTracking: Taking PDG from MC (keeping backup default opt = " << fDefaultPdgCode << ")" << endl;
261  else cout << "-I- PndSttMvdGemTracking: using default PDG " << fDefaultPdgCode << endl;
262 
263  cout << "-I- PndSttMvdGemTracking: Intialisation successfull" << endl;
264  return kSUCCESS;
265 }
int fVerbose
Definition: poormantracks.C:24
TClonesArray * fTrackCandArray
TClonesArray * fCompleteTrackCandArray
TClonesArray * fCompleteTrackArray
TClonesArray * fMvdPixelHitArray
TClonesArray * fMvdStripHitArray
TClonesArray * FillTubeArray()
Bool_t PndSttMvdGemTracking::IntersectionFinder ( Double_t  xc,
Double_t  yc,
Double_t  radius,
PndSttHit stthit,
TVector3 &  xyz,
TVector3 &  dxyz 
)

Definition at line 2452 of file PndSttMvdGemTracking.cxx.

References Double_t, fabs(), fTubeArray, fVerbose, PndSttHit::GetIsochrone(), PndSttHit::GetIsochroneError(), PndSttTube::GetPosition(), PndSttHit::GetTubeID(), PndSttTube::GetWireDirection(), m, point, and sqrt().

Referenced by Prefit().

2453  {
2454 
2455  // tubeID CHECK added
2456  Int_t tubeID = stthit->GetTubeID();
2457  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
2458  TVector3 wiredirection = tube->GetWireDirection();
2459 
2460  if(wiredirection != TVector3(0.,0.,1.)) {
2461  // cout << "wire skewed" << endl;
2462  return false;
2463  }
2464 
2465  // [xp, yp] point = coordinates xy of the centre of the firing tube
2466  TVector2 point;
2467  point.Set(tube->GetPosition().X(), tube->GetPosition().Y());
2468  Double_t isochrone = stthit->GetIsochrone();
2469 
2470  // the coordinates of the point are taken from the intersection
2471  // between the circumference from the drift time and the R radius of
2472  // curvature. -------------------------------------------------------
2473  // 2. find the intersection between the little circle and the line // R
2474  // 2.a
2475  // find the line passing throught [xc, yc] (centre of curvature) and [xp, yp] (first wire)
2476  // y = mx + q
2477  Double_t m = (point.Y() - yc)/(point.X() - xc);
2478  Double_t q = point.Y() - m*point.X();
2479 
2480  // cut on radius CHECK
2481  // if the simulated radius is too small, the stthit
2482  // is not used because rouning errors may occur
2483  // if(isochrone < 0.7e-3) { cout << "isochrone < 0.7e-3" << endl; return false; } // CHECK
2484 
2485  Double_t x1 = 0, y1 = 0,
2486  x2 = 0, y2 = 0,
2487  xb1 = 0, yb1 = 0,
2488  xb2 = 0, yb2 = 0;
2489 
2490  // CHECK the vertical track
2491  if(fabs(point.X() - xc) < 1e-6) {
2492 
2493  // 2.b
2494  // intersection little circle and line --> [x1, y1]
2495  // + and - refer to the 2 possible intersections
2496  // +
2497  x1 = point.X();
2498  y1 = point.Y() + sqrt(isochrone * isochrone - (x1 - point.X()) * (x1 - point.X()));
2499  // -
2500  x2 = x1;
2501  y2 = point.Y() - sqrt(isochrone * isochrone - (x2 - point.X()) * (x2 - point.X()));
2502 
2503  // 2.c intersection between line and circle
2504  // +
2505  xb1 = xc;
2506  yb1 = yc + sqrt(radius * radius - (xb1 - xc) * (xb1 - xc));
2507  // -
2508  xb2 = xb1;
2509  yb2 = yc - sqrt(radius * radius - (xb2 - xc) * (xb2 - xc));
2510 
2511  } // END CHECK
2512  else {
2513 
2514  // 2.b
2515  // intersection little circle and line --> [x1, y1]
2516  // + and - refer to the 2 possible intersections
2517  // +
2518 
2519  if(((m*(q - point.Y()) - point.X())*(m*(q - point.Y()) - point.X()) - (m*m + 1)*((q - point.Y())*(q - point.Y()) + point.X()*point.X() - isochrone*isochrone)) < 0.) { if(fVerbose > 0) cout << "IntersectionFinder round errors: " << isochrone << endl;
2520  return false; }
2521 
2522  x1 = (-(m*(q - point.Y()) - point.X()) + sqrt((m*(q - point.Y()) - point.X())*(m*(q - point.Y()) - point.X()) - (m*m + 1)*((q - point.Y())*(q - point.Y()) + point.X()*point.X() - isochrone*isochrone))) / (m*m + 1);
2523  y1 = m*x1 + q;
2524  // -
2525  x2 = (-(m*(q - point.Y()) - point.X()) - sqrt((m*(q - point.Y()) - point.X())*(m*(q - point.Y()) - point.X()) - (m*m + 1)*((q - point.Y())*(q - point.Y()) + point.X()*point.X() - isochrone*isochrone))) / (m*m + 1);
2526  y2 = m*x2 + q;
2527 
2528  // 2.c intersection between line and circle
2529  // +
2530  xb1 = (-(m*(q - yc) - xc) + sqrt((m*(q - yc) - xc)*(m*(q - yc) - xc) - (m*m + 1)*((q - yc)*(q - yc) + xc*xc - (radius) *(radius) ))) / (m*m + 1);
2531  yb1 = m*xb1 + q;
2532  // -
2533  xb2 = (-(m*(q - yc) - xc) - sqrt((m*(q - yc) - xc)*(m*(q - yc) - xc) - (m*m + 1)*((q - yc)*(q - yc) + xc*xc - (radius) *(radius)))) / (m*m + 1);
2534  yb2 = m*xb2 + q;
2535  }
2536 
2537  // calculation of the distance between [xb, yb] and [xp, yp]
2538  Double_t distb1 = sqrt((yb1 - point.Y())*(yb1 - point.Y()) + (xb1 - point.X())*(xb1 - point.X()));
2539  Double_t distb2 = sqrt((yb2 - point.Y())*(yb2 - point.Y()) + (xb2 - point.X())*(xb2 - point.X()));
2540 
2541  // choice of [xb, yb]
2542  TVector2 xyb;
2543  if(distb1 > distb2) xyb.Set(xb2, yb2);
2544  else xyb.Set(xb1, yb1);
2545 
2546  // calculation of the distance between [x, y] and [xb. yb]
2547  Double_t dist1 = sqrt((xyb.Y() - y1)*(xyb.Y() - y1) + (xyb.X() - x1)*(xyb.X() - x1));
2548  Double_t dist2 = sqrt((xyb.Y() - y2)*(xyb.Y() - y2) + (xyb.X() - x2)*(xyb.X() - x2));
2549 
2550  // choice of [x, y]
2551  if(dist1 > dist2) xyz.SetXYZ(x2, y2, stthit->GetZ());
2552  else xyz.SetXYZ(x1, y1, stthit->GetZ()); // <========= THIS IS THE NEW POINT to be used for the fit
2553 
2554  Double_t sigr = stthit->GetIsochroneError();
2555  Double_t sigx = sigr; // fabs(sigr * TMath::Cos(m));
2556  Double_t sigy = sigr; // fabs(sigr * TMath::Sin(m));
2557  dxyz.SetXYZ(sigx, sigy, 0);
2558 
2559  return kTRUE;
2560 }
int fVerbose
Definition: poormantracks.C:24
__m128 m
Definition: P4_F32vec4.h:28
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
Int_t GetTubeID() const
Definition: PndSttHit.h:75
Double_t GetIsochroneError() const
Definition: PndSttHit.h:63
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
Double_t PndSttMvdGemTracking::IsAssignable ( FairTrackParP *  gempar,
PndGemHit gemhit 
)

Definition at line 1675 of file PndSttMvdGemTracking.cxx.

References CAMath::ATan2(), CAMath::Cos(), display, Double_t, dPhi, fDisplayOn, fMaxDistance, fOrdering, fOrderingIterator, fTimes, fTurn, fVerbose, PndGemHit::GetDp(), PndGemHit::GetDr(), GetPosIndex(), PndGemHit::GetPosition(), hchosen, hchosen2, hdist, hdist2, hsigma, hsigma2, IDEAL, phi, CAMath::Sin(), and CAMath::Sqrt().

Referenced by AssignHits().

1675  {
1676 
1677  TVector3 gemErr(gempar->GetDV(),
1678  gempar->GetDW(),
1679  0.0); // gempar->GetDZ());
1680 
1681  // error of gemhit
1682  // dr = radial
1683  // dp = perpendicular to radius
1684  // dp ~ dl = dphi * R --> dphi ~ dp / R
1685  Double_t phi = TMath::ATan2(gemhit->GetY(), gemhit->GetX());
1686  Double_t sinp = TMath::Sin(phi);
1687  Double_t cosp = TMath::Cos(phi);
1688  //Double_t raddist = gemhit->GetPosition().Perp(); //[R.K. 01/2017] unused variable?
1689  Double_t dR, dP;
1690  if(IDEAL == true) {
1691  dR = 0.01;
1692  dP = 0.08;
1693  }
1694  else {
1695  dR = gemhit->GetDr();
1696  dP = gemhit->GetDp();
1697  }
1698  Double_t dPhi = dP; // / raddist;
1699  Double_t errx = TMath::Sqrt(cosp * cosp * dR * dR + sinp * sinp * dPhi * dPhi);
1700  Double_t erry = TMath::Sqrt(sinp * sinp * dR * dR + cosp * cosp * dPhi * dPhi);
1701 
1702 
1703  TVector3 hitErr(errx, erry, 0.0); // CHECK
1704  if(fVerbose > 0) {
1705  cout << "errors " << endl;
1706  gemErr.Print();
1707  hitErr.Print();
1708  }
1709  TVector3 distVec = gempar->GetPosition() - gemhit->GetPosition();
1710  double distance = distVec.Perp(); // CHECK using xy distance ON plane
1711  // error on distance: sqrt((x - x0)**2 * (sigx**2 + sigx0**2) +
1712  // ... (y) ...) / distance
1713  // with x, y, z = extrap pos/err
1714  // x0, y0, z0 = gem hit pos/err
1715  double distErr = (1./distance) * TMath::Sqrt((distVec.X() * distVec.X()) *
1716  (gemErr.X() * gemErr.X() +
1717  hitErr.X() * hitErr.X()) +
1718  (distVec.Y() * distVec.Y()) *
1719  (gemErr.Y() * gemErr.Y() +
1720  hitErr.Y() * hitErr.Y()) ); // CHECK if distance == 0
1721 
1722  if(fVerbose > 0) cout << "distance " << distance << " err " << distErr << endl;
1723 
1724  // maximum distance **************************
1725  double maxdistance = 0;
1726  if(fMaxDistance != -1) maxdistance = fMaxDistance;
1727  else {
1728  maxdistance = fTimes * distErr; // distErr; // CHECK if it has to be distErr ; // CHECK put this as a parameter and tune it
1729  if(maxdistance < 1) maxdistance = 1;
1730  if(maxdistance > 10) maxdistance = 10;
1731  }
1732  if(fTurn == 2 && gempar->GetMomentum().Mag() >= 0.5) maxdistance = 5; // CHECK NOW
1733 
1734 
1735  // cout << "-> " << fTimes << " " << fTurn << " " << maxdistance << endl;
1736 
1737  int posindex = GetPosIndex(gemhit);
1738  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
1739  int ipos = fOrderingIterator - fOrdering.begin();
1740  if(fTurn == 1) {
1741  hdist[ipos]->Fill(distance);
1742  hsigma[ipos]->Fill(maxdistance);
1743  }
1744  else if(fTurn == 2) {
1745  hdist2[ipos]->Fill(distance);
1746  hsigma2[ipos]->Fill(maxdistance);
1747  }
1748 
1749  if(fDisplayOn) {
1750  // int posindex = GetPosIndex(gemhit);
1751 // fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
1752 // int ipos = fOrderingIterator - fOrdering.begin();
1753 
1754  TArc *earc = new TArc(gempar->GetPosition().X(), gempar->GetPosition().Y(), maxdistance);
1755  earc->SetFillStyle(0);
1756  if(fTurn == 1) earc->SetLineColor(3); // green
1757  else if(fTurn == 2) earc->SetLineColor(2); // red
1758  display->cd(ipos + 1);
1759  earc->Draw("SAME");
1760  display->Update();
1761  display->Modified();
1762  }
1763 
1764 
1765 
1766 
1767  if(distance <= maxdistance) {
1768  if(fTurn == 1) hchosen[ipos]->Fill(distance);
1769  else if(fTurn == 2) hchosen2[ipos]->Fill(distance);
1770  return distance;
1771  }
1772  return -1;
1773 
1774 }
std::vector< int >::iterator fOrderingIterator
int fVerbose
Definition: poormantracks.C:24
Int_t GetPosIndex(PndGemHit *hit)
#define IDEAL
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
static T Sin(const T &x)
Definition: PndCAMath.h:42
static T Cos(const T &x)
Definition: PndCAMath.h:43
std::vector< int > fOrdering
Double_t
static T ATan2(const T &y, const T &x)
Double_t GetDp() const
Definition: PndGemHit.h:76
TVector3 GetPosition() const
Definition: PndGemHit.h:71
Double_t GetDr() const
Definition: PndGemHit.h:75
double dPhi
Definition: RiemannTest.C:17
void PndSttMvdGemTracking::Kalman ( TMatrixT< double >  extrap,
TMatrixT< double >  measurement,
TMatrixT< double >  H,
TMatrixT< double >  extrap_cov,
TMatrixT< double >  measurement_cov,
TMatrixT< double > &  kalman,
TMatrixT< double > &  kalman_cov 
)

Definition at line 1832 of file PndSttMvdGemTracking.cxx.

Referenced by Retrack().

1836  {
1837 
1838 
1839  // gain
1840  // calculate sum of covariances = sum_cov = m_cov + H e_cov HT
1841  TMatrixT<double> extrapHT_cov(extrap_cov, TMatrixT<double>::kMultTranspose, H);
1842  TMatrixT<double> HextrapHT_cov(H, TMatrixT<double>::kMult, extrapHT_cov);
1843  TMatrixT<double> sum_cov = measurement_cov + HextrapHT_cov;
1844  //double det = 0; // CHECK fill it //[R.K. 01/2017] unused variable?
1845  sum_cov.Invert();
1846  // calculate gain
1847  TMatrixT<double> Hsum_cov(H, TMatrixT<double>::kTransposeMult, sum_cov);
1848  TMatrixT<double> gain(extrap_cov, TMatrixT<double>::kMult, Hsum_cov);
1849 
1850  // kalman state
1851  // extrap + gain * (measurement - H extrap)
1852  TMatrixT<double> Hextrap(H, TMatrixT<double>::kMult, extrap);
1853  TMatrixT<double> meas_Hextrap = measurement - Hextrap;
1854  TMatrixT<double> gainmeas_Hextrap(gain, TMatrixT<double>::kMult, meas_Hextrap);
1855  kalman = extrap + gainmeas_Hextrap;
1856 
1857  // kalman cov matrix
1858  // kalman_cov = (identity - gain H) extrap_cov
1859  TMatrixT<double> gainH(gain, TMatrixT<double>::kMult, H);
1860  TMatrixT<double> I(5, 5);
1861  for(int imat = 0; imat < 5; imat++) {
1862  for(int jmat = 0; jmat < 5; jmat++) {
1863  I[imat][jmat] = 0;
1864  if(imat == jmat) I[imat][jmat] = 1;
1865  }
1866  }
1867  TMatrixT<double> I_gainH = I - gainH;
1868  kalman_cov = TMatrixT<double>(I_gainH, TMatrixT<double>::kMult, extrap_cov);
1869 
1870  // kalman
1871  // cout << "kalman" << endl;
1872  // kalman.Print();
1873  // kalman_cov.Print();
1874 }
KalmanTask * kalman
void PndSttMvdGemTracking::OnlyOneHitToEachTrack ( Int_t  nhits,
Int_t  ntracks 
)

Definition at line 1162 of file PndSttMvdGemTracking.cxx.

References CountTracks(), DeleteHitFromTrack(), distancemap, fGemHitArray, fNPositions, fOrdering, fOrderingIterator, fProTracks, fVerbose, GetHitsAssociatedToTrack(), GetPosIndex(), GetTrackIndex(), and i.

Referenced by Exec().

1163 {
1164  if(fVerbose > 0) cout << "ONLY ONE HIT FOR EACH TRACK : DELETING HITS" << endl;
1165 
1166  for(Int_t i = 0; i < CountTracks(); i++) {
1167  int itrk = GetTrackIndex(i);
1168  if(fProTracks[itrk] == false) continue;
1169  std::vector<int> thistrackhits = GetHitsAssociatedToTrack(itrk);
1170 
1171  std::vector<double> tmpposdistance; // CHECK needs init?
1172  std::vector<int> tmpposhitindex; // CHECK needs init?
1173 
1174  for(Int_t ipos = 0; ipos < fNPositions; ipos++) {
1175  tmpposdistance.push_back(-1);
1176  tmpposhitindex.push_back(-1);
1177  }
1178 
1179  for(size_t j = 0; j < thistrackhits.size(); j++) {
1180  int ihit = thistrackhits[j];
1181 
1182  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(ihit);
1183  if(!gemhit) continue;
1184 
1185  int posindex = GetPosIndex(gemhit);
1186 
1187  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
1188  int index = fOrderingIterator - fOrdering.begin();
1189 
1190  if(tmpposdistance[index] == -1){
1191  tmpposdistance[index] = distancemap[itrk][ihit]; // CHECK
1192  tmpposhitindex[index] = ihit;
1193  }
1194  else {
1195  // if more than one hit on the same plane
1196  // is assigned to the same track, choose!
1197  if(fVerbose > 0){
1198  cout << "hit " << ihit << " and " << tmpposhitindex[index] << " on the same plane" << endl;
1199  cout << "distance " << ihit << " " << distancemap[itrk][ihit] << endl;
1200  cout << "distance " << tmpposhitindex[index] << " " << tmpposdistance[index] << endl;
1201  }
1202  // if this distance is better...
1203  if(distancemap[itrk][ihit] < tmpposdistance[index]) {
1204  // ... delete the old hit and...
1205  DeleteHitFromTrack(tmpposhitindex[index], itrk); // CHECK
1206  // ... update
1207  tmpposdistance[index] = distancemap[itrk][ihit]; // CHECK
1208  tmpposhitindex[index] = ihit;
1209  }
1210  else {
1211  // otherwise, delete this hit!
1212  DeleteHitFromTrack(ihit, itrk);
1213  }
1214  if(fVerbose > 0) cout << tmpposhitindex[index] << " wins" << endl;
1215 
1216  }
1217  }
1218  }
1219 }
std::vector< int >::iterator fOrderingIterator
int fVerbose
Definition: poormantracks.C:24
Int_t GetPosIndex(PndGemHit *hit)
Int_t i
Definition: run_full.C:25
std::vector< int > GetHitsAssociatedToTrack(Int_t itrk)
std::vector< int > fOrdering
TMatrixT< double > distancemap
std::map< int, bool > fProTracks
void DeleteHitFromTrack(Int_t ihit, Int_t itrk)
void PndSttMvdGemTracking::OrderGemHits ( Int_t  nhits)

Definition at line 431 of file PndSttMvdGemTracking.cxx.

References display, fDisplayOn, fGemHitArray, fNPositions, fOrdering, fOrderingIterator, fVerbose, GetPosIndex(), PndGemHit::GetSensorNr(), PndGemHit::GetStationNr(), hit, hitcounter, hitmap, idx, nhits, and towhichplane.

Referenced by Exec().

431  {
432  if(nhits == 0) return;
433 
434  // loop on GEM hits to fill hitcounter
435  // and hitmap and be able to order them
436  for(int ihit = 0; ihit < nhits; ihit++) {
437  PndGemHit *hit = (PndGemHit*) fGemHitArray->At(ihit);
438  if(!hit) continue;
439  if(fVerbose > 0) cout << "hit "
440  << ihit << " @ "
441  << hit->GetX() << " "
442  << hit->GetY() << " "
443  << hit->GetZ() << " "
444  << hit->GetStationNr() << " "
445  << hit->GetSensorNr() << " "
446  << endl;
447 
448  // translate posindex to index in ordering
449  // 11 ==> 0
450  // 12 ==> 1
451  // ...
452  int posindex = GetPosIndex(hit);
453 
454 
455 
456  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
457  int index = fOrderingIterator - fOrdering.begin();
458  int idx = (int) hitcounter(index, 0);
459  if(fVerbose > 0) cout << posindex << " " << idx << " " << index << endl;
460 
461 
462  if(fDisplayOn) {
463  TMarker *gmrk = new TMarker(hit->GetX(), hit->GetY(), 7);
464  display->cd(index + 1);
465  gmrk->Draw("SAME");
466  display->Update();
467  display->Modified();
468  }
469 
470 
471 
472 
473  hitmap(index, idx) = ihit; // CHECK how to make a TMatrixT of int
474  hitcounter(index, 0)++;
475  // set the corresponding position in towhichplane to true
476  towhichplane[index] = true;
477  }
478  // hitmap.Print();
479  // hitcounter.Print();
480 
481  if(fVerbose > 0) for(int ipos = 0; ipos < fNPositions; ipos++) cout << "pos " << ipos << " " << towhichplane[ipos] << endl; // CHECK delete it
482 
483 }
std::vector< int >::iterator fOrderingIterator
int fVerbose
Definition: poormantracks.C:24
TMatrixT< float > hitmap
Int_t GetPosIndex(PndGemHit *hit)
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
std::vector< int > fOrdering
int idx[MAX]
Definition: autocutx.C:38
TMatrixT< float > hitcounter
Int_t GetStationNr() const
Definition: PndGemHit.h:81
PndSdsMCPoint * hit
Definition: anasim.C:70
std::map< int, bool > towhichplane
Bool_t PndSttMvdGemTracking::Prefit ( PndTrack sttmvdTrack,
PndTrackCand sttmvdCand,
TVector3 &  lastpos,
TVector3 &  lastmom 
)

Definition at line 2253 of file PndSttMvdGemTracking.cxx.

References CAMath::ATan2(), Bool_t, CalculatePhi(), CompareToPreviousPhi(), CAMath::Cos(), d, Double_t, fGemBranchName, Fit(), fMvdPixelBranchName, fMvdPixelHitArray, fMvdStripBranchName, fMvdStripHitArray, fSttBranchName, fSttHitArray, fTubeArray, PndSdsHit::GetCov(), PndTrackCandHit::GetDetId(), PndTrackCandHit::GetHitId(), GetInitialParams(), PndSttHit::GetIsochrone(), PndSttHit::GetIsochroneError(), PndTrackCand::GetNHits(), PndTrack::GetParamFirst(), PndSttTube::GetPosition(), PndTrackCand::GetSortedHit(), PndSttHit::GetTubeID(), PndSttTube::GetWireDirection(), IntersectionFinder(), nhits, p1, p2, phi, pz, CAMath::Sin(), CAMath::Sqrt(), v, x0, y0, z, and ZFit().

Referenced by SetStartParameters().

2254 {
2255 
2256  Int_t nhits = sttmvdCand->GetNHits();
2257  // 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
2258  // hitId | detid | x | y | z | dx | dy | dz | iso | isoerr | useforfit (0 = only xy, 1 = only z, 2 = both)
2259  TMatrixT<double> points(nhits, 11);
2260  Double_t xc, yc, radius, fitm, fitp;
2261 
2262  GetInitialParams(sttmvdTrack, xc, yc, radius, fitm, fitp);
2263  // cout << " PR " << xc << " " << yc << " " << radius << " " << fitm << " " << fitp << endl;
2264 
2265  Int_t charge = sttmvdTrack->GetParamFirst().GetQ();
2266  Int_t firsthitid = -1, lasthitid = -1;
2267  for(int ihit = 0; ihit < nhits; ihit++)
2268  {
2269  // get hit
2270  PndTrackCandHit candhit = sttmvdCand->GetSortedHit(ihit);
2271  Int_t hitId = candhit.GetHitId();
2272  Int_t detId = candhit.GetDetId();
2273 
2274  points[ihit][0] = -1;
2275  points[ihit][1] = detId;
2276  points[ihit][2] = 0;
2277  points[ihit][3] = 0;
2278  points[ihit][4] = 0;
2279  points[ihit][5] = 0;
2280  points[ihit][6] = 0;
2281  points[ihit][7] = 0;
2282  points[ihit][8] = 0;
2283  points[ihit][9] = 0;
2284  points[ihit][10] = -1;
2285 
2286  if(detId == FairRootManager::Instance()->GetBranchId(fGemBranchName)) continue;
2287  else if(detId == FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName) || detId == FairRootManager::Instance()->GetBranchId(fMvdStripBranchName)) {
2288  PndSdsHit *mvdhit;
2289  if(detId == FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName)) mvdhit = (PndSdsHit*) fMvdPixelHitArray->At(hitId);
2290  else mvdhit = (PndSdsHit*) fMvdStripHitArray->At(hitId);
2291  if(!mvdhit) { cout << "mvd hit " << ihit << " " << hitId << " does not exist" << endl; continue;}
2292  points[ihit][0] = hitId;
2293  points[ihit][2] = mvdhit->GetX();
2294  points[ihit][3] = mvdhit->GetY();
2295  points[ihit][4] = mvdhit->GetZ();
2296  points[ihit][5] = TMath::Sqrt(mvdhit->GetCov()[0][0]);
2297  points[ihit][6] = TMath::Sqrt(mvdhit->GetCov()[1][1]);
2298  points[ihit][7] = TMath::Sqrt(mvdhit->GetCov()[2][2]);
2299  points[ihit][10] = 2;
2300 
2301  lasthitid = ihit;
2302  if(firsthitid == -1) firsthitid = ihit;
2303  // cout << "MVD " << ihit << " " << hitId << " " << points[ihit][2] << " " << points[ihit][3] << " " << points[ihit][4] << endl;
2304 
2305  }
2306  else if(detId == FairRootManager::Instance()->GetBranchId(fSttBranchName)) {
2307  PndSttHit *stthit = (PndSttHit*) fSttHitArray->At(hitId);
2308  if(!stthit) { cout << "stt hit " << ihit << " " << hitId << " does not exist" << endl; continue; }
2309  TVector3 xyz(0, 0, 0);
2310  TVector3 dxyz(0, 0, 0);
2311 
2312  Int_t tubeID = stthit->GetTubeID();
2313  PndSttTube *tube = (PndSttTube* ) fTubeArray->At(tubeID);
2314 
2315  TVector3 wireDirection = tube->GetWireDirection();
2316  if(wireDirection == TVector3(0., 0., 1.)) { // is parallel
2317  Bool_t intfin = IntersectionFinder(xc, yc, radius, stthit, xyz, dxyz);
2318  if(intfin == false) {
2319  // cout << "intfin false" << endl;
2320  continue;
2321  }
2322 
2323  // CHECK get only parallel tubes?
2324  lasthitid = ihit;
2325  if(firsthitid == -1) firsthitid = ihit;
2326  // cout << "STT " << ihit << " " << hitId << " " << points[ihit][2] << " " << points[ihit][3] << " " << points[ihit][4] << "iso " << points[ihit][8] << endl;
2327 
2328  points[ihit][10] = 0;
2329  }
2330  else { // is skewed
2331  // continue;
2332  xyz = tube->GetPosition();
2333  points[ihit][10] = 1;
2334  }
2335 
2336  points[ihit][0] = hitId;
2337  points[ihit][2] = xyz.X();
2338  points[ihit][3] = xyz.Y();
2339  points[ihit][4] = xyz.Z();
2340  points[ihit][5] = dxyz.X();
2341  points[ihit][6] = dxyz.Y();
2342  points[ihit][7] = dxyz.Z();
2343  points[ihit][8] = stthit->GetIsochrone();
2344  points[ihit][9] = stthit->GetIsochroneError();
2345 
2346  }
2347  }
2348 
2349 
2350  for(int iteration = 0; iteration < 1; iteration++) {
2351  Bool_t fitting = Fit(points, xc, yc, radius);
2352  if(fitting == false) {
2353  // cout << "fit false" << endl;
2354  break;
2355  }
2356  // cout << "FIT " << iteration << " " << xc << " " << yc << " " << radius << endl;
2357 
2358  for(int ihit = 0; ihit < nhits; ihit++) {
2359 
2360  Int_t hitId = (Int_t) points[ihit][0];
2361  Int_t detId = (Int_t) points[ihit][1];
2362 
2363  if(hitId == -1) continue;
2364  Int_t fitflag = (Int_t) points[ihit][10];
2365  if(fitflag != 0) continue;
2366  if(detId != FairRootManager::Instance()->GetBranchId(fSttBranchName)) continue;
2367  PndSttHit *stthit = (PndSttHit*) fSttHitArray->At(hitId);
2368  if(!stthit) continue;
2369 
2370  TVector3 xyz(0, 0, 0);
2371  TVector3 dxyz(0, 0, 0);
2372  Bool_t intfin = IntersectionFinder(xc, yc, radius, stthit, xyz, dxyz);
2373  if(intfin == false) {
2374  // cout << "intfin false2" << endl;
2375  continue;
2376  }
2377  points[ihit][2] = xyz.X();
2378  points[ihit][3] = xyz.Y();
2379  points[ihit][4] = xyz.Z();
2380  points[ihit][5] = dxyz.X();
2381  points[ihit][6] = dxyz.Y();
2382  points[ihit][7] = dxyz.Z();
2383  points[ihit][8] = stthit->GetIsochrone();
2384  points[ihit][9] = stthit->GetIsochroneError();
2385  lasthitid = ihit;
2386  if(firsthitid == -1) firsthitid = ihit;
2387 
2388 
2389  }
2390  }
2391 
2392  // z ---------------
2393  // ZFind(nhits, points, xc, yc, radius); // CHECK // to use zfinder
2394 
2395 
2396  Bool_t zfitting = ZFit(points, charge, xc, yc, radius, fitm, fitp);
2397  if(zfitting == false) {
2398  // cout << "zfit false " << endl;
2399  return false;
2400  }
2401 
2402 
2403  // cout << "ZFIT " << fitm << " " << fitp << endl;
2404  if(firsthitid == -1 || lasthitid == -1) return kFALSE; // CHECK
2405 
2406  // x0 y0
2407  Double_t d = TMath::Sqrt(xc * xc + yc * yc) - radius;
2408  Double_t phi = TMath::ATan2(yc, xc);
2409 
2410  Double_t x0 = d * TMath::Cos(phi);
2411  Double_t y0 = d * TMath::Sin(phi);
2412 
2413  // z
2414  Double_t Phi0 = TMath::ATan2((y0 - yc),(x0 - xc));
2415 
2416  TVector2 v(x0 - xc, y0 - yc);
2417  Double_t alpha1 = TMath::ATan2(points[firsthitid][3] - y0 + radius * TMath::Sin(Phi0), points[firsthitid][2] - x0 + radius * TMath::Cos(Phi0));
2418  TVector2 p1(points[firsthitid][2] - xc, points[firsthitid][3] - yc);
2419  Double_t Fi1 = CalculatePhi(v, p1, alpha1, Phi0, charge);
2420 
2421  Double_t alpha2 = TMath::ATan2(points[lasthitid][3] - y0 + radius * TMath::Sin(Phi0), points[lasthitid][2] - x0 + radius * TMath::Cos(Phi0));
2422  TVector2 p2(points[lasthitid][2] - xc, points[lasthitid][3] - yc);
2423  Double_t Fi2 = CalculatePhi(v, p2, alpha2, Phi0, charge);
2424  Fi2 = CompareToPreviousPhi(Fi2, Fi1, charge);
2425 
2426  Double_t scos = - charge * radius * Fi2; // scos = -q * R * phi CHECK :-)GOOD!
2427 
2428  double z = (fitm * scos + fitp);
2429 
2430  double versor[2];
2431  versor[0] = xc - points[lasthitid][2];
2432  versor[1] = yc - points[lasthitid][3];
2433 
2434  Double_t Distance = TMath::Sqrt(versor[0] * versor[0] + versor[1] * versor[1]);
2435  versor[0] /= Distance;
2436  versor[1] /= Distance;
2437 
2438  double px = - charge * radius * 0.006 * versor[1];
2439  double py = charge * radius * 0.006 * versor[0];
2440  double pz = TMath::Sqrt(px * px + py * py) * fitm; // CHECK :-)GOOD!
2441 
2442  lastpos.SetXYZ(points[lasthitid][2], points[lasthitid][3], z);
2443  lastmom.SetXYZ(px, py, pz);
2444 
2445 // cout << "PREFIT " << endl;
2446 // lastpos.Print();
2447 // lastmom.Print();
2448 
2449  return kTRUE;
2450 }
Double_t CompareToPreviousPhi(Double_t Fi, Double_t Fi_pre, int charge)
Double_t x0
Definition: checkhelixhit.C:70
TObjArray * d
Bool_t Fit(TMatrixT< double > points, Double_t &outxc, Double_t &outyc, Double_t &outradius)
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
static T Sin(const T &x)
Definition: PndCAMath.h:42
Double_t CalculatePhi(TVector2 v, TVector2 p, double alpha, double Phi0, int charge)
Bool_t IntersectionFinder(Double_t xc, Double_t yc, Double_t radius, PndSttHit *stthit, TVector3 &xyz, TVector3 &dxyz)
Bool_t ZFit(TMatrixT< double > points, Int_t charge, Double_t xc, Double_t yc, Double_t radius, Double_t &fitm, Double_t &fip)
static T Cos(const T &x)
Definition: PndCAMath.h:43
__m128 v
Definition: P4_F32vec4.h:4
TMatrixD GetCov() const
Definition: PndSdsHit.h:98
TClonesArray * fMvdPixelHitArray
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
TClonesArray * fMvdStripHitArray
Double_t y0
Definition: checkhelixhit.C:71
static T ATan2(const T &y, const T &x)
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
Int_t GetTubeID() const
Definition: PndSttHit.h:75
Double_t GetIsochroneError() const
Definition: PndSttHit.h:63
Double_t z
TPad * p2
Definition: hist-t7.C:117
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
TPad * p1
Definition: hist-t7.C:116
Int_t GetHitId() const
Bool_t GetInitialParams(PndTrack *sttmvd, Double_t &xc, Double_t &yc, Double_t &radius, Double_t &fitm, Double_t &fitp)
Int_t GetDetId() const
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
double pz[39]
Definition: pipisigmas.h:14
Bool_t PndSttMvdGemTracking::PropagateToGemPlane ( FairTrackParP *  tmppar,
FairTrackParP *  gempar,
Int_t  ipos 
)

Definition at line 992 of file PndSttMvdGemTracking.cxx.

References Bool_t, fDj, fDk, fPdgCode, fPro, fSensPositions, and fVerbose.

Referenced by Exec(), FillTrueDistances(), and Retrack().

992  {
993 
994  TVector3 *sensorpos = (TVector3*) fSensPositions->At(ipos);
995  TVector3 dj = tmppar->GetJVer();
996  TVector3 dk = tmppar->GetKVer();
997  fPro->PropagateFromPlane(dj, dk);
998 // fPro->PropagateFromPlane(fDj, fDk);
999  fPro->PropagateToPlane(*sensorpos, fDj, fDk);
1000 
1001  if(fVerbose > 0) {
1002  cout << "propagation from " << endl;
1003  tmppar->GetPosition().Print();
1004  tmppar->GetMomentum().Print();
1005  }
1006 
1007  // last z position
1008  // get 1st sensor in 1st station
1009  //PndGemStation *station = fGemParameters->GetStation(0); //[R.K. 01/2017] unused variable?
1010  //PndGemSensor *sensor = station->GetSensor(0); //[R.K. 01/2017] unused variable?
1011  if(tmppar->GetPosition().Z() > sensorpos->Z()) {
1012  if(fVerbose > 0) cout << "-> Z OUT OF BOUNDS: backpropagation" << endl;
1013  fPro->setBackProp();
1014  }
1015 
1016 
1017  Bool_t prop = fPro->Propagate(tmppar, gempar, fPdgCode);
1018  if(fVerbose > 0) {
1019  cout << "propagation to " << prop << endl;
1020  gempar->GetPosition().Print();
1021  gempar->GetMomentum().Print();
1022  }
1023  return prop;
1024 }
int fVerbose
Definition: poormantracks.C:24
Bool_t PndSttMvdGemTracking::PropagateToGemPlaneAsHelix ( PndTrack sttmvd,
FairTrackParP *  gempar,
Int_t  ipos 
)

Definition at line 1027 of file PndSttMvdGemTracking.cxx.

References alpha, CAMath::ATan2(), CalculatePhi(), CAMath::Cos(), d, Double_t, fabs(), fDj, fDk, fGemParameters, fOrdering, fOrderingIterator, fSensPositions, fVerbose, GetInitialParams(), PndGemStation::GetNSensors(), PndGemDigiPar::GetNStations(), PndGemSensor::GetOuterRadius(), PndTrack::GetParamFirst(), PndTrack::GetParamLast(), PndGemStation::GetSensor(), PndGemDigiPar::GetStation(), p, phi, pz, sensor, CAMath::Sin(), CAMath::Sqrt(), v, x, x0, y, y0, and z.

Referenced by Exec(), and FillTrueDistances().

1027  {
1028 
1029  TVector3 *sensorpos = (TVector3*) fSensPositions->At(ipos);
1030  Double_t z = sensorpos->Z();
1031 
1032  Double_t xc, yc, radius, fitm, fitp;
1033  GetInitialParams(sttmvd, xc, yc, radius, fitm, fitp);
1034  Int_t charge = sttmvd->GetParamFirst().GetQ();
1035 
1036  // z = fitp + scos * fitm
1037  if(fitm == 0) return kFALSE;
1038  Double_t scos = (z - fitp) / fitm; // CHECK :-)GOOD!
1039  double Fi = - charge * scos / radius; // CHECK :-)GOOD!
1040  // cout << "z1 " << z << " " << scos << " " << fitm << " " << fitp << endl;
1041 
1042  // x0 y0
1043  Double_t d = TMath::Sqrt(xc * xc + yc * yc) - radius;
1044  Double_t phi = TMath::ATan2(yc, xc);
1045 
1046  Double_t x0 = d * TMath::Cos(phi);
1047  Double_t y0 = d * TMath::Sin(phi);
1048  Double_t Phi0 = TMath::ATan2((y0 - yc),(x0 - xc));
1049 
1050  Double_t x = x0 + radius * (TMath::Cos(Phi0 + Fi) - TMath::Cos(Phi0));
1051  Double_t y = y0 + radius * (TMath::Sin(Phi0 + Fi) - TMath::Sin(Phi0));
1052 
1053  // ---
1054  // count sensors and stations and fill the sensor positions
1055  Int_t nstations = fGemParameters->GetNStations();
1056  Int_t nsensors = 0;
1057  double sens_rad = -1;
1058  Int_t posindex = 0;
1059  for(int istat = 0; istat < nstations; istat++) {
1060  PndGemStation *station = fGemParameters->GetStation(istat);
1061  for(int isens = 0; isens < station->GetNSensors(); isens++) {
1062  PndGemSensor *sensor = station->GetSensor(isens);
1063  nsensors++;
1064  posindex = (istat + 1) * 10 + (isens + 1);
1065  fOrderingIterator = find(fOrdering.begin(), fOrdering.end(), posindex);
1066  int jpos = fOrderingIterator - fOrdering.begin();
1067  if(ipos == jpos) {
1068  sens_rad = sensor->GetOuterRadius();
1069  break;
1070  }
1071  }
1072  }
1073  if(fabs(x) > sens_rad || fabs(y) > sens_rad) {
1074  if(fVerbose > 0) cout << "OUT OF SENSOR " << x << " " << y << " " << sens_rad << endl;
1075  return kFALSE;
1076  }
1077 
1078  // ---
1079 
1080  TVector2 v(x0 - xc, y0 - yc);
1081  Double_t alpha = TMath::ATan2(y - y0 + radius * TMath::Sin(Phi0), x - x0 + radius * TMath::Cos(Phi0));
1082  TVector2 p(x - xc, y - yc);
1083  Fi = CalculatePhi(v, p, alpha, Phi0, charge);
1084  scos = - charge * Fi * radius; // CHECK :-)GOOD!
1085  z = fitm * scos + fitp; // CHECK
1086 
1087  // cout << "z2 " << z << " " << scos << " " << fitm << " " << fitp << endl;
1088 
1089  double versor[2];
1090  versor[0] = xc - x;
1091  versor[1] = yc - y;
1092 
1093  Double_t Distance = TMath::Sqrt(versor[0] * versor[0] + versor[1] * versor[1]);
1094  versor[0] /= Distance;
1095  versor[1] /= Distance;
1096 
1097  double px = - charge * radius * 0.006 * versor[1];
1098  double py = charge * radius * 0.006 * versor[0];
1099  double pz = TMath::Sqrt(px * px + py * py) * fitm; // CHECK ?
1100 
1101  // error approx
1102  Double_t covMat[15];
1103  sttmvd->GetParamLast().GetCov(covMat); // CHECK it is not propagated
1104 
1105  gempar->SetTrackPar(x, y, z,
1106  px, py, pz,
1107  charge,
1108  covMat,
1109  *sensorpos, fDj.Cross(fDk), fDj, fDk);
1110 
1111 
1112 // cout << "x0, y0 " << x0 << " " << y0 << endl;
1113 // cout << "x y z px py pz " << endl;
1114 // cout << x << " " << y << " " << z << endl;
1115 // cout << px << " " << py << " " << pz << endl;
1116 
1117 
1118  return kTRUE;
1119 }
std::vector< int >::iterator fOrderingIterator
Double_t x0
Definition: checkhelixhit.C:70
Double_t p
Definition: anasim.C:58
int fVerbose
Definition: poormantracks.C:24
PndGemDigiPar * fGemParameters
TObjArray * d
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
static T Sin(const T &x)
Definition: PndCAMath.h:42
Double_t CalculatePhi(TVector2 v, TVector2 p, double alpha, double Phi0, int charge)
TGeoVolume * sensor
static T Cos(const T &x)
Definition: PndCAMath.h:43
__m128 v
Definition: P4_F32vec4.h:4
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
std::vector< int > fOrdering
Int_t GetNSensors() const
Definition: PndGemStation.h:60
Double_t
Double_t y0
Definition: checkhelixhit.C:71
static T ATan2(const T &y, const T &x)
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
Double_t z
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
Double_t x
Double_t GetOuterRadius() const
Definition: PndGemSensor.h:103
Double_t y
double alpha
Definition: f_Init.h:9
Bool_t GetInitialParams(PndTrack *sttmvd, Double_t &xc, Double_t &yc, Double_t &radius, Double_t &fitm, Double_t &fitp)
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
double pz[39]
Definition: pipisigmas.h:14
PndGemStation * GetStation(Int_t iStation)
void PndSttMvdGemTracking::Reset ( Int_t  nhits,
Int_t  ntracks 
)

Definition at line 371 of file PndSttMvdGemTracking.cxx.

References display, distancemap, fCombiMap, fDisplayOn, fGemParameters, fNPositions, fVerbose, PndGemStation::GetNSensors(), PndGemDigiPar::GetNStations(), PndGemSensor::GetOuterRadius(), PndGemStation::GetSensor(), PndGemDigiPar::GetStation(), h, hitcounter, hitmap, nhits, notassignedhits, sensor, towhichplane, trackindexes, trackvector, TString, and usabletracks.

Referenced by Exec().

371  {
372 
373  if(fVerbose > 0) cout << "******* PndSttMvdGemTracking::Reset *******" << endl;
374  if(fDisplayOn) {
375  char goOnChar;
376  cout << "press any key" << endl;
377  cin >> goOnChar;
378  cout << "GOING ON" << endl;
379 
380  display = new TCanvas("display", "display", 0, 0, 800, 500);
381  display->Divide(3, 2); // fNPositions/2, 2);
382  // GEM
383  int ipos = 0;
384  Int_t nstations = fGemParameters->GetNStations();
385  for(int istat = 0; istat < nstations; istat++) {
386  PndGemStation *station = fGemParameters->GetStation(istat);
387  for(int isens = 0; isens < station->GetNSensors(); isens++) {
388  PndGemSensor *sensor = station->GetSensor(isens);
389  TString iposname = "hipos"; TString ipostitle = "position ";
390  iposname += ipos; ipostitle += ipos;
391 
392  h[ipos] = new TH2F(iposname, ipostitle, 100, -sensor->GetOuterRadius(), sensor->GetOuterRadius(),
393  100, -sensor->GetOuterRadius(), sensor->GetOuterRadius());
394  display->cd(ipos + 1); h[ipos]->Draw();
395  display->Update();
396  display->Modified();
397  ipos++;
398  }
399  }
400  }
401 
402  hitmap.ResizeTo(fNPositions, nhits);
403  distancemap.ResizeTo(ntracks, nhits);
404 
405  notassignedhits.clear();
406  for(int ihit = 0; ihit < nhits; ihit++) {
407  for(int itrk = 0; itrk < ntracks; itrk++) distancemap[itrk][ihit] = -1;
408  notassignedhits.push_back(ihit);
409  }
410  if(fVerbose > 0) cout << "not assig " << nhits << " " << notassignedhits.size() << endl;
411 
412  for(int ipos = 0; ipos < fNPositions; ipos++) {
413  towhichplane[ipos] = false;
414 
415  // initialize all to 0
416  hitcounter(ipos, 0) = 0;
417 
418  for(int ihit = 0; ihit < nhits; ihit++) hitmap(ipos, ihit) = -1;
419  }
420 
421  if(fVerbose > 0) cout << "npositions/nhits " << fNPositions << " " << nhits << endl;
422  trackvector.clear();
423  usabletracks.clear(); // CHECK 4 PERFORMANCE delete this!
424  trackindexes.clear();
425 
426  // initialize the combimap to all combinatorials
427  fCombiMap.clear();
428  for(int ihit = 0; ihit < nhits; ihit++) fCombiMap[ihit] = 1;
429 }
int fVerbose
Definition: poormantracks.C:24
TMatrixT< float > hitmap
PndGemDigiPar * fGemParameters
std::vector< int > trackindexes
TGeoVolume * sensor
std::vector< std::pair< int, int > > trackvector
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
Int_t GetNSensors() const
Definition: PndGemStation.h:60
TMatrixT< double > distancemap
TMatrixT< float > hitcounter
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
Double_t GetOuterRadius() const
Definition: PndGemSensor.h:103
std::vector< int > notassignedhits
std::map< int, bool > towhichplane
PndGemStation * GetStation(Int_t iStation)
std::map< int, int > fCombiMap
std::vector< int > usabletracks
void PndSttMvdGemTracking::Retrack ( )

Definition at line 1472 of file PndSttMvdGemTracking.cxx.

References AddHitToTrack(), AssignHits(), CAMath::ATan2(), Bool_t, CHECKCOMBI, CAMath::Cos(), countprop, CountTracks(), DeleteHitFromTrack(), distancemap, Double_t, dPhi, fCombiMap, fCompleteTrackArray, fGemHitArray, fNPositions, fProTracks, fTurn, fVerbose, PndGemHit::GetDp(), PndGemHit::GetDr(), GetHitsAssociatedToTrack(), GetHitsAssociatedToTrackOnPlane(), PndTrack::GetParamLast(), PndTrack::GetTrackCandPtr(), GetTrackIndex(), IDEAL, kalman, Kalman(), par, phi, PropagateToGemPlane(), SetStartParameters(), CAMath::Sin(), CAMath::Sqrt(), and towhichplane.

Referenced by Exec().

1472  {
1473 
1474  if(fVerbose > 0) cout << "RETRACKING" << endl;
1475 
1476  for(Int_t k = 0; k < CountTracks(); k++) {
1477  int itrk = GetTrackIndex(k);
1478  if(fProTracks[itrk] == false) continue;
1479 
1480  std::vector<int> alreadyassociatedhits = GetHitsAssociatedToTrack(itrk);
1481  if(fVerbose > 0) {
1482  cout << "TRK " << itrk << " has hits ";
1483  for(size_t ihit = 0; ihit < alreadyassociatedhits.size(); ihit++) cout << " " << alreadyassociatedhits[ihit];
1484  cout << endl;
1485  }
1486 
1487  std::vector<int> hitvector = GetHitsAssociatedToTrack(itrk);
1488  //Int_t nhitinthistrack = hitvector.size(); //[R.K. 01/2017] unused variable?
1489 
1490  PndTrack* completeTrack = (PndTrack*) fCompleteTrackArray->At(itrk);
1491  if(!completeTrack) continue;
1492 
1493  // parameters on last stt
1494  FairTrackParP par = completeTrack->GetParamLast();
1495 
1496  FairTrackParP *gempar = new FairTrackParP();
1497 
1498  FairTrackParP tmppar = SetStartParameters(completeTrack, completeTrack->GetTrackCandPtr());
1499 
1500  for(int ipos = 0; ipos < fNPositions; ipos++) {
1501  // if no hit then continue
1502  if(towhichplane[ipos] == false) continue;
1503  std::vector<int> alreadyassociatedhitsonplane = GetHitsAssociatedToTrackOnPlane(itrk, ipos);
1504  if(alreadyassociatedhitsonplane.size() > 1) cout << "ERROR!! more than 1 hit on a plane assigned! " << alreadyassociatedhitsonplane.size() << endl; // CHECK delete this
1505 
1506  // ...else propagate here
1507  Bool_t prop = PropagateToGemPlane(&tmppar, gempar, ipos);
1508 
1509  if (prop == kFALSE) {
1510  countprop[fTurn-1][ipos]++;
1511  continue; // CHECK (continue or break?)
1512  }
1513 
1514 
1515  // cout << "propagation successful" << endl;
1516 
1517  // if the propagation was successful ...
1518  // assign the hits to track itrk, extrapolated to gempar on ipos
1519  std::vector<int> assignedhits = AssignHits(itrk, gempar, ipos);
1520 
1521  // cout << "on plane " << ipos << " there are " << assignedhits.size() << "(there were " << alreadyassociatedhitsonplane.size() << ")" << endl;
1522  // if there is no associated hit, continue
1523  if(assignedhits.size() == 0) {
1524  // if there was no hit associated before,
1525  // continue, since no kalman is possible
1526  if(alreadyassociatedhitsonplane.size() == 0) continue;
1527  else { // if there was a hit, delete it!
1528  int oldhit = alreadyassociatedhitsonplane[0];
1529  if(fVerbose > 0) cout << "delete old hit " << oldhit << endl;
1530  DeleteHitFromTrack(oldhit, itrk);
1531  continue;
1532  }
1533  }
1534 
1535  // if there is more than one hit, take the
1536  // minimum distance
1537  if(assignedhits.size() > 1) {
1538 
1539  std::vector<int>::iterator iter;
1540  Int_t tmphit = -1;
1541  Double_t tmpdistance = 1000;
1542 
1543  // loop over the hits
1544  for(iter = assignedhits.begin(); iter != assignedhits.end(); iter++) {
1545  std::vector<int>::iterator iter2;
1546  Int_t ihit = *iter;
1547  if(CHECKCOMBI > 0 && fCombiMap[ihit] != 0) { assignedhits.erase(iter); iter--; continue; }
1548  if(distancemap[itrk][ihit] < tmpdistance) {
1549  iter2 = std::find(assignedhits.begin(), assignedhits.end(), tmphit);
1550  int where = iter2 - assignedhits.begin();
1551  if(where != (int)assignedhits.size()) {
1552  if(fVerbose > 0) cout << "deleting old " << *iter2 << endl;
1553  assignedhits.erase(iter2); // remove from assigned list
1554  iter--;
1555  }
1556  tmpdistance = distancemap[itrk][ihit];
1557  tmphit = ihit;
1558  }
1559  else {
1560  if(fVerbose > 0) cout << "deleting new " << *iter << endl;
1561  assignedhits.erase(iter); // remove from assigned list
1562  iter--;
1563  }
1564  }
1565  }
1566  if(assignedhits.size() > 1) cout << "ERROR!! still more than 1 hit assigned!" << endl; // CHECK delete this
1567 
1568  // now we have only one hit assigned:
1569  // is it the one already assigned?
1570  std::vector<int>::iterator iter;
1571  iter = std::find(alreadyassociatedhits.begin(), alreadyassociatedhits.end(), assignedhits[0]);
1572  int where = iter - alreadyassociatedhits.begin();
1573 
1574  // if not...
1575  if(where == (int)alreadyassociatedhits.size()) {
1576  // if there was a previosly assigned hit
1577  // here, then make the substitution
1578  if(alreadyassociatedhitsonplane.size() > 0) {
1579  // there is only one hit on this plane
1580  int oldhit = alreadyassociatedhitsonplane[0];
1581  if(fVerbose > 0) cout << "deleting old hit2 " << oldhit << endl;
1582  DeleteHitFromTrack(oldhit, itrk);
1583  }
1584  if(fVerbose > 0) cout << "adding " << assignedhits[0] << " to track " << itrk << endl;
1585  AddHitToTrack(assignedhits[0], itrk);
1586  }
1587 
1588  if(GetHitsAssociatedToTrackOnPlane(itrk, ipos).size() == 0) continue;
1589  // retrieve hit
1590  Int_t finalhit = GetHitsAssociatedToTrackOnPlane(itrk, ipos)[0];
1591 
1592  PndGemHit *gemhit = (PndGemHit*) fGemHitArray->At(finalhit);
1593  if(!gemhit) continue;
1594 
1595  // kalman filter ----------------------
1596 
1597  // extrapolated point
1598  TMatrixT<double> extrap(5, 1);
1599  extrap[0][0] = gempar->GetQp();
1600  extrap[1][0] = gempar->GetTV();
1601  extrap[2][0] = gempar->GetTW();
1602  extrap[3][0] = gempar->GetV();
1603  extrap[4][0] = gempar->GetW();
1604 
1605  TMatrixT<double> extrap_cov(5, 5);
1606  Double_t covMat[15], covMat55[5][5];
1607  gempar->GetCov(covMat); // CHECK
1608  FairGeaneUtil util;
1609  util.FromVec15ToMat25(covMat, covMat55);
1610  for(int icov = 0; icov < 5; icov++) for(int jcov = 0; jcov < 5; jcov++) extrap_cov(icov, jcov) = covMat55[icov][jcov];
1611 
1612  // measurement matrix
1613  TMatrixT<double> H(2, 5);
1614  H[0][0] = 0.;
1615  H[0][1] = 0.;
1616  H[0][2] = 0.;
1617  H[0][3] = 1.;
1618  H[0][4] = 0.;
1619 
1620  H[1][0] = 0.;
1621  H[1][1] = 0.;
1622  H[1][2] = 0.;
1623  H[1][3] = 0.;
1624  H[1][4] = 1.;
1625 
1626  TMatrixT<double> measurement(2, 1);
1627  measurement[0][0] = gemhit->GetX();
1628  measurement[1][0] = gemhit->GetY();
1629 
1630  Double_t phi = TMath::ATan2(gemhit->GetY(), gemhit->GetX());
1631  Double_t sinp = TMath::Sin(phi);
1632  Double_t cosp = TMath::Cos(phi);
1633  //Double_t raddist = gemhit->GetPosition().Perp(); //[R.K. 01/2017] unused variable?
1634  Double_t dR, dP;
1635  if(IDEAL == true) {
1636  dR = 0.01;
1637  dP = 0.08;
1638  }
1639  else {
1640  dR = gemhit->GetDr();
1641  dP = gemhit->GetDp();
1642  }
1643  Double_t dPhi = dP; // / raddist;
1644  Double_t errx = TMath::Sqrt(cosp * cosp * dR * dR + sinp * sinp * dPhi * dPhi);
1645  Double_t erry = TMath::Sqrt(sinp * sinp * dR * dR + cosp * cosp * dPhi * dPhi);
1646 
1647  // CHECK covariances??
1648  TMatrixT<double> measurement_cov(2, 2);
1649  measurement_cov[0][0] = errx * errx;
1650  measurement_cov[1][1] = erry * erry;
1651 
1652  TMatrixT<double> kalman(5, 1);
1653  TMatrixT<double> kalman_cov(5, 5);
1654 
1655  // kalman
1656  Kalman(extrap, measurement, H, extrap_cov, measurement_cov, kalman, kalman_cov);
1657 
1658  Double_t kalmanCov15[15], kalmanCov55[5][5];
1659  for(int icov = 0; icov < 5; icov++) for(int jcov = 0; jcov < 5; jcov++) kalmanCov55[icov][jcov] = kalman_cov(icov, jcov);
1660  util.FromMat25ToVec15(kalmanCov55, kalmanCov15);
1661 
1662  tmppar.SetTrackPar(kalman[3][0], kalman[4][0],
1663  kalman[1][0], kalman[2][0], kalman[0][0],
1664  // covMat,
1665  kalmanCov15, // CHECK NOW
1666  gempar->GetOrigin(), gempar->GetIVer(), gempar->GetJVer(), gempar->GetKVer(),
1667  gempar->GetSPU()); // CHECK recalculate spu
1668 
1669 
1670  }
1671  }
1672 }
int fVerbose
Definition: poormantracks.C:24
#define IDEAL
std::vector< int > GetHitsAssociatedToTrackOnPlane(Int_t itrk, Int_t ipos)
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
static T Sin(const T &x)
Definition: PndCAMath.h:42
Double_t par[3]
KalmanTask * kalman
static T Cos(const T &x)
Definition: PndCAMath.h:43
TClonesArray * fCompleteTrackArray
std::vector< int > GetHitsAssociatedToTrack(Int_t itrk)
Bool_t PropagateToGemPlane(FairTrackParP *tmppar, FairTrackParP *gempar, Int_t ipos)
Double_t
TMatrixT< double > distancemap
static T ATan2(const T &y, const T &x)
Double_t GetDp() const
Definition: PndGemHit.h:76
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
void Kalman(TMatrixT< double > extrap, TMatrixT< double > measurement, TMatrixT< double > H, TMatrixT< double > extrap_cov, TMatrixT< double > measurement_cov, TMatrixT< double > &kalman, TMatrixT< double > &kalman_cov)
std::map< int, bool > fProTracks
void DeleteHitFromTrack(Int_t ihit, Int_t itrk)
#define CHECKCOMBI
void AddHitToTrack(Int_t ihit, Int_t itrk)
FairTrackParP SetStartParameters(PndTrack *sttmvd, PndTrackCand *sttmvdCand)
std::vector< int > AssignHits(Int_t itrk, FairTrackParP *gempar, Int_t ipos)
Double_t GetDr() const
Definition: PndGemHit.h:75
PndTrackCand * GetTrackCandPtr()
Definition: PndTrack.h:48
double dPhi
Definition: RiemannTest.C:17
std::map< int, bool > towhichplane
std::map< int, int > fCombiMap
Int_t PndSttMvdGemTracking::SelectPdgCode ( PndTrack sttmvd)

Definition at line 2057 of file PndSttMvdGemTracking.cxx.

References Bool_t, fPro, PndTrack::GetParamFirst(), and PndTrack::GetParamLast().

2057  {
2058 
2059  FairTrackParP first = sttmvd->GetParamFirst();
2060  FairTrackParP last = sttmvd->GetParamLast();
2061 
2062  // first plane
2063  // TVector3 o1 = firs.GetOrigin();
2064  TVector3 dj1 = first.GetJVer();
2065  TVector3 dk1 = first.GetKVer();
2066 
2067  // last plane
2068  TVector3 o2 = last.GetOrigin();
2069  TVector3 dj2 = last.GetJVer();
2070  TVector3 dk2 = last.GetKVer();
2071 
2072  fPro->PropagateFromPlane(dj1,dk1);
2073  fPro->PropagateToPlane(o2, dj2, dk2);
2074 
2075  FairTrackParP *propag = new FairTrackParP();
2076  Bool_t prop = kFALSE;
2077  int charge = first.GetQ();
2078  int pdg[5] = {-11 * charge, -13 * charge, 211 * charge, 321 * charge, 2212 * charge};
2079  int tmppdg = 0;
2080  double tmpdistance = 100000;
2081  for(int ipdg = 0; ipdg < 5; ipdg++) {
2082  prop = fPro->Propagate(&first, propag, pdg[ipdg]);
2083  if(prop) {
2084  TVector3 proppos = propag->GetPosition();
2085  TVector3 lastpos = last.GetPosition();
2086  double distance = (lastpos - proppos).Mag();
2087  cout << "this pdg " << pdg[ipdg] << " " << distance << endl;
2088 
2089  if(distance < tmpdistance) {
2090  tmpdistance = distance;
2091  tmppdg = pdg[ipdg];
2092  }
2093  }
2094  }
2095  cout << "CHOSEN PDG " << tmppdg << " " << tmpdistance << endl;
2096  return tmppdg;
2097 }
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
void PndSttMvdGemTracking::SetBranchNames ( TString  mvdpixel,
TString  mvdstrip,
TString  stt,
TString  gem 
)

Definition at line 3513 of file PndSttMvdGemTracking.cxx.

References fGemBranchName, fMvdPixelBranchName, fMvdStripBranchName, and fSttBranchName.

3513  {
3514  fMvdPixelBranchName = mvdpixel;
3515  fMvdStripBranchName = mvdstrip;
3516  fSttBranchName = stt;
3517  fGemBranchName = gem;
3518 }
void PndSttMvdGemTracking::SetCombinatorialDistance ( Double_t  combidistance)
inline

Definition at line 107 of file PndSttMvdGemTracking.h.

References fCombiDistance.

107 { fCombiDistance = combidistance; }
void PndSttMvdGemTracking::SetDefaultPdg ( int  pdg)
inline

Definition at line 127 of file PndSttMvdGemTracking.h.

References fDefaultPdgCode.

127 { fDefaultPdgCode = pdg; }
void PndSttMvdGemTracking::SetEvaluateFlag ( Bool_t  flag)
inline

Definition at line 117 of file PndSttMvdGemTracking.h.

References fEvaluate.

117 { fEvaluate = flag; }
void PndSttMvdGemTracking::SetMaxDistance ( Double_t  maxdistance)
inline

Definition at line 50 of file PndSttMvdGemTracking.h.

References fMaxDistance.

50 {fMaxDistance = maxdistance;}
void PndSttMvdGemTracking::SetParContainers ( )
virtual

Definition at line 268 of file PndSttMvdGemTracking.cxx.

References fGemParameters, fSttParameters, and rtdb.

268  {
269 
270  FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
271 
272  // get STT parameters
273  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
274 
275  // get GEM parameters
276  fGemParameters = (PndGemDigiPar*) rtdb->getContainer("PndGemDetectors");
277 
278 
279 }
PndGemDigiPar * fGemParameters
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
void PndSttMvdGemTracking::SetPdgFromMC ( )
inline

Definition at line 122 of file PndSttMvdGemTracking.h.

References fPdgFromMC.

122 { fPdgFromMC = kTRUE; }
void PndSttMvdGemTracking::SetPdgFromMC ( TString  trackid)
inline

Definition at line 123 of file PndSttMvdGemTracking.h.

References fPdgFromMC, and fStartTrackIDBranchName.

123  {
124  fStartTrackIDBranchName = trackid;
125  fPdgFromMC = kTRUE;
126  }
void PndSttMvdGemTracking::SetPersistenc ( Bool_t  persistence)
inline

set persistence flag

Definition at line 41 of file PndSttMvdGemTracking.h.

References PndPersistencyTask::SetPersistency().

41 { SetPersistency(persistence); }
void SetPersistency(Bool_t val=kTRUE)
void PndPersistencyTask::SetPersistency ( Bool_t  val = kTRUE)
inlineinherited

Definition at line 31 of file PndPersistencyTask.h.

References PndPersistencyTask::fPersistency, and val.

Referenced by barrelTrackFinder(), digi_complete(), digi_complete_newSTT(), digiOnly_complete(), PndBarrelTrackFinder::PndBarrelTrackFinder(), PndCATracking::PndCATracking(), PndDrcHitFinder::PndDrcHitFinder(), PndEmc2DLocMaxFinder::PndEmc2DLocMaxFinder(), PndEmcExpClusterSplitter::PndEmcExpClusterSplitter(), PndEmcFullDigiTask::PndEmcFullDigiTask(), PndEmcFWEndcapDigi::PndEmcFWEndcapDigi(), PndEmcFWEndcapTimebasedWaveforms::PndEmcFWEndcapTimebasedWaveforms(), PndEmcHitProducer::PndEmcHitProducer(), PndEmcHitsToWaveform::PndEmcHitsToWaveform(), PndEmcMakeBump::PndEmcMakeBump(), PndEmcMakeCluster::PndEmcMakeCluster(), PndEmcMakeClusterOnline::PndEmcMakeClusterOnline(), PndEmcMakeDigi::PndEmcMakeDigi(), PndEmcMakeRecoHit::PndEmcMakeRecoHit(), PndEmcMultiWaveformToCalibratedDigi::PndEmcMultiWaveformToCalibratedDigi(), PndEmcPhiBumpSplitter::PndEmcPhiBumpSplitter(), PndEmcTmpWaveformToDigi::PndEmcTmpWaveformToDigi(), PndEmcWaveformToCalibratedDigi::PndEmcWaveformToCalibratedDigi(), PndEmcWaveformToDigi::PndEmcWaveformToDigi(), PndFtofHitProducerIdeal::PndFtofHitProducerIdeal(), PndFtsCATracking::PndFtsCATracking(), PndFtsHitProducerIdeal::PndFtsHitProducerIdeal(), PndFtsHitProducerRealFast::PndFtsHitProducerRealFast(), PndFtsHitProducerRealFull::PndFtsHitProducerRealFull(), PndFtsHoughTrackerTask::PndFtsHoughTrackerTask(), PndGemDigitize::PndGemDigitize(), PndGemFindHits::PndGemFindHits(), PndIdealTrackFinder::PndIdealTrackFinder(), PndLmdPixelClusterTask::PndLmdPixelClusterTask(), PndLmdPixelHitProducerFast::PndLmdPixelHitProducerFast(), PndMdtClusterTask::PndMdtClusterTask(), PndMdtDigitization::PndMdtDigitization(), PndMdtHitProducerIdeal::PndMdtHitProducerIdeal(), PndMdtPointsToWaveform::PndMdtPointsToWaveform(), PndMdtTrkProducer::PndMdtTrkProducer(), PndMissingPzCleanerTask::PndMissingPzCleanerTask(), PndMvdGemTrackFinderOnHits::PndMvdGemTrackFinderOnHits(), PndMvdHitProducerIdeal::PndMvdHitProducerIdeal(), PndMvdPixelClusterTask::PndMvdPixelClusterTask(), PndMvdTimeWalkCorrTask::PndMvdTimeWalkCorrTask(), PndMvdToPix4ClusterTask::PndMvdToPix4ClusterTask(), PndRecoKalmanTask::PndRecoKalmanTask(), PndRecoKalmanTask2::PndRecoKalmanTask2(), PndRecoMultiKalmanTask::PndRecoMultiKalmanTask(), PndRecoMultiKalmanTask2::PndRecoMultiKalmanTask2(), PndRichHitFinder::PndRichHitFinder(), PndRichHitProducer::PndRichHitProducer(), PndSciTDigiTask::PndSciTDigiTask(), PndSciTHitProducerIdeal::PndSciTHitProducerIdeal(), PndSdsHitProducerIdeal::PndSdsHitProducerIdeal(), PndSdsHybridHitProducer::PndSdsHybridHitProducer(), PndSdsIdealClusterTask::PndSdsIdealClusterTask(), PndSdsIdealRecoTask::PndSdsIdealRecoTask(), PndSdsNoiseProducer::PndSdsNoiseProducer(), PndSdsPixelClusterTask::PndSdsPixelClusterTask(), PndSdsStripClusterTask::PndSdsStripClusterTask(), PndSdsStripHitProducer::PndSdsStripHitProducer(), PndSdsTimeWalkCorrTask::PndSdsTimeWalkCorrTask(), PndSttFindTracks::PndSttFindTracks(), PndSttHelixHitProducer::PndSttHelixHitProducer(), PndSttHitProducerIdeal::PndSttHitProducerIdeal(), PndSttHitProducerRealFast::PndSttHitProducerRealFast(), PndSttHitProducerRealFull::PndSttHitProducerRealFull(), PndSttMatchTracks::PndSttMatchTracks(), PndSttMvdGemTracking(), PndTrackSmearTask::PndTrackSmearTask(), PndTrkTracking2::PndTrkTracking2(), reco(), reco_complete(), reco_complete_gf2(), reco_complete_newSTT(), reco_complete_sec(), recoideal_complete(), PndMvdClusterTask::SetPersistance(), PndMvdDigiTask::SetPersistance(), PndLmdPixelHitProducerFast::SetPersistance(), PndSdsHitProducerIdeal::SetPersistance(), SetPersistenc(), PndMdtClusterTask::SetPersistence(), PndSttHelixHitProducer::SetPersistence(), PndMissingPzCleanerTask::SetPersistence(), PndFtsHitProducerRealFast::SetPersistence(), PndFtsHitProducerRealFull::SetPersistence(), PndSttHitProducerIdeal::SetPersistence(), PndSttHitProducerRealFull::SetPersistence(), PndSttHitProducerRealFast::SetPersistence(), PndFtsHitProducerIdeal::SetPersistence(), PndTrackSmearTask::SetPersistence(), PndSciTHitProducerIdeal::SetPersistence(), PndIdealTrackFinder::SetPersistence(), PndSttMatchTracks::SetPersistence(), PndSttFindTracks::SetPersistence(), PndFtsHoughTrackerTask::SetPersistence(), PndTrkTracking2::SetPersistence(), PndEmcMakeRecoHit::SetStorageOfData(), PndEmcMakeClusterOnline::SetStorageOfData(), PndEmcFWEndcapDigi::SetStorageOfData(), PndEmcFWEndcapTimebasedWaveforms::SetStorageOfData(), PndEmcMakeDigi::SetStorageOfData(), PndMdtPointsToWaveform::SetStorageOfData(), PndEmc2DLocMaxFinder::SetStorageOfData(), PndEmcMakeCluster::SetStorageOfData(), PndEmcHitsToWaveform::SetStorageOfData(), PndEmcMakeBump::SetStorageOfData(), PndEmcTmpWaveformToDigi::SetStorageOfData(), PndEmcWaveformToDigi::SetStorageOfData(), PndEmcWaveformToCalibratedDigi::SetStorageOfData(), PndEmcMultiWaveformToCalibratedDigi::SetStorageOfData(), PndEmcExpClusterSplitter::SetStorageOfData(), PndEmcPhiBumpSplitter::SetStorageOfData(), standard_tracking(), and PndEmcFullDigiTask::StoreDigi().

31 { fPersistency = val; }
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
FairTrackParP PndSttMvdGemTracking::SetStartParameters ( PndTrack sttmvd,
PndTrackCand sttmvdCand 
)

Definition at line 1876 of file PndSttMvdGemTracking.cxx.

References fMCTrackArray, fMvdPixelBranchName, fMvdPixelHitArray, fMvdPointArray, fMvdStripBranchName, fMvdStripHitArray, fSttBranchName, fSttHitArray, fSttPointArray, fUseMC, fVerbose, PndTrackCandHit::GetDetId(), PndTrackCandHit::GetHitId(), PndTrackCand::getMcTrackId(), PndTrackCand::GetNHits(), PndTrack::GetParamLast(), PndMCTrack::GetPdgCode(), PndTrackCand::GetSortedHit(), and Prefit().

Referenced by Exec(), FillTrueDistances(), and Retrack().

1876  {
1877 
1878  // in the following cuts, if it does not work, return this!
1879  FairTrackParP lastpar = sttmvd->GetParamLast();
1880 
1881  if(fVerbose > 0) {
1882  cout << "start from " << endl;
1883  lastpar.GetPosition().Print();
1884  cout << "lastpar error " << lastpar.GetDX() << " " << lastpar.GetDY() << " " << lastpar.GetDZ() << endl;
1885  // lastpar.GetMomentum().Print();
1886  }
1887 
1888  FairTrackParP startpar;
1889  if(fUseMC) {
1890  // --------------------------------------------------
1891  Int_t mcIndex = sttmvdCand->getMcTrackId();
1892  if(mcIndex == -1) return lastpar;
1893  PndMCTrack *mctrk = (PndMCTrack*) fMCTrackArray->At(mcIndex);
1894  // FairTrackParP tmppar(mctrk->GetStartVertex(), mctrk->GetMomentum(),
1895  // TVector3(0.1, 0.1, 0.1), TVector3(0.1, 0.1, 0.1),
1896  // TDatabasePDG::Instance()->GetParticle(mctrk->GetPdgCode())->Charge()/3.,
1897  // mctrk->GetStartVertex(), TVector3(1., 0., 0.), TVector3(0., 1., 0.));
1898 
1899  // start from MC last point on STT
1900  Int_t nsttmvdhits = sttmvdCand->GetNHits();
1901  PndTrackCandHit candhit = sttmvdCand->GetSortedHit(nsttmvdhits - 1);
1902  Int_t iHit = candhit.GetHitId();
1903  Int_t detId = candhit.GetDetId();
1904  FairHit *fhit = NULL;
1905  FairMCPoint *fpnt = NULL;
1906 
1907 // cout << "LAST HIT " << detId
1908 // << " " << FairRootManager::Instance()->GetBranchId(fMvdStripBranchName)
1909 // << " " << FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName)
1910 // << " " << FairRootManager::Instance()->GetBranchId(fSttBranchName)
1911 // << endl;
1912 
1913  if(detId != FairRootManager::Instance()->GetBranchId(fSttBranchName)) return lastpar;
1914 
1915  if(detId == FairRootManager::Instance()->GetBranchId(fMvdStripBranchName))
1916  {
1917  fhit = (FairHit*) fMvdStripHitArray->At(iHit);
1918  if(!fhit) return lastpar;
1919  fpnt = (FairMCPoint*) fMvdPointArray->At(fhit->GetRefIndex());
1920  if(!fpnt) return lastpar;
1921  }
1922  else if(detId == FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName)) {
1923  fhit = (FairHit*) fMvdPixelHitArray->At(iHit);
1924  if(!fhit) return lastpar;
1925  fpnt = (FairMCPoint*) fMvdPointArray->At(fhit->GetRefIndex());
1926  if(!fpnt) return lastpar;
1927  }
1928  else if(detId == FairRootManager::Instance()->GetBranchId(fSttBranchName)) {
1929  fhit = (FairHit*) fSttHitArray->At(iHit);
1930  if(!fhit) return lastpar;
1931  fpnt = (FairMCPoint*) fSttPointArray->At(fhit->GetRefIndex());
1932  if(!fpnt) return lastpar;
1933 
1934  }
1935 
1936  startpar = FairTrackParP(TVector3(fpnt->GetX(), fpnt->GetY(), fpnt->GetZ()),
1937  TVector3(fpnt->GetPx(), fpnt->GetPy(), fpnt->GetPz()),
1938  TVector3(0.1, 0.1, 0.1), TVector3(0.1, 0.1, 0.1),
1939  (int) TMath::Sign(1., TDatabasePDG::Instance()->GetParticle(mctrk->GetPdgCode())->Charge()/3.),
1940  TVector3(fpnt->GetX(), fpnt->GetY(), fpnt->GetZ()),
1941  TVector3(1., 0., 0.), TVector3(0., 1., 0.));
1942  }
1943  else {
1944  // startpar = lastpar;
1945  // cout << "from PR " << endl;
1946  // lastpar.GetPosition().Print();
1947  // lastpar.GetMomentum().Print();
1948 
1949  bool startpoint = false;
1950  //PndGemStation *station = fGemParameters->GetStation(0); //[R.K. 01/2017] unused variable?
1951  //PndGemSensor *sensor = station->GetSensor(0); //[R.K. 01/2017] unused variable?
1952  TVector3 position;
1953  TVector3 momentum;
1954  // if(fabs(lastpar.GetPosition().X()) > 42. || fabs(lastpar.GetPosition().Y()) > 42. ||
1955  // lastpar.GetMomentum().Mag() < 0.15 ||
1956  // lastpar.GetPosition().Z() > sensor->GetZ0())
1957  {
1958  startpoint = Prefit(sttmvd, sttmvdCand, position, momentum);
1959  }
1960 
1961  if(startpoint == true) {
1962  // find plane orthogonal to mom
1963  TVector3 iver = momentum; iver.SetMag(1.);
1964  TVector3 jver = momentum.Orthogonal(); jver.SetMag(1.);
1965  TVector3 kver = iver.Cross(jver); kver.SetMag(1.);
1966  // by hand
1967  startpar = FairTrackParP(position,
1968  momentum,
1969  TVector3(lastpar.GetDX(), lastpar.GetDY(), 3.),
1970  TVector3(lastpar.GetDPx(), lastpar.GetDPy(), lastpar.GetDPz()),
1971  lastpar.GetQ(),
1972  position, jver, kver);
1973  }
1974  else {
1975  // by hand
1976  startpar = FairTrackParP(lastpar.GetPosition(),
1977  lastpar.GetMomentum(),
1978  TVector3(lastpar.GetDX(), lastpar.GetDY(), 3.),
1979  TVector3(lastpar.GetDPx(), lastpar.GetDPy(), lastpar.GetDPz()),
1980  lastpar.GetQ(),
1981  lastpar.GetOrigin(), lastpar.GetJVer(), lastpar.GetKVer());
1982  }
1983 
1984  // kalman filter between the lastpar and the last hit
1985 
1986  // Int_t nsttmvdhits = sttmvdCand->GetNHits();
1987  // PndTrackCandHit candhit = sttmvdCand->GetSortedHit(nsttmvdhits - 1);
1988  // Int_t iHit = candhit.GetHitId();
1989  // Int_t detId = candhit.GetDetId();
1990 
1991  // // CHECK
1992  // if(detId != FairRootManager::Instance()->GetBranchId(fSttBranchName)) return lastpar;
1993 
1994  // PndSttHit *stthit = (PndSttHit*) fSttHitArray->At(iHit);
1995  // if(!stthit) return lastpar;
1996 
1997  // // extrapolated point
1998  // TMatrixT<double> extrap(5, 1);
1999  // extrap[0][0] = lastpar.GetQp();
2000  // extrap[1][0] = lastpar.GetTV();
2001  // extrap[2][0] = lastpar.GetTW();
2002  // extrap[3][0] = lastpar.GetV();
2003  // extrap[4][0] = lastpar.GetW();
2004 
2005  // TMatrixT<double> extrap_cov(5, 5);
2006  // Double_t covMat[15], covMat55[5][5];
2007  // lastpar.GetCov(covMat); // CHECK
2008  // FairGeaneUtil util;
2009  // util.FromVec15ToMat25(covMat, covMat55);
2010  // for(int icov = 0; icov < 5; icov++) for(int jcov = 0; jcov < 5; jcov++) extrap_cov(icov, jcov) = covMat55[icov][jcov];
2011 
2012  // // measurement matrix
2013  // TMatrixT<double> H(1, 5);
2014  // H[0][0] = 0.;
2015  // H[0][1] = 0.;
2016  // H[0][2] = 0.;
2017  // H[0][3] = 1.;
2018  // H[0][4] = 0.;
2019 
2020  // TMatrixT<double> measurement(1, 1);
2021  // measurement[0][0] = stthit->GetIsochrone();
2022 
2023  // TMatrixT<double> measurement_cov(1, 1);
2024  // measurement_cov[0][0] = stthit->GetIsochroneError() * stthit->GetIsochroneError();
2025 
2026  // TMatrixT<double> kalman(5, 1);
2027  // TMatrixT<double> kalman_cov(5, 5);
2028 
2029  // // kalman
2030  // Kalman(extrap, measurement, H, extrap_cov, measurement_cov, kalman, kalman_cov);
2031 
2032  // Double_t kalmanCov15[15], kalmanCov55[5][5];
2033  // for(int icov = 0; icov < 5; icov++) for(int jcov = 0; jcov < 5; jcov++) kalmanCov55[icov][jcov] = kalman_cov(icov, jcov);
2034  // util.FromMat25ToVec15(kalmanCov55, kalmanCov15);
2035 
2036 
2037  // startpar.SetTrackPar(kalman[3][0], kalman[4][0],
2038  // kalman[1][0], kalman[2][0], kalman[0][0],
2039  // kalmanCov15, // CHECK!!
2040  // lastpar.GetOrigin(), lastpar.GetIVer(), lastpar.GetJVer(), lastpar.GetKVer(),
2041  // lastpar.GetSPU()); // CHECK recalculate spu
2042  // cout << "from prefit " << endl;
2043  // startpar.GetPosition().Print();
2044  // startpar.GetMomentum().Print();
2045 
2046  }
2047 
2048  if(fVerbose) {
2049  cout << "from prefit " << endl;
2050  startpar.GetPosition().Print();
2051  startpar.GetMomentum().Print();
2052  }
2053 
2054  return startpar;
2055 }
int fVerbose
Definition: poormantracks.C:24
int getMcTrackId() const
Definition: PndTrackCand.h:60
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
TClonesArray * fMvdPixelHitArray
TClonesArray * fMvdStripHitArray
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
Bool_t Prefit(PndTrack *sttmvdTrack, PndTrackCand *sttmvdCand, TVector3 &lastpos, TVector3 &lastmom)
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
Int_t GetHitId() const
Int_t GetDetId() const
void PndSttMvdGemTracking::SetTimes ( Int_t  times)
inline

Definition at line 49 of file PndSttMvdGemTracking.h.

References fTimes.

49 {fTimes = times;}
void PndSttMvdGemTracking::SetTrackBranchNames ( TString  startingtrack,
TString  startingtrackcand 
)

Definition at line 3520 of file PndSttMvdGemTracking.cxx.

References fStartTrackBranchName, and fStartTrackCandBranchName.

3521 {
3522  fStartTrackBranchName = starttrack;
3523  fStartTrackCandBranchName = starttrackcand;
3524 
3525 }
void PndSttMvdGemTracking::SetupGEMPlanes ( )

Definition at line 281 of file PndSttMvdGemTracking.cxx.

References a, fDj, fDk, fGemParameters, fNPositions, fOrdering, fSensPositions, fVerbose, PndGemStation::GetNSensors(), PndGemDigiPar::GetNStations(), PndGemStation::GetSensor(), PndGemDigiPar::GetStation(), PndGemSensor::GetX0(), PndGemSensor::GetY0(), PndGemSensor::GetZ0(), hchosen, hchosen2, hdist, hdist2, hitcounter, hmcdist, hmcx, hmcy, hsigma, hsigma2, i, sensor, and TString.

Referenced by Init().

281  {
282 
283 
284  if(fVerbose > 0) cout << "********* SPndSttMvdGemTracking::SetupGEMPlanes ************" << endl;
285 
286  // geometry of GEM stations
287  Int_t nstations = fGemParameters->GetNStations();
288  Int_t nsensors = 0;
289  Int_t posindex = 0;
290 
291  if(fVerbose > 0) cout << "NSTAT " << nstations << endl;
292 
293  if(nstations == 0) cerr << "PndSttMvdGemTracking::SetupGEMPlanes: NO GEOMETRY FOR GEMs!" << endl;
294 
295  // position of the planes
296  fSensPositions = new TClonesArray("TVector3", 100);
297 
298  // order of the sensors:
299  // stat = 1, sens = 1 ==> posindex = 11 ==> ordering = 0
300  // stat = 1, sens = 2 ==> posindex = 12 ==> ordering = 1
301  // stat = 2, sens = 1 ==> posindex = 21 ==> ordering = 2
302  // stat = 2, sens = 2 ==> posindex = 22 ==> ordering = 3
303  // ...
304 
305  // count sensors and stations and fill the sensor positions
306  for(int istat = 0; istat < nstations; istat++) {
307  PndGemStation *station = fGemParameters->GetStation(istat);
308  for(int isens = 0; isens < station->GetNSensors(); isens++) {
309  PndGemSensor *sensor = station->GetSensor(isens);
310  nsensors++;
311  posindex = (istat + 1) * 10 + (isens + 1);
312  fOrdering.push_back(posindex);
313  if(fVerbose > 0) cout << "istat/isens/posindex/ordering " << istat + 1 << " " << isens + 1 << " " << posindex << " " << fOrdering.size() << " " << fOrdering[fOrdering.size() - 1] << endl;
314 
315  // fill the TClonesArray with the origin
316  // of the planes corresponding to each
317  // sensor of GEM
318  int size = fOrdering.size() - 1;
319  TClonesArray& clref = *fSensPositions;
320  new(clref[size]) TVector3(sensor->GetX0(), sensor->GetY0(), sensor->GetZ0());
321 
322  }
323  }
324 
325  // the sensor planes are parallel to xy // CHECK is this true?
326  fDj.SetXYZ(1., 0., 0.);
327  fDk.SetXYZ(0., 1., 0.);
328 
329  // CHECK delete this ---
330  if(fVerbose > 0) {
331  cout << "entries " << fSensPositions->GetEntriesFast() << endl;
332  for(int i = 0; i < fSensPositions->GetEntriesFast(); i++) {
333  TVector3 *a = (TVector3*) fSensPositions->At(i);
334  cout << i << endl;
335  if(!a) continue;
336  a->Print();
337  }
338  }
339  // ---
340 
341 
342  fNPositions = nsensors; // * nstations;
343  for(int ipos = 0; ipos < fNPositions; ipos++) {
344  TString hname = "hdist"; hname += ipos;
345  hdist[ipos] = new TH1F(hname, hname, 200, 0, 100);
346  hname += " turn 2";
347  hdist2[ipos] = new TH1F(hname, hname, 200, 0, 100);
348  hname = "hsigma"; hname += ipos;
349  hsigma[ipos] = new TH1F(hname, hname, 200, 0, 20);
350  hname += " turn 2";
351  hsigma2[ipos] = new TH1F(hname, hname, 200, 0, 20);
352  hname = "hchosen"; hname += ipos;
353  hchosen[ipos] = new TH1F(hname, hname, 200, 0, 100);
354  hname += " turn 2";
355  hchosen2[ipos] = new TH1F(hname, hname, 200, 0, 100);
356  hname = "hmcdist"; hname += ipos;
357  hmcdist[ipos] = new TH1F(hname, hname, 200, 0, 100);
358  hname = "hmcx"; hname += ipos;
359  hmcx[ipos] = new TH1F(hname, hname, 200, 0, 100);
360  hname = "hmcy"; hname += ipos;
361  hmcy[ipos] = new TH1F(hname, hname, 200, 0, 100);
362 
363  }
364 
365  // resize hitcounter
366  hitcounter.ResizeTo(fNPositions, 1);
367 
368 }
int fVerbose
Definition: poormantracks.C:24
PndGemDigiPar * fGemParameters
Int_t i
Definition: run_full.C:25
TGeoVolume * sensor
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
std::vector< int > fOrdering
Int_t a
Definition: anaLmdDigi.C:126
Int_t GetNSensors() const
Definition: PndGemStation.h:60
Double_t GetX0() const
Definition: PndGemSensor.h:98
TMatrixT< float > hitcounter
Double_t GetY0() const
Definition: PndGemSensor.h:99
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
Double_t GetZ0() const
Definition: PndGemSensor.h:100
PndGemStation * GetStation(Int_t iStation)
void PndSttMvdGemTracking::SwitchOnDisplay ( )
inline

Definition at line 51 of file PndSttMvdGemTracking.h.

References fDisplayOn.

51 {fDisplayOn = true;}
void PndSttMvdGemTracking::UpdateMCTrackId ( PndTrackCand completeCand)

Definition at line 3390 of file PndSttMvdGemTracking.cxx.

References counter, fGemBranchName, fGemHitArray, fGemPointArray, fMvdPixelBranchName, fMvdPixelHitArray, fMvdPointArray, fMvdStripBranchName, fMvdStripHitArray, fSttBranchName, fSttHitArray, fSttPointArray, PndTrackCandHit::GetDetId(), PndTrackCandHit::GetHitId(), PndTrackCand::getMcTrackId(), PndTrackCand::GetNHits(), PndTrackCand::GetSortedHit(), refIndex, and PndTrackCand::setMcTrackId().

3390  {
3391 
3392 
3393  cout << "UPDATE MC TRACK ID" << endl;
3394  cout << "ORIGINAL MC " << completeCand->getMcTrackId() << endl;
3395 
3396  // track ID <---> counter
3397  std::map<int, int> mctrackids;
3398  // std::vector<int> newtracks;
3399 
3400  for (size_t ihit = 0; ihit < completeCand->GetNHits(); ihit++) {
3401  PndTrackCandHit candhit = completeCand->GetSortedHit(ihit);
3402  Int_t iHit = candhit.GetHitId();
3403  Int_t detId = candhit.GetDetId();
3404 
3405  FairHit *fhit;
3406  FairMCPoint *fpnt;
3407  if(detId == FairRootManager::Instance()->GetBranchId(fMvdStripBranchName))
3408  {
3409  fhit = (FairHit*) fMvdStripHitArray->At(iHit);
3410  if(!fhit) continue;
3411  Int_t refIndex = fhit->GetRefIndex();
3412  if(refIndex == -1) continue;
3413  fpnt = (FairMCPoint*) fMvdPointArray->At(refIndex);
3414  if(!fpnt) continue;
3415  }
3416  else if(detId == FairRootManager::Instance()->GetBranchId(fMvdPixelBranchName)) {
3417  fhit = (FairHit*) fMvdPixelHitArray->At(iHit);
3418  if(!fhit) continue;
3419  Int_t refIndex = fhit->GetRefIndex();
3420  if(refIndex == -1) continue;
3421  fpnt = (FairMCPoint*) fMvdPointArray->At(refIndex);
3422  if(!fpnt) continue;
3423  }
3424  else if(detId == FairRootManager::Instance()->GetBranchId(fSttBranchName)) {
3425  fhit = (FairHit*) fSttHitArray->At(iHit);
3426  if(!fhit) continue;
3427  Int_t refIndex = fhit->GetRefIndex();
3428  if(refIndex == -1) continue;
3429  fpnt = (FairMCPoint*) fSttPointArray->At(refIndex);
3430  if(!fpnt) continue;
3431  }
3432  else if (detId == FairRootManager::Instance()->GetBranchId(fGemBranchName)) {
3433  fhit = (FairHit*) fGemHitArray->At(iHit);
3434  if(!fhit) continue;
3435  Int_t refIndex = fhit->GetRefIndex();
3436  if(refIndex == -1) continue;
3437  fpnt = (FairMCPoint*) fGemPointArray->At(refIndex);
3438  }
3439  else continue; // CHECK
3440 
3441  int trackID = fpnt->GetTrackID();
3442  cout << "HIT ID " << iHit << " DET ID " << detId << " TRACK ID " << trackID << endl;
3443  // // is the track ID new?
3444  // std::vector<int>::iterator iter;
3445  // iter = find(newtracks.begin(), newtracks.end(), trackID);
3446  // int index = iter - newtracks.begin();
3447  // // if no
3448  // if(index == newtracks.size()) newtracks.push_back(trackID);
3449 
3450  mctrackids[trackID]++;
3451  }
3452 
3453  int counter = 0;
3454  int tmptrackID = -1;
3455  for(size_t itrk = 0; itrk < mctrackids.size(); itrk++) {
3456  if(counter < mctrackids[itrk]) {
3457  counter = mctrackids[itrk];
3458  tmptrackID = itrk;
3459  }
3460  }
3461  completeCand->setMcTrackId(tmptrackID);
3462  cout << "NEW MC " << completeCand->getMcTrackId() << endl;
3463 
3464 }
int getMcTrackId() const
Definition: PndTrackCand.h:60
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
void setMcTrackId(int i)
Definition: PndTrackCand.h:72
TClonesArray * fMvdPixelHitArray
Int_t refIndex
Definition: checkmomentum.C:30
int counter
Definition: ZeeAnalysis.C:59
TClonesArray * fMvdStripHitArray
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
Int_t GetHitId() const
Int_t GetDetId() const
void PndSttMvdGemTracking::UseMonteCarlo ( )
inline

Definition at line 113 of file PndSttMvdGemTracking.h.

References fUseMC.

113 { fUseMC = kTRUE; }
void PndSttMvdGemTracking::WriteHistograms ( )

Definition at line 2099 of file PndSttMvdGemTracking.cxx.

References file, fNPositions, hchosen, hchosen2, hdist, hdist2, hmcdist, hmcx, hmcy, hsigma, and hsigma2.

2099  {
2100 
2101  TFile* file = FairRootManager::Instance()->GetOutFile();
2102  file->cd();
2103  file->mkdir("PndSttMvdGemTracking");
2104  file->cd("PndSttMvdGemTracking");
2105  for(int ipos = 0; ipos < fNPositions; ipos++) {
2106  hdist[ipos]->Write();
2107  delete hdist[ipos];
2108  hdist2[ipos]->Write();
2109  delete hdist2[ipos];
2110  hsigma[ipos]->Write();
2111  delete hsigma[ipos];
2112  hsigma2[ipos]->Write();
2113  delete hsigma2[ipos];
2114  hchosen[ipos]->SetFillColor(3);
2115  hchosen[ipos]->Write();
2116  delete hchosen[ipos];
2117  hchosen2[ipos]->SetFillColor(2);
2118  hchosen2[ipos]->Write();
2119  delete hchosen2[ipos];
2120  hmcdist[ipos]->SetFillColor(4);
2121  hmcdist[ipos]->Write();
2122  delete hmcdist[ipos];
2123  hmcx[ipos]->SetFillColor(4);
2124  hmcx[ipos]->Write();
2125  delete hmcx[ipos];
2126  hmcy[ipos]->SetFillColor(4);
2127  hmcy[ipos]->Write();
2128  delete hmcy[ipos]; }
2129 
2130 
2131 }
TFile * file
Bool_t PndSttMvdGemTracking::ZFind ( Int_t  nhits,
TMatrixT< double >  points,
Double_t  xc,
Double_t  yc,
Double_t  radius 
)

Definition at line 3293 of file PndSttMvdGemTracking.cxx.

References Double_t, fabs(), fSttBranchName, fSttHitArray, fTubeArray, PndSttTube::GetHalfLength(), PndSttTube::GetPosition(), PndSttHit::GetTubeID(), PndSttTube::GetWireDirection(), hit, m, nhits, pos, and CAMath::Sqrt().

3294 {
3295  // cout << "Z FINDER" << endl;
3296  if(nhits == 0) return kFALSE;
3297 
3298  for(int ihit = 0; ihit < nhits; ihit++)
3299  {
3300  Int_t detId = (Int_t) points[ihit][1];
3301  Int_t hitId = (Int_t) points[ihit][0];
3302  // cout << "hitId " << hitId << " detId " << detId << endl;
3303  if(hitId == -1) continue;
3304  Int_t fitflag = (Int_t) points[ihit][10];
3305  if(fitflag != 1) continue;
3306  if(detId != FairRootManager::Instance()->GetBranchId(fSttBranchName)) continue;
3307 
3308  // intersection: tube line with circle trajectory in xy
3309 
3310  // find the tube line
3311  // y = mx + q
3312  PndSttHit *hit = (PndSttHit* ) fSttHitArray->At(hitId);
3313  if(!hit) continue;
3314 
3315  Int_t tubeID = hit->GetTubeID();
3316 
3317  PndSttTube *tube = (PndSttTube* ) fTubeArray->At(tubeID);
3318  TVector3 pos = tube->GetPosition();
3319  TVector3 wireDirection = tube->GetWireDirection();
3320  Double_t halflength = tube->GetHalfLength();
3321  if(wireDirection == TVector3(0., 0., 1.)) continue;
3322 
3323 // pos.Print();
3324 // wireDirection.Print();
3325  // cout << "hl " << halflength << endl;
3326  TVector3 first = pos + wireDirection * halflength; // CHECK
3327  TVector3 second = pos - wireDirection * halflength; // CHECK
3328  // first.Print();
3329 // second.Print();
3330 
3331  double xint, yint, x1, y1, x2, y2, delta;
3332 
3333  // when tube is vertical
3334  if(fabs(second.X() - first.X()) < 1.e-5) {
3335  x1 = first.X();
3336  x2 = x1;
3337 
3338  delta = radius * radius - (x1 - xc) * (x1 - xc);
3339  if(delta < 0) continue;
3340  y1 = yc + TMath::Sqrt(delta);
3341  y2 = yc - TMath::Sqrt(delta);
3342 
3343  }
3344  else {
3345 
3346  Double_t m = (second.Y() - first.Y())/(second.X() - first.X());
3347  Double_t q = first.Y() - m * first.X();
3348 
3349  // center of trajectory xc, yc, radius
3350  delta = (m * (q - yc) - xc) - (m * m + 1) * ((q - yc) * (q - yc) + xc * xc - radius * radius);
3351  if(delta < 0) continue;
3352 
3353  x1 = (- (m * (q - yc) - xc) + delta) / (m * m + 1);
3354  y1 = m * x1 + q;
3355  x2 = (- (m * (q - yc) - xc) - delta) / (m * m + 1);
3356  y2 = m * x2 + q;
3357  }
3358 
3359  double d1 = 0, d2 = 0;
3360  d1 = TMath::Sqrt((y1 - first.Y()) * (y1 - first.Y()) + (x1 - first.X()) * (x1 - first.X()));
3361  d2 = TMath::Sqrt((y2 - first.Y()) * (y2 - first.Y()) + (x2 - first.X()) * (x2 - first.X()));
3362 
3363  if(d1 < d2) {
3364  xint = x1;
3365  yint = y1;
3366  }
3367  else {
3368  xint = x2;
3369  yint = y2;
3370  }
3371 
3372  // APPROXIMATION: take the z of this intersection point // CHECK
3373  Double_t ll = ((xint - first.X()) + (yint - first.Y())) / (wireDirection.X() + wireDirection.Y());
3374  double zint = first.Z() + wireDirection.Z() * ll;
3375 
3376  points[ihit][2] = xint;
3377  points[ihit][3] = yint;
3378  points[ihit][4] = zint;
3379  points[ihit][5] = 1.; // CHECK
3380  points[ihit][6] = 1.; // CHECK
3381  points[ihit][7] = 1.; // CHECK
3382 
3383  // cout << "inters " << xint << " " << yint << " " << zint << endl;
3384 
3385  }
3386 
3387  return kTRUE;
3388 }
TVector3 pos
__m128 m
Definition: P4_F32vec4.h:28
Double_t GetHalfLength()
Definition: PndSttTube.cxx:99
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
Double_t
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
Int_t GetTubeID() const
Definition: PndSttHit.h:75
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndSdsMCPoint * hit
Definition: anasim.C:70
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
Bool_t PndSttMvdGemTracking::ZFit ( TMatrixT< double >  points,
Int_t  charge,
Double_t  xc,
Double_t  yc,
Double_t  radius,
Double_t fitm,
Double_t fip 
)

Definition at line 2727 of file PndSttMvdGemTracking.cxx.

References alpha, CAMath::ATan2(), CalculatePhi(), CompareToPreviousPhi(), CAMath::Cos(), d, Double_t, fGemBranchName, fSttBranchName, nhits, p, phi, CAMath::Sin(), CAMath::Sqrt(), v, x0, and y0.

Referenced by Prefit().

2728 {
2729 
2730  // recalculate z - s fit only from MVD
2731  Double_t Sxx, Sx, Sz, Sxz, S1z;
2732  Double_t Detz = 0.;
2733 
2734  Sx = 0.;
2735  Sz = 0.;
2736  Sxx = 0.;
2737  Sxz = 0.;
2738  S1z = 0.;
2739 
2740  // x0 y0
2741  Double_t d = TMath::Sqrt(xc * xc + yc * yc) - radius;
2742  Double_t phi = TMath::ATan2(yc, xc);
2743 
2744  Double_t x0 = d * TMath::Cos(phi);
2745  Double_t y0 = d * TMath::Sin(phi);
2746 
2747  Double_t Phi0 = TMath::ATan2((y0 - yc),(x0 - xc));
2748  Double_t scos = 0;
2749 
2750  Int_t nhits = points.GetRowUpb() + 1;
2751 
2752  Double_t Fi_pre = 0.;
2753  int zcounter = 0;
2754  for(int ihit = 0; ihit < nhits; ihit++)
2755  {
2756  Int_t detId = (Int_t) points[ihit][1];
2757  Int_t hitId = (Int_t) points[ihit][0];
2758  // cout << "hitId " << hitId << " detId " << detId << endl;
2759  if(hitId == -1) continue;
2760  Int_t fitflag = (Int_t) points[ihit][10];
2761  // if(fitflag == 0 || fitflag == -1) continue; // CHECK // to use zfinder
2762  if(fitflag != 2) continue; // CHECK // to use zfinder
2763 
2764 
2765  if(detId == FairRootManager::Instance()->GetBranchId(fSttBranchName) ||
2766  detId == FairRootManager::Instance()->GetBranchId(fGemBranchName)) continue;
2767 
2768  TVector2 v(x0 - xc, y0 - yc);
2769  Double_t alpha = TMath::ATan2(points[ihit][3] - y0 + radius * TMath::Sin(Phi0), points[ihit][2] - x0 + radius * TMath::Cos(Phi0));
2770  TVector2 p(points[ihit][2] - xc, points[ihit][3] - yc);
2771  Double_t Fi = CalculatePhi(v, p, alpha, Phi0, charge);
2772  if(ihit > 0) Fi = CompareToPreviousPhi(Fi, Fi_pre, charge);
2773  Fi_pre = Fi;
2774  scos = - charge * radius * Fi; // scos = -q * R * phi CHECK :-)GOOD!
2775  // cout << charge << " zfit " << scos << endl;
2776 
2777 
2778  Double_t sigz2 = points[ihit][7] * points[ihit][7]; // CHECK
2779 
2780  if(sigz2 == 0) sigz2 = 1e-5; // CHECK MVD covariance
2781  // cout << "scosl " << scos << " " << points[ihit][4] << " " << sigz2 << " " << points[ihit][7] << endl;
2782  zcounter++;
2783  Sx = Sx + (scos /(sigz2));
2784  Sz = Sz + (points[ihit][4]/(sigz2));
2785  Sxz = Sxz + ((scos * points[ihit][4])/(sigz2));
2786  Sxx = Sxx + ((scos * scos)/(sigz2));
2787  S1z = S1z + 1/(sigz2);
2788 
2789  }
2790 
2791  if(zcounter <= 1) return kFALSE; // CHECK
2792  Detz = S1z*Sxx - Sx*Sx;
2793  if(Detz == 0) {
2794  cout << "DET Z = 0" << endl;
2795  return kFALSE; } // CHECK
2796  fitp = (1/Detz)*(Sxx*Sz - Sx*Sxz);
2797  fitm = (1/Detz)*(S1z*Sxz - Sx*Sz);
2798  // cout << "z fit " << fitm << " " << fitp << endl;
2799 
2800 
2801  return true;
2802 }
Double_t CompareToPreviousPhi(Double_t Fi, Double_t Fi_pre, int charge)
Double_t x0
Definition: checkhelixhit.C:70
Double_t p
Definition: anasim.C:58
TObjArray * d
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
static T Sin(const T &x)
Definition: PndCAMath.h:42
Double_t CalculatePhi(TVector2 v, TVector2 p, double alpha, double Phi0, int charge)
static T Cos(const T &x)
Definition: PndCAMath.h:43
__m128 v
Definition: P4_F32vec4.h:4
Double_t
Double_t y0
Definition: checkhelixhit.C:71
static T ATan2(const T &y, const T &x)
double alpha
Definition: f_Init.h:9

Member Data Documentation

Int_t PndSttMvdGemTracking::countbad
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), and Init().

Int_t PndSttMvdGemTracking::countdoubt
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), and Init().

Int_t PndSttMvdGemTracking::countgood
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), and Init().

Int_t PndSttMvdGemTracking::countnohitonplane[8]
private

Definition at line 179 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), and Init().

Int_t PndSttMvdGemTracking::countplane[8]
private

Definition at line 179 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), and Init().

Int_t PndSttMvdGemTracking::countprinotassigned
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), and Init().

Int_t PndSttMvdGemTracking::countprop[2][8]
private

Definition at line 179 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), Init(), and Retrack().

Int_t PndSttMvdGemTracking::countreconstructablehit
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), and Init().

Int_t PndSttMvdGemTracking::countreconstructablepoint
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances().

Int_t PndSttMvdGemTracking::countsecnotassigned
private

Definition at line 178 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), and Init().

Int_t PndSttMvdGemTracking::countsttmvd
private

Definition at line 179 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), and Init().

Int_t PndSttMvdGemTracking::countsttmvdusable
private

Definition at line 179 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), and Init().

TCanvas* PndSttMvdGemTracking::display
private
TMatrixT<double> PndSttMvdGemTracking::distancemap
private
Int_t PndSttMvdGemTracking::evt
private

Definition at line 179 of file PndSttMvdGemTracking.h.

Referenced by Exec(), and Init().

Double_t PndSttMvdGemTracking::fCombiDistance
private
std::map<int, int> PndSttMvdGemTracking::fCombiMap
private

combimap: hitID <-> 1/0 whether it is combinatorial or not

Definition at line 246 of file PndSttMvdGemTracking.h.

Referenced by AddRemainingHits(), AssignHits(), CheckCombinatorial(), ConsiderCombinatorialEffect(), GetClosestOnFirst(), Reset(), and Retrack().

TClonesArray* PndSttMvdGemTracking::fCompleteTrackArray
private

Output array of PndTracks

Definition at line 164 of file PndSttMvdGemTracking.h.

Referenced by Copy(), Exec(), Init(), and Retrack().

TClonesArray* PndSttMvdGemTracking::fCompleteTrackCandArray
private

Output array of PndTrackCands

Definition at line 162 of file PndSttMvdGemTracking.h.

Referenced by Copy(), EvaluatePerformances(), Exec(), and Init().

Int_t PndSttMvdGemTracking::fDefaultPdgCode
private
Bool_t PndSttMvdGemTracking::fDisplayOn
private
TVector3 PndSttMvdGemTracking::fDj
private

plane of the sensor

Definition at line 194 of file PndSttMvdGemTracking.h.

Referenced by PropagateToGemPlane(), PropagateToGemPlaneAsHelix(), and SetupGEMPlanes().

TVector3 PndSttMvdGemTracking::fDk
private
Bool_t PndSttMvdGemTracking::fEvaluate
private

Definition at line 250 of file PndSttMvdGemTracking.h.

Referenced by Exec(), PndSttMvdGemTracking(), and SetEvaluateFlag().

TString PndSttMvdGemTracking::fGemBranchName
private
TClonesArray* PndSttMvdGemTracking::fGemHitArray
private
PndGemDigiPar* PndSttMvdGemTracking::fGemParameters
private
TClonesArray* PndSttMvdGemTracking::fGemPointArray
private

Input array of GEMPoints

Definition at line 137 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), FillTrueDistances(), Init(), and UpdateMCTrackId().

Double_t PndSttMvdGemTracking::fMaxDistance
private

Definition at line 227 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), PndSttMvdGemTracking(), and SetMaxDistance().

TClonesArray* PndSttMvdGemTracking::fMCTrackArray
private

Input array of mctracks

Definition at line 153 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), FillTrueDistances(), GetPdgFromMC(), Init(), and SetStartParameters().

TString PndSttMvdGemTracking::fMvdPixelBranchName
private
TClonesArray* PndSttMvdGemTracking::fMvdPixelHitArray
private

Input array of MVDHitsPixel

Definition at line 140 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), Init(), Prefit(), SetStartParameters(), and UpdateMCTrackId().

TClonesArray* PndSttMvdGemTracking::fMvdPointArray
private

Input array of MVDPoints

Definition at line 147 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), Init(), SetStartParameters(), and UpdateMCTrackId().

TString PndSttMvdGemTracking::fMvdStripBranchName
private
TClonesArray* PndSttMvdGemTracking::fMvdStripHitArray
private

Input array of MVDHitsStrip

Definition at line 142 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), Init(), Prefit(), SetStartParameters(), and UpdateMCTrackId().

Int_t PndSttMvdGemTracking::fNPositions
private
std::vector<int> PndSttMvdGemTracking::fOrdering
private
std::vector<int>::iterator PndSttMvdGemTracking::fOrderingIterator
private
Int_t PndSttMvdGemTracking::fPdgCode
private

PDG Code

Definition at line 197 of file PndSttMvdGemTracking.h.

Referenced by Exec(), FillTrueDistances(), and PropagateToGemPlane().

Bool_t PndSttMvdGemTracking::fPdgFromMC
private
FairGeanePro* PndSttMvdGemTracking::fPro
private

GEANE propagator

Definition at line 176 of file PndSttMvdGemTracking.h.

Referenced by Init(), PropagateToGemPlane(), and SelectPdgCode().

std::map<int, bool> PndSttMvdGemTracking::fProTracks
private

map which tracks can(not) be extrapolated to GEM

Definition at line 209 of file PndSttMvdGemTracking.h.

Referenced by AddRemainingHits(), CheckCombinatorial(), Exec(), FillTrueDistances(), OnlyOneHitToEachTrack(), and Retrack().

TClonesArray* PndSttMvdGemTracking::fSensPositions
private

position of the planes x, y, z

Definition at line 188 of file PndSttMvdGemTracking.h.

Referenced by PropagateToGemPlane(), PropagateToGemPlaneAsHelix(), and SetupGEMPlanes().

TString PndSttMvdGemTracking::fStartTrackBranchName
private

Definition at line 242 of file PndSttMvdGemTracking.h.

Referenced by Init(), PndSttMvdGemTracking(), and SetTrackBranchNames().

TString PndSttMvdGemTracking::fStartTrackCandBranchName
private

Definition at line 242 of file PndSttMvdGemTracking.h.

Referenced by Init(), PndSttMvdGemTracking(), and SetTrackBranchNames().

TString PndSttMvdGemTracking::fStartTrackIDBranchName
private

Definition at line 242 of file PndSttMvdGemTracking.h.

Referenced by Init(), PndSttMvdGemTracking(), and SetPdgFromMC().

TString PndSttMvdGemTracking::fSttBranchName
private
TClonesArray* PndSttMvdGemTracking::fSttHitArray
private

Input array of STTHits

Definition at line 144 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), Init(), Prefit(), SetStartParameters(), UpdateMCTrackId(), and ZFind().

PndGeoSttPar* PndSttMvdGemTracking::fSttParameters
private

Definition at line 172 of file PndSttMvdGemTracking.h.

Referenced by Init(), and SetParContainers().

TClonesArray* PndSttMvdGemTracking::fSttPointArray
private

Input array of STTPoints

Definition at line 149 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), Init(), SetStartParameters(), and UpdateMCTrackId().

Int_t PndSttMvdGemTracking::fTimes
private

Definition at line 228 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), PndSttMvdGemTracking(), and SetTimes().

TClonesArray* PndSttMvdGemTracking::fTrackArray
private

Input array of mvd + stt tracks

Definition at line 155 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), FillTrueDistances(), and Init().

TClonesArray* PndSttMvdGemTracking::fTrackCandArray
private

Input array of mvd + stt track cand

Definition at line 159 of file PndSttMvdGemTracking.h.

Referenced by Init().

TClonesArray* PndSttMvdGemTracking::fTrackIDArray
private

Input array of mvd + stt trackID

Definition at line 157 of file PndSttMvdGemTracking.h.

Referenced by GetPdgFromMC(), and Init().

TClonesArray* PndSttMvdGemTracking::fTubeArray
private

from parameters array of PndSttTube

Definition at line 167 of file PndSttMvdGemTracking.h.

Referenced by Init(), IntersectionFinder(), Prefit(), and ZFind().

Int_t PndSttMvdGemTracking::fTurn
private

Definition at line 229 of file PndSttMvdGemTracking.h.

Referenced by Exec(), IsAssignable(), PndSttMvdGemTracking(), and Retrack().

Bool_t PndSttMvdGemTracking::fUseMC
private

Definition at line 239 of file PndSttMvdGemTracking.h.

Referenced by PndSttMvdGemTracking(), SetStartParameters(), and UseMonteCarlo().

TH2F* PndSttMvdGemTracking::h[8]
private

Definition at line 234 of file PndSttMvdGemTracking.h.

Referenced by Reset().

TH1F * PndSttMvdGemTracking::hchosen[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), SetupGEMPlanes(), and WriteHistograms().

TH1F * PndSttMvdGemTracking::hchosen2[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), SetupGEMPlanes(), and WriteHistograms().

TH1F* PndSttMvdGemTracking::hdist[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), SetupGEMPlanes(), and WriteHistograms().

TH1F * PndSttMvdGemTracking::hdist2[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), SetupGEMPlanes(), and WriteHistograms().

TMatrixT<float> PndSttMvdGemTracking::hitcounter
private

hitcounter: vector of number of hits for each position in ordering (e.g. 4 hits on 1st plane, 0 on 2nd, 3 on 3rd, ...)

Definition at line 221 of file PndSttMvdGemTracking.h.

Referenced by AddRemainingHits(), AssignHits(), GetClosestOnFirst(), GetHitsAssociatedToTrackOnPlane(), OrderGemHits(), Reset(), and SetupGEMPlanes().

TMatrixT<float> PndSttMvdGemTracking::hitmap
private

hitmap: rows = position in ordering (plane: 0, 1, 2, ...) cols = number of hits in that position (if I have 3 hits on the i-th plane there will be 3 cols filled for row i-th) ... filled with hitID of each hit

Definition at line 216 of file PndSttMvdGemTracking.h.

Referenced by AddRemainingHits(), AssignHits(), GetClosestOnFirst(), GetHitsAssociatedToTrackOnPlane(), OrderGemHits(), and Reset().

TH1F * PndSttMvdGemTracking::hmcdist[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), SetupGEMPlanes(), and WriteHistograms().

TH1F * PndSttMvdGemTracking::hmcx[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), SetupGEMPlanes(), and WriteHistograms().

TH1F * PndSttMvdGemTracking::hmcy[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by FillTrueDistances(), SetupGEMPlanes(), and WriteHistograms().

TH2F * PndSttMvdGemTracking::hnotskewed
private

Definition at line 234 of file PndSttMvdGemTracking.h.

TH1F * PndSttMvdGemTracking::hsigma[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), SetupGEMPlanes(), and WriteHistograms().

TH1F * PndSttMvdGemTracking::hsigma2[8]
private

Definition at line 235 of file PndSttMvdGemTracking.h.

Referenced by IsAssignable(), SetupGEMPlanes(), and WriteHistograms().

TH2F * PndSttMvdGemTracking::hskewed
private

Definition at line 234 of file PndSttMvdGemTracking.h.

std::vector<int> PndSttMvdGemTracking::notassignedhits
private
std::vector<int> PndSttMvdGemTracking::notassignedtracks
private

Definition at line 203 of file PndSttMvdGemTracking.h.

std::map<int, bool> PndSttMvdGemTracking::towhichplane
private

hit to plane map

Definition at line 206 of file PndSttMvdGemTracking.h.

Referenced by AddRemainingHits(), Exec(), OrderGemHits(), Reset(), and Retrack().

std::vector<int> PndSttMvdGemTracking::trackindexes
private

Definition at line 201 of file PndSttMvdGemTracking.h.

Referenced by CountTracks(), Exec(), GetTrackIndex(), and Reset().

std::vector<std::pair<int, int> > PndSttMvdGemTracking::trackvector
private
std::vector<int> PndSttMvdGemTracking::usabletracks
private

Definition at line 181 of file PndSttMvdGemTracking.h.

Referenced by EvaluatePerformances(), Exec(), and Reset().


The documentation for this class was generated from the following files: