FairRoot/PandaRoot
PndPidMdtInfo.cxx
Go to the documentation of this file.
1 #include "PndPidCorrelator.h"
2 #include "PndMdtHit.h"
3 #include "PndMdtPoint.h"
4 #include "PndMdtTrk.h"
5 
6 #include <cmath>
7 
8 
9 //_________________________________________________________________
11  //---
12  FairTrackParP par = track->GetParamLast();
13  Int_t ierr = 0;
14  FairTrackParH *helix = new FairTrackParH(&par, ierr);
15 
16  map<Int_t, Int_t>mapMdtTrk;
17 
18  if (fMdtMode == 3)
19  {
20  for (Int_t tt = 0; tt<fMdtTrk->GetEntriesFast(); tt++)
21  {
22  PndMdtTrk *mdtTrk = (PndMdtTrk*)fMdtTrk->At(tt);
23  mapMdtTrk[mdtTrk->GetHitIndex(0)] = tt;
24  }
25  }
26  PndMdtHit *mdtHit = NULL;
27  Int_t mdtEntries = fMdtHit->GetEntriesFast();
28  Int_t mdtIndex = -1, mdtMod = 0, mdtLayer = 0, mdtHits = 0;
29  Float_t mdtGLength = -1000;
30  Float_t mdtQuality = 1000000;
31  Float_t mdtIron = 0., mdtMom = 0, mdtTempMom = 0;
32 
33  //Float_t chi2 = 0; //[R.K. 01/2017] unused variable
34  TVector3 vertex(0., 0., 0.);
35  TVector3 vertexD(0., 0., 0.);
36  TVector3 mdtPos(0., 0., 0.);
37  TVector3 momentum(0., 0., 0.);
38  for (Int_t mm = 0; mm<mdtEntries; mm++)
39  {
40  mdtHit = (PndMdtHit*)fMdtHit->At(mm);
41  if ( fIdeal && ( ((PndMdtPoint*)fMdtPoint->At(mdtHit->GetRefIndex()))->GetTrackID() !=pidCand->GetMcIndex()) ) continue;
42  if (mdtHit->GetLayerID()!=0) continue;
43  if (mdtHit->GetModule()>2) continue;
44  mdtHit->Position(mdtPos);
45  if (fGeanePro) // Overwrites vertex if Geane is used
46  {
47  fGeanePropagator->SetPoint(mdtPos);
48  fGeanePropagator->PropagateToPCA(1, 1);
49  vertex.SetXYZ(-10000, -10000, -10000); // reset vertex
50  vertexD.SetXYZ(-10000, -10000, -10000); // reset vertex
51  FairTrackParH *fRes= new FairTrackParH();
52  Bool_t rc = fGeanePropagator->Propagate(helix, fRes, fPidHyp*pidCand->GetCharge());
53  if (!rc) continue;
54  mdtTempMom = fRes->GetMomentum().Mag();
55  vertex.SetXYZ(fRes->GetX(), fRes->GetY(), fRes->GetZ());
56  vertexD.SetXYZ(fRes->GetDX(), fRes->GetDY(), fRes->GetDZ());
57  mdtGLength = fGeanePropagator->GetLengthAtPCA();
58  }
59 
60  Float_t dist;
61  if (mdtHit->GetModule()==1)
62  {
63  dist = (mdtPos-vertex).Mag2();
64  }
65  else
66  {
67  dist = (vertex.X()-mdtPos.X())*(vertex.X()-mdtPos.X())+(vertex.Y()-mdtPos.Y())*(vertex.Y()-mdtPos.Y());
68  }
69 
70  if ( mdtQuality > dist)
71  {
72  mdtIndex = mm;
73  mdtQuality = dist;
74  mdtMod = mdtHit->GetModule();
75  mdtMom = mdtTempMom;
76  mdtLayer = 1;
77  if (fMdtMode==3)
78  {
79  PndMdtTrk *mdtTrk = (PndMdtTrk*)fMdtTrk->At(mapMdtTrk[mdtIndex]);
80  mdtIndex = mapMdtTrk[mm];
81  mdtLayer = mdtTrk->GetLayerCount();
82  mdtIron = mdtTrk->GetIronDist();
83  mdtMod = mdtTrk->GetModule();
84  mdtHits = 0;
85  for (Int_t iLayer=0; iLayer<mdtLayer; iLayer++)
86  {
87  mdtHits = mdtHits + mdtTrk->GetHitMult(iLayer);
88  //std::cout << iLayer << "\t" << mdtTrk->GetHitMult(iLayer) << "\t" << mdtHits << std::endl;
89  }
90  }
91  }
92  if (fDebugMode)
93  {
94  Float_t ntuple[] = {static_cast<Float_t>(vertex.X()), static_cast<Float_t>(vertex.Y()), static_cast<Float_t>(vertex.Z()),
95  static_cast<Float_t>(vertexD.X()), static_cast<Float_t>(vertexD.Y()), static_cast<Float_t>(vertexD.Z()), static_cast<Float_t>(vertex.Phi()),
96  static_cast<Float_t>(helix->GetMomentum().Mag()), static_cast<Float_t>(helix->GetQ()), static_cast<Float_t>(helix->GetMomentum().Theta()), static_cast<Float_t>(helix->GetZ()),
97  static_cast<Float_t>(mdtPos.X()), static_cast<Float_t>(mdtPos.Y()), static_cast<Float_t>(mdtPos.Z()), static_cast<Float_t>(mdtPos.Phi()),mdtTempMom,
98  dist, static_cast<Float_t>(mdtHit->GetModule()), static_cast<Float_t>(vertex.DeltaPhi(mdtPos)), mdtGLength, static_cast<Float_t>(mdtLayer), static_cast<Float_t>(mdtHits) };
99  // Float_t ntuple[] = {vertex.X(), vertex.Y(), vertex.Z(),
100  // vertexD.X(), vertexD.Y(), vertexD.Z(), vertex.Phi(),
101  // helix->GetMomentum().Mag(), helix->GetQ(), helix->GetMomentum().Theta(), helix->GetZ(),
102  // mdtPos.X(), mdtPos.Y(), mdtPos.Z(), mdtPos.Phi(), mdtTempMom,
103  // dist, mdtHit->GetModule(), vertex.DeltaPhi(mdtPos), mdtGLength, mdtLayer, mdtHits};
104  mdtCorr->Fill(ntuple);
105  }
106  }
107 
108  if ((mdtQuality<fCorrPar->GetMdtCut()) || ( fIdeal && mdtIndex!=-1))
109  {
110  pidCand->SetMuoIndex(mdtIndex);
111  pidCand->SetMuoQuality(mdtQuality);
112  pidCand->SetMuoIron(mdtIron);
113  pidCand->SetMuoMomentumIn(mdtMom);
114  pidCand->SetMuoModule(mdtMod);
115  pidCand->SetMuoNumberOfLayers(mdtLayer);
116  pidCand->SetMuoHits(mdtHits);
117  }
118 
119  // if (fMdtRefit && (mdtIndex!=-1) && (mdtMom>0.) )
120  // {
121  // PndMdtTrk *mdtTrk = (PndMdtTrk*)fMdtTrk->At(mdtIndex);
122  // PndTrack *mdtTrack = new PndTrack(*track);
123  // PndTrackCand *oldCand = track->GetTrackCandPtr();
124  // PndTrackCand *newCand = mdtTrk->AddTrackCand(oldCand);
125  // mdtTrack->SetTrackCand(*newCand);
126  // Int_t fCharge= mdtTrack->GetParamFirst().GetQ();
127  // Int_t PDGCode = fPidHyp*fCharge;
128 
129  // PndTrack *fitTrack = new PndTrack();
130  // fitTrack = fFitter->Fit(mdtTrack, PDGCode);
131  // PndTrack* pndTrack = new PndTrack(fitTrack->GetParamFirst(), fitTrack->GetParamLast(), fitTrack->GetTrackCand(),
132  // fitTrack->GetFlag(), fitTrack->GetChi2(), fitTrack->GetNDF(), fitTrack->GetPidHypo(), fitTrack->GetRefIndex(), kLheTrack);
133  // AddMdtTrack(pndTrack);
134  // }
135  return kTRUE;
136 }
137 
138 //_________________________________________________________________
140 
141  Int_t nHits = fMdtHit->GetEntriesFast();
142  if (nHits==0) return kFALSE;
143 
144  PndMdtHit *mdtHit = NULL;
145  for (Int_t iHit=0; iHit<nHits; iHit++)
146  {
147  mdtHit = (PndMdtHit*) fMdtHit->At(iHit);
148  Int_t mdtLayer = -1;// mdtModule = -1, //[R.K.03/2017] unused variable
149  //std::cout << mdtHit->GetModule() << "\t" << mdtHit->GetLayerID() << "\t" << iHit << std::endl;
150  switch (mdtHit->GetModule())
151  {
152  case 1:
153  //mdtModule = 1; //[R.K.03/2017] unused variable
154  mdtLayer = mdtHit->GetLayerID();
155  mapMdtBarrel[mdtLayer].push_back(iHit);
156  break;
157 
158  case 2:
159  //mdtModule = 2; //[R.K.03/2017] unused variable
160  mdtLayer = mdtHit->GetLayerID();
161  mapMdtEndcap[mdtLayer].push_back(iHit);
162  break;
163 
164  case 3:
165  //mdtModule = 2; //[R.K.03/2017] unused variable
166  mdtLayer = mdtHit->GetLayerID()+5;
167  mapMdtEndcap[mdtLayer].push_back(iHit);
168  break;
169 
170  case 4:
171  //mdtModule = 3; //[R.K.03/2017] unused variable
172  mdtLayer = mdtHit->GetLayerID();
173  mapMdtForward[mdtLayer].push_back(iHit);
174  break;
175 
176  default:
177  std::cout << "-E- PndPidCorrelator::MdtMapping: Wrong MDT Module" << std::endl;
178  return kFALSE;
179  }
180  }
181 
182  return kTRUE;
183 }
184 
185 //_________________________________________________________________
187 {
188  //fVerbose = 3;
189  // Thichness of barrel iron layers [cm]
190  mdtIronThickness[0][0] = 6.;
191  for (Int_t ll=1; ll<12;ll++) mdtIronThickness[0][ll] = 3;
192  mdtIronThickness[0][12] = 6.;
193  // Thichness of endcap+mf iron layers [cm]
194  for (Int_t ll=0; ll<10;ll++) mdtIronThickness[1][ll] = 6;
195  // Thichness of forward iron layers [cm]
196  for (Int_t ll=0; ll<17;ll++) mdtIronThickness[2][ll] = 6;
197 
198  Text_t buffer[250];
199 
200  for (Int_t ll=0; ll<5; ll++)
201  {
202  sprintf(buffer,"cave_1/Mdt_1/MdtEndcap_1/MdtEndcapLayer0%i_1/MDT2s0l%ib0w0_%i",ll,ll,200+(8*ll));
203  gGeoManager->cd(buffer);
204  Double_t local[3] = {0., 0., 0.};
205  Double_t master[3];
206  gGeoManager->LocalToMaster(local, master);
207  mdtLayerPos[1][ll] = master[2];
208  if (fVerbose>1) std::cout << "mdtLayerPos[1][" << ll << "] = " << mdtLayerPos[1][ll] << ";" << std::endl;
209  }
210  for (Int_t ll=0; ll<6; ll++)
211  {
212  sprintf(buffer,"cave_1/Mdt_1/MdtMuonFilter_1/MdtMuonFilterLayer0%i_1/MDT3s0l%ib0w0_%i",ll,ll,300+(8*ll));
213  gGeoManager->cd(buffer);
214  Double_t local[3] = {0., 0., 0.};
215  Double_t master[3];
216  gGeoManager->LocalToMaster(local, master);
217  mdtLayerPos[1][ll+5] = master[2];
218  if (fVerbose>1) std::cout << "mdtLayerPos[1][" << ll+5 << "] = " << mdtLayerPos[1][ll+5] << ";" << std::endl;
219  }
220  for (Int_t ll=0; ll<17; ll++)
221  {
222  sprintf(buffer,"cave_1/Mdt_1/MdtForward_1/MDT4s0l%ib0w0_%i",ll,ll);
223  gGeoManager->cd(buffer);
224  Double_t local[3] = {0., 0., 0.};
225  Double_t master[3];
226  gGeoManager->LocalToMaster(local, master);
227  mdtLayerPos[2][ll] = master[2];
228  if (fVerbose>1) std::cout << "mdtLayerPos[2][" << ll << "] = " << mdtLayerPos[2][ll] << ";" << std::endl;
229  }
230 
231  return kTRUE;
232 }
233 
234 //_________________________________________________________________
235 Bool_t PndPidCorrelator::GetMdt2Info(FairTrackParH* , PndPidCandidate* ) { // helix pidCand //[R.K.03/2017] unused variable(s)
236 
237  return kTRUE;
238 }
239 
240 //_________________________________________________________________
241 Bool_t PndPidCorrelator::GetFMdtInfo(FairTrackParP* helix, PndPidCandidate* pidCand) {
242 
243  if(!helix)
244  {
245  std::cerr << "<Error> PndPidCorrelator::GetFMdtInfo: FairTrackParH NULL pointer parameter."<<std::endl;
246  return kFALSE;
247  }
248  if(!pidCand)
249  {
250  std::cerr << "<Error> PndPidCorrelator::GetFMdtInfo: pidCand NULL pointer parameter."<<std::endl;
251  return kFALSE;
252  }
253 
254  if(helix->GetZ() < fCorrPar->GetZLastPlane())
255  {
256  if (fVerbose>0) std::cout << "-W- PndPidCorrelator::GetFMdtInfo: Skipping tracks not reaching the last FTS layer" << std::endl;
257  return kFALSE;
258  }
259 
260  if(helix->GetPz() <= 0.)
261  {
262  std::cout << "-W- PndPidCorrelator::GetFMdtInfo: Skipping tracks going backward" << std::endl;
263  return kFALSE;
264  }
265 
266  FairGeanePro *fProMdt = new FairGeanePro();
267  if (!fCorrErrorProp) fProMdt->PropagateOnlyParameters();
268 
269  PndMdtHit *mdtHit = NULL;
270 
271  Int_t mdtIndex = -1;
272  Float_t mdtGLength = -1000;
273  Float_t mdtQuality = 1000000;
274 
275  //Float_t chi2 = 0; //[R.K. 01/2017] unused variable
276  TVector3 mdtPos(0., 0., 0.);
277  TVector3 momentum(0., 0., 0.);
278  TVector3 vertex(0., 0., 0.);
279  TVector3 vertexD(0., 0., 0.);
280  Float_t prop0Z = mdtLayerPos[2][0] ;
281  Float_t prop0X = helix->GetX() + (prop0Z - helix->GetZ()) * helix->GetPx() / helix->GetPz();
282  Float_t prop0Y = helix->GetY() + (prop0Z - helix->GetZ()) * helix->GetPy() / helix->GetPz();
283 
284  vertex.SetXYZ(prop0X, prop0Y, prop0Z);
285  mdtGLength = (vertex-helix->GetPosition()).Mag();
286 
287  if (fGeanePro && !fIdeal) // Overwrites vertex if Geane is used
288  {
289  TVector3 jj(1,0,0), kk(0,1,0);
290  fProMdt->PropagateFromPlane(jj, kk);
291  fProMdt->PropagateToPlane(vertex, jj, kk);
292  //fProMdt->PropagateToVolume("MDT4s0l0b0w0",0,1);
293  FairTrackParP *fRes= new FairTrackParP();
294  Bool_t rc = fProMdt->Propagate(helix, fRes, fPidHyp*pidCand->GetCharge());
295  if (rc)
296  {
297  mdtGLength = fProMdt->GetLengthAtPCA();
298  vertex.SetXYZ(fRes->GetX(), fRes->GetY(), fRes->GetZ());
299  momentum.SetXYZ(fRes->GetPx(), fRes->GetPy(), fRes->GetPz());
300  vertexD.SetXYZ(fRes->GetDX(), fRes->GetDY(), fRes->GetDZ());
301  }
302  else
303  {
304  std::cout << "-W- PndPidCorrelator::GetFMdtInfo: Propagation failed" << std::endl;
305  }
306  }
307 
308  vector<Int_t>vecMdt0 = mapMdtForward[0]; // hit in layer0
309  TVector3 oldPos(0., 0., 0.);
310  for (size_t mm0 = 0; mm0< vecMdt0.size(); mm0++)
311  {
312  mdtHit = (PndMdtHit*)fMdtHit->At(vecMdt0[mm0]);
313  if ( fIdeal && ( ((PndMdtPoint*)fMdtPoint->At(mdtHit->GetRefIndex()))->GetTrackID() !=pidCand->GetMcIndex()) ) continue;
314  mdtHit->Position(mdtPos);
315  Float_t dist = (mdtPos-vertex).Mag2();
316 
317  if ( mdtQuality > dist)
318  {
319  mdtIndex = mm0;
320  mdtQuality = dist;
321  oldPos = mdtPos;
322  }
323 
324  if (fDebugMode)
325  {
326  Int_t mdtLayer = 0, mdtHits = 1; // dummy values to fill the ntuple
327  Float_t ntuple[] = {static_cast<Float_t>(vertex.X()), static_cast<Float_t>(vertex.Y()), static_cast<Float_t>(vertex.Z()),
328  static_cast<Float_t>(vertexD.X()), static_cast<Float_t>(vertexD.Y()), static_cast<Float_t>(vertexD.Z()),
329  static_cast<Float_t>(vertex.Phi()),
330  static_cast<Float_t>(helix->GetMomentum().Mag()), static_cast<Float_t>(helix->GetQ()), static_cast<Float_t>(helix->GetMomentum().Theta()), static_cast<Float_t>(helix->GetZ()),
331  static_cast<Float_t>(mdtPos.X()), static_cast<Float_t>(mdtPos.Y()), static_cast<Float_t>(mdtPos.Z()), static_cast<Float_t>(mdtPos.Phi()), static_cast<Float_t>(momentum.Mag()),
332  dist, static_cast<Float_t>(mdtHit->GetModule()), static_cast<Float_t>(vertex.DeltaPhi(mdtPos)), mdtGLength, static_cast<Float_t>(mdtLayer), static_cast<Float_t>(mdtHits) };
333  // Float_t ntuple[] = {vertex.X(), vertex.Y(), vertex.Z(),
334  // vertexD.X(), vertexD.Y(), vertexD.Z(),
335  // vertex.Phi(),
336  // helix->GetMomentum().Mag(), helix->GetQ(), helix->GetMomentum().Theta(), helix->GetZ(),
337  // mdtPos.X(), mdtPos.Y(), mdtPos.Z(), mdtPos.Phi(), momentum.Mag(),
338  // dist, mdtHit->GetModule(), vertex.DeltaPhi(mdtPos), mdtGLength, mdtLayer, mdtHits};
339  mdtCorr->Fill(ntuple);
340  }
341  } // end of layer0 loop
342 
343  if ( (mdtQuality<fCorrPar->GetFMdtCut()) || (fIdeal && mdtIndex!=-1) )
344  {
345  pidCand->SetMuoQuality(mdtQuality);
346  pidCand->SetMuoIndex(mdtIndex);
347  pidCand->SetMuoModule(4);
348  pidCand->SetMuoMomentumIn(momentum.Mag());
349 
350  // MUON TRACKING
351  Int_t layerCount = 0, lastLayer = 0, mdtHits = 1;
352  Float_t ironDist = 0.;
353 
354  for (Int_t ll=1; ll<17; ll++)
355  {
356  vector<Int_t>vecMdt = mapMdtForward[ll]; // hit in layer ll
357  if (vecMdt.size()==0) continue;
358 
359  // extrapolation w/o genfit
360  Float_t propZ = mdtLayerPos[2][ll] ;
361  Float_t propX = helix->GetX() + (propZ - helix->GetZ()) * helix->GetPx() / helix->GetPz();
362  Float_t propY = helix->GetY() + (propZ - helix->GetZ()) * helix->GetPy() / helix->GetPz();
363  vertex.SetXYZ(propX, propY, propZ);
364 
365  Float_t corrDist = -1, corrHitDist = -1; //mdtLQuality = 1000000, //[R.K. 01/2017] unused variable
366  Int_t layerMult = 0; // mdtLIndex = -1, hitLCounts = 0, //[R.K. 01/2017] unused variable
367  TVector3 corrPos(0., 0., 0.);
368  for (size_t mm = 0; mm< vecMdt.size(); mm++)
369  {
370  mdtHit = (PndMdtHit*)fMdtHit->At(vecMdt[mm]);
371  if ( fIdeal && ( ((PndMdtPoint*)fMdtPoint->At(mdtHit->GetRefIndex()))->GetTrackID() !=pidCand->GetMcIndex()) ) continue;
372  mdtHit->Position(mdtPos);
373  Float_t dist = (mdtPos-vertex).Mag2();
374  Float_t hitDist = (mdtPos-oldPos).Mag2();
375  if ( (corrDist<0) || (dist<corrDist) )
376  {
377  //mdtLIndex = vecMdt[mm]; //[R.K.03/2017] unused variable
378  corrDist = dist;
379  corrHitDist = hitDist;
380  corrPos = mdtPos;
381  }
382  if ( (dist>0.) && (TMath::Sqrt(dist) < 20) ) layerMult++;
383  } // end of mdtHit loop per layer
384 
385  if ( (corrDist>0.) && (TMath::Sqrt(corrDist) < 20) )
386  {
387  ironDist = ironDist + mdtIronThickness[2][ll] * TMath::Sqrt(corrHitDist)/(mdtLayerPos[2][ll] - mdtLayerPos[2][lastLayer]);
388  mdtHits = mdtHits + layerMult;
389  layerCount++;
390  lastLayer = ll;
391  oldPos = corrPos; // reset position for next mdt layer
392  }
393 
394  } // end of internal layers loop
395 
396  //pidCand->SetMuoNumberOfLayers(layerCount);
397  pidCand->SetMuoNumberOfLayers(lastLayer);
398  pidCand->SetMuoIron(ironDist);
399  pidCand->SetMuoHits(mdtHits);
400  }
401 
402  return kTRUE;
403 }
404 
Bool_t GetMdt2Info(FairTrackParH *helix, PndPidCandidate *pid)
Float_t GetIronDist() const
Definition: PndMdtTrk.h:42
int fVerbose
Definition: poormantracks.C:24
Int_t GetCharge() const
Short_t GetLayerID() const
Definition: PndMdtHit.h:34
void GetHitMult(Int_t *hit)
Definition: PndMdtTrk.h:34
map< Int_t, vector< Int_t > > mapMdtBarrel
TClonesArray * fMdtHit
PndMdtPoint TCA.
PndRiemannTrack track
Definition: RiemannTest.C:33
Bool_t GetFMdtInfo(FairTrackParP *helix, PndPidCandidate *pid)
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
Float_t mdtLayerPos[3][20]
Double_t par[3]
PndPidCorrPar * fCorrPar
PndRichHit TCA.
Bool_t GetMdtInfo(PndTrack *track, PndPidCandidate *pid)
void SetMuoNumberOfLayers(Int_t val)
TGeoManager * gGeoManager
Int_t GetMcIndex() const
TClonesArray * fMdtTrk
PndMdtHit TCA.
Int_t GetModule() const
Definition: PndMdtTrk.h:48
void SetMuoQuality(Double_t val)
int nHits
Definition: RiemannTest.C:16
Double_t
void SetMuoIron(Double_t val)
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
Int_t GetLayerCount() const
Definition: PndMdtTrk.h:44
map< Int_t, vector< Int_t > > mapMdtForward
FairGeanePro * fGeanePropagator
Refitter for MDT tracks.
TClonesArray * fMdtPoint
PndEmcDigi TCA.
void SetMuoMomentumIn(Double_t val)
Float_t mdtIronThickness[3][20]
Int_t GetHitIndex(Int_t lay) const
Definition: PndMdtTrk.h:37
Double_t const mm
void SetMuoModule(Int_t val)
Short_t GetModule() const
Definition: PndMdtHit.h:32
ClassImp(PndAnaContFact)
void SetMuoIndex(Int_t val)
map< Int_t, vector< Int_t > > mapMdtEndcap
Float_t GetZLastPlane()
Definition: PndPidCorrPar.h:22
void SetMuoHits(Int_t val)