FairRoot/PandaRoot
PndLmdBPtestTask.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndLmdBPtestTask source file -----
3 // ----- Created 12/12/13 by A.Karavdina -----
4 // -------------------------------------------------------------------------
5 // libc includes
6 #include <iostream>
7 
8 // Root includes
9 #include <TMatrixDSym.h>
10 #include "TClonesArray.h"
11 #include "TROOT.h"
12 #include "TTree.h"
13 #include "TVector3.h"
14 
15 // framework includes
16 #include "FairBaseParSet.h"
17 #include "FairHit.h"
18 #include "FairRootManager.h"
19 #include "FairRun.h"
20 #include "FairRunAna.h"
21 #include "FairRuntimeDb.h"
22 #include "FairTrackParH.h"
23 #include "FairTrackParP.h"
24 #include "PndLmdBPtestTask.h"
25 #include "PndMCTrack.h"
26 #include "PndMultiField.h"
27 #include "PndTrack.h"
28 #include "TDatabasePDG.h"
29 #include "TGeant3.h"
30 // PndSds includes
31 #include "PndSdsHit.h"
32 #include "PndSdsMCPoint.h"
33 #include "PndSdsMergedHit.h"
34 // PndLmd includes
35 #include <map>
36 #include <vector>
37 #include "PndLinTrack.h"
38 
39 // ----- Default constructor -------------------------------------------
40 PndLmdBPtestTask::PndLmdBPtestTask()
41  : FairTask("test with Geane Task for PANDA Lmd"),
42  fEventNr(0),
43  fUseMVDPoint(false) {
44  // tprop = new TNtuple();
45 }
46 // -------------------------------------------------------------------------
47 
48 PndLmdBPtestTask::PndLmdBPtestTask(Double_t pBeam, TVector3 IP)
49  : FairTask("Geane Task for PANDA Lmd"), fEventNr(0), fUseMVDPoint(false) {
50  fPDGid = -2212;
51  fPbeam = pBeam;
52  cout << "Beam Momentum for particle with PDGid#" << fPDGid << " this run is "
53  << fPbeam << endl;
54  vtx = IP;
55  cout << "Interaction Point:" << endl;
56  vtx.Print();
57  // tprop = new TNtuple("tprop","forward MC vs. backward
58  // rec","xrec:yrec:zrec:prec:thetarec:phirec:xmc:ymc:zmc:pmc:thetamc:phimc") ;
59  tprop = new TTree("tprop",
60  "forward MC vs. backward rec with Mag.Field comp@MC point");
61  tprop->Branch("xrec", &fxrec);
62  tprop->Branch("yrec", &fyrec);
63  tprop->Branch("zrec", &fzrec);
64  tprop->Branch("xmc", &fxmc);
65  tprop->Branch("ymc", &fymc);
66  tprop->Branch("zmc", &fzmc);
67  tprop->Branch("xmclmd", &fxmclmd);
68  tprop->Branch("ymclmd", &fymclmd);
69  tprop->Branch("zmclmd", &fzmclmd);
70  tprop->Branch("thetarec", &fthetarec);
71  tprop->Branch("phirec", &fphirec);
72  tprop->Branch("prec", &fprec);
73  tprop->Branch("thetamc", &fthetamc);
74  tprop->Branch("phimc", &fphimc);
75  tprop->Branch("pmc", &fpmc);
76  tprop->Branch("thetamclmd", &fthetamclmd);
77  tprop->Branch("phimclmd", &fphimclmd);
78  tprop->Branch("pmclmd", &fpmclmd);
79  tprop->Branch("xrec_err", &fxrec_err);
80  tprop->Branch("yrec_err", &fyrec_err);
81  tprop->Branch("xmc_err", &fxmc_err);
82  tprop->Branch("ymc_err", &fymc_err);
83  tprop->Branch("xmclmd_err", &fxmclmd_err);
84  tprop->Branch("ymclmd_err", &fymclmd_err);
85 
86  tprop->Branch("vrec", &fvrec);
87  tprop->Branch("wrec", &fwrec);
88  tprop->Branch("tvrec", &ftvrec);
89  tprop->Branch("twrec", &ftwrec);
90  tprop->Branch("vrec_err", &fvrec_err);
91  tprop->Branch("wrec_err", &fwrec_err);
92  tprop->Branch("tvrec_err", &ftvrec_err);
93  tprop->Branch("twrec_err", &ftwrec_err);
94  tprop->Branch("vmc", &fvmc);
95  tprop->Branch("wmc", &fwmc);
96  tprop->Branch("tvmc", &ftvmc);
97  tprop->Branch("twmc", &ftwmc);
98  tprop->Branch("vmc_err", &fvmc_err);
99  tprop->Branch("wmc_err", &fwmc_err);
100  tprop->Branch("tvmc_err", &ftvmc_err);
101  tprop->Branch("twmc_err", &ftwmc_err);
102  tprop->Branch("twmclmd_err", &ftwmclmd_err);
103  tprop->Branch("vmclmd", &fvmclmd);
104  tprop->Branch("wmclmd", &fwmclmd);
105  tprop->Branch("tvmclmd", &ftvmclmd);
106  tprop->Branch("twmclmd", &ftwmclmd);
107  tprop->Branch("vmclmd_err", &fvmclmd_err);
108  tprop->Branch("wmclmd_err", &fwmclmd_err);
109  tprop->Branch("tvmclmd_err", &ftvmclmd_err);
110  tprop->Branch("twmclmd_err", &ftwmclmd_err);
111 
112  tprop->Branch("Bx", &fbx);
113  tprop->Branch("By", &fby);
114  tprop->Branch("Bz", &fbz);
115 }
116 
117 // ----- Destructor ----------------------------------------------------
118 PndLmdBPtestTask::~PndLmdBPtestTask() {}
119 
120 // ----- Public method Init --------------------------------------------
121 InitStatus PndLmdBPtestTask::Init() {
122  // Get RootManager
123  FairRootManager *ioman = FairRootManager::Instance();
124  if (!ioman) {
125  std::cout << "-E- PndLmdBPtestTask::Init: "
126  << "RootManager not instantiated!" << std::endl;
127  return kFATAL;
128  }
129 
130  fMCHits = (TClonesArray *)ioman->GetObject("LMDPoint");
131  if (!fMCHits) {
132  std::cout << "-W- PndLmdBPtestTask::Init: "
133  << "No LMDPoint"
134  << " array!" << std::endl;
135  return kERROR;
136  }
137 
138  // Get MC tracks
139  fMCTracks = (TClonesArray *)ioman->GetObject("MCTrack");
140  if (!fMCTracks) {
141  std::cout << "-W- PndLmdTrkQTask::Init: "
142  << "No MCTrack"
143  << " array!" << std::endl;
144  return kERROR;
145  }
146 
147  // fTracks = (TClonesArray*) ioman->GetObject("LMDTrack");
148  fTracks = (TClonesArray *)ioman->GetObject("LMDPndTrackFilt");
149  if (!fTracks) {
150  std::cout << "-W- PndLmdBPtestTask::Init: "
151  << "No Track"
152  << " array!" << std::endl;
153  return kERROR;
154  }
155 
156  // Get trk cand [needed only for tests!!!]
157  fRecCandTracks = (TClonesArray *)ioman->GetObject("LMDTrackCand");
158  if (!fRecCandTracks) {
159  std::cout << "-W- PndLmdBPtestTask::Init: "
160  << "No LMDTrackCand [needed only for tests!!!]"
161  << " array!" << std::endl;
162  return kERROR;
163  }
164 
165  // Get rec. hits [needed only for tests!!!]
166  fRecHits = (TClonesArray *)ioman->GetObject("LMDHitsMerged");
167  if (!fRecHits) {
168  std::cout << "-W- PndLmdBPtestTask::Init: "
169  << "No LMDHitsMerged [needed only for tests!!!]"
170  << " array!" << std::endl;
171  return kERROR;
172  }
173 
174  fTrackParGeane = new TClonesArray("FairTrackParH");
175  ioman->Register("GeaneTrackPar", "Geane", fTrackParGeane, kTRUE);
176 
177  fTrackParIni = new TClonesArray("FairTrackParH");
178  ioman->Register("GeaneTrackIni", "Geane", fTrackParIni, kTRUE);
179 
180  fTrackParFinal = new TClonesArray("FairTrackParH");
181  ioman->Register("GeaneTrackFinal", "Geane", fTrackParFinal, kTRUE);
182 
183  fDetName = new TClonesArray("TObjString");
184  ioman->Register("DetName", "Geane", fDetName, kTRUE);
185 
186  fPro = new FairGeanePro();
187  fGeoH = PndGeoHandling::Instance();
188  // FairRun* fRun = FairRun::Instance(); //[R.K. 01/2017] unused variable?
189  // FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); //[R.K. 01/2017] unused
190  // variable
191 
192  pndField = FairRunAna::Instance()->GetField();
193 
194  return kSUCCESS;
195 }
196 // -------------------------------------------------------------------------
197 void PndLmdBPtestTask::SetParContainers() {
198  // Get Base Container
200  // FairRuntimeDb* rtdb=ana->GetRuntimeDb();
201 }
202 
203 // ----- Public method Exec --------------------------------------------
204 void PndLmdBPtestTask::Exec(Option_t *) {
205  // if(fVerbose>5){
206  // if((fTracks->GetEntries())!=(fMCTracks->GetEntries()))
207  // return;
208  // }
209  // cout<<"PndLmdBPtestTask::Exec starts!"<<endl;
210  std::map<int, std::vector<int> > mcHitMap; // Track -> MCHits
211  fTrackParGeane->Delete();
212  fTrackParIni->Delete();
213  fTrackParFinal->Delete();
214  fDetName->Delete();
215 
216  mcHitMap = AssignHitsToTracks();
217 
218  if (fVerbose > 2) {
219  cout << " ---- Info: " << fEventNr << endl;
220  }
221  fEventNr++;
222  // Charge & mass of particle
223  Int_t PDGCode = -2212; // antiproton
224  // Int_t PDGCode = fPDGid;
225  TDatabasePDG *fdbPDG = TDatabasePDG::Instance();
226  TParticlePDG *fParticle = fdbPDG->GetParticle(PDGCode);
227  Double_t fCharge = fParticle->Charge() / 3.;
228  // cout<<"fCharge = "<<fCharge<<endl;
229 
230  // go through all tracks
231  int glI = fTracks->GetEntriesFast();
232  // cout<<"glI = "<<glI<<endl;
233  Int_t counterGeaneTrk = 0;
234 
235  for (Int_t i = 0; i < glI; i++) {
236  // cout<<"PndLmdBPtestTask::Exec for track#"<<i<<endl;
237  TVector3 StartPos, StartPosErr, StartMom, StartMomErr, StartO, StartU,
238  StartV;
239  // int p = 0; //[R.K. 01/2017] unused variable
240 
241  PndTrack *recTrack = (PndTrack *)(fTracks->At(i));
242  FairTrackParP fFittedTrkP = recTrack->GetParamFirst();
243 
244  TVector3 PosRecLMD(fFittedTrkP.GetX(), fFittedTrkP.GetY(),
245  fFittedTrkP.GetZ());
246  if (fFittedTrkP.GetZ() > 1130)
247  continue; // TEST: skip trks from 2nd plane. TODO: check are they fine???
248  TVector3 MomRecLMD(fFittedTrkP.GetPx(), fFittedTrkP.GetPy(),
249  fFittedTrkP.GetPz());
250  MomRecLMD *= fPbeam / MomRecLMD.Mag(); // external assumption about mom
251  // magnitude //<- commented as TEST
252  fFittedTrkP.SetPx(MomRecLMD.X());
253  fFittedTrkP.SetPy(MomRecLMD.Y());
254  fFittedTrkP.SetPz(MomRecLMD.Z());
255  double covMARS[6][6];
256  fFittedTrkP.GetMARSCov(covMARS);
257  TVector3 errMomRecLMD(sqrt(covMARS[0][0]), sqrt(covMARS[1][1]),
258  sqrt(covMARS[2][2]));
259  TVector3 errPosRecLMD(sqrt(covMARS[3][3]), sqrt(covMARS[4][4]),
260  sqrt(covMARS[5][5]));
261 
262  StartPos = PosRecLMD;
263  StartPosErr = errPosRecLMD;
264  StartMom = MomRecLMD;
265  StartMomErr = errMomRecLMD;
266 
267  // Make true2 track from MC info on 1st plane of LMD-----------
268  FairTrackParP *fStartMCLMD = new FairTrackParP();
269  if (fVerbose > 5 && (fTracks->GetEntries() == fMCTracks->GetEntries())) {
270  int candID = recTrack->GetRefIndex();
271  PndTrackCand *trkcand = (PndTrackCand *)fRecCandTracks->At(candID);
272  PndTrackCandHit candhit =
273  (PndTrackCandHit)(trkcand->GetSortedHit(0)); // 1st hit info
274  Int_t hitID = candhit.GetHitId();
275  PndSdsMergedHit *myHit = (PndSdsMergedHit *)(fRecHits->At(hitID));
276  int mcrefbot = myHit->GetSecondMCHit();
277  int mcreftop = myHit->GetRefIndex();
278  PndSdsMCPoint *MCPointHit;
279  if (mcreftop >= 0) {
280  MCPointHit = (PndSdsMCPoint *)(fMCHits->At(mcreftop));
281  } else {
282  MCPointHit = (PndSdsMCPoint *)(fMCHits->At(mcrefbot));
283  }
284 
285  TVector3 PosMClmd = MCPointHit->GetPosition();
286  double pxTrue = MCPointHit->GetPx();
287  double pyTrue = MCPointHit->GetPy();
288  double pzTrue = MCPointHit->GetPz();
289  TVector3 MomMClmd(pxTrue, pyTrue, pzTrue);
290  TVector3 dirMClmd = MomMClmd;
291  dirMClmd *= 1. / (MomMClmd.Mag());
292  MomMClmd *= fPbeam / (MomMClmd.Mag()); //<- commented as TEST
293  double xMClmdNew = PosMClmd.X() - (dirMClmd.X() * 0.010);
294  double yMClmdNew = PosMClmd.Y() - (dirMClmd.Y() * 0.010);
295  double zMClmdNew = PosMClmd.Z() - (dirMClmd.Z() * 0.010);
296  PosMClmd.SetXYZ(xMClmdNew, yMClmdNew, zMClmdNew); // shift 100 mkm to
297  // adoid adding MC
298  // errors during back
299  // propagation;
300  TVector3 ocMC(0, 0, 0); // define plane perpendicular to z-axis in IP
301  TVector3 djMC(1., 0., 0.);
302  TVector3 dkMC(0., 1., 0.);
303  TVector3 PosMClmderr(0., 0., 0.);
304  TVector3 MomMClmderr(0., 0., 0.);
305  fStartMCLMD = new FairTrackParP(PosMClmd, MomMClmd, PosMClmderr,
306  MomMClmderr, fCharge, ocMC, djMC, dkMC);
307  }
308  //-------------------------------------------------------------------------------
309 
312  // Comment: seems back propagation in one step (for 11 m) is too much
313  // //try smoothing it by small steps, which follow mag.field
314  // const int nstep=7;
315  // double zbend[nstep]={661, 660.5, 660., 659, 319, 316, 220};//entarance
316  // and exit mag.field
317  // //TEST for backward and forward propagation: more steps!
318  // const int nstep=50;
319  // double zbend[nstep];
320  // const double z0=fFittedTrkP.GetZ();
321  // const double z1=1;
322  // const double zstep=(z0-z1)/nstep;
323  // for(int js=0;js<nstep;js++){
324  // zbend[js]=z0-zstep*js;
325  // }
326 
327  const int nstep0 = 10;
328  // double zbend0[nstep0]={fFittedTrkP.GetZ(), 661, 660.5, 660., 659,
329  // 319, 316, 220,10};//entarance and exit mag.field
330  double zbend0[nstep0] = {fFittedTrkP.GetZ(),
331  660,
332  602,
333  450,
334  342,
335  283,
336  248,
337  180,
338  100,
339  1}; // entarance and exit mag.field
340  // TEST for backward and forward propagation: more steps!
341  const int nintstep = 100;
342  const int nstep = nstep0 * nintstep;
343  double zbend[nstep];
344  for (int is = 1; is <= nstep0; is++) {
345  const double z0 = zbend0[is - 1];
346  const double z1 = zbend0[is];
347  const double zstep = (z0 - z1) / nintstep;
348  // cout<<"is = "<<is<<": z0 = "<<z0<<" z1 = "<<z1<<" zstep =
349  //"<<zstep<<endl;
350  for (int js = 0; js < nintstep; js++) {
351  int curint = (is - 1) * nintstep + js;
352  zbend[curint] = z0 - zstep * js;
353  // cout<<"zbend["<<curint<<"]="<<zbend[curint]<<endl;
354  }
355  }
356 
357  vector<double> vxmc(nstep), vymc(nstep), vxmc_err(nstep), vymc_err(nstep),
358  vzmc(nstep), vthetamc(nstep), vphimc(nstep), vpmc(nstep), vvmc(nstep),
359  vwmc(nstep), vtvmc(nstep), vtwmc(nstep), vvmc_err(nstep),
360  vwmc_err(nstep), vtvmc_err(nstep), vtwmc_err(nstep);
361  FairTrackParP *fStartMC = new FairTrackParP();
362  if (fVerbose > 2) {
363  cout << "------------------------------------------" << endl;
364  cout << "StartPos:" << endl;
365  StartPos.Print();
366  cout << "StartPosErr:" << endl;
367  StartPosErr.Print();
368 
369  cout << "StartMom: " << StartMom.Mag() << endl;
370  StartMom.Print();
371  cout << "StartMomErr: " << StartMomErr.Mag() << endl;
372  StartMomErr.Print();
373  cout << "" << endl;
374 
375  if (fVerbose > 5 && (fTracks->GetEntries() ==
376  fMCTracks->GetEntries())) { // check results for MC
377  // track propagation to
378  // LMD plane
379  FairTrackParP *fStartPst = new FairTrackParP(fFittedTrkP);
380  TVector3 istLMD(fStartPst->GetIVer());
381  TVector3 jstLMD(fStartPst->GetJVer());
382  TVector3 kstLMD(fStartPst->GetKVer());
383  // cout<<"i,j,k of LMD trk:"<<endl;
384  // istLMD+=StartPos;
385  istLMD = StartPos;
386  // istLMD.Print();
387  // jstLMD.Print();
388  // kstLMD.Print();
389  TVector3 ocMC(0, 0, 0); // define plane perpendicular to z-axis in IP
390  TVector3 djMC(1., 0., 0.);
391  TVector3 dkMC(0., 1., 0.);
392  djMC.SetMag(1);
393  dkMC.SetMag(1);
394 
395  PndMCTrack *mctrk = (PndMCTrack *)(fMCTracks->At(i));
396  TVector3 MomMC = mctrk->GetMomentum();
397  TVector3 PosMC = mctrk->GetStartVertex();
398  TVector3 MomMCerr(0., 0., 0.);
399  TVector3 PosMCerr(0., 0., 0.);
400  // if(PosMC.Z()>0){
401  if (PosMC.Z() > 1) {
402  cout << "!!! Achtung: " << PosMC.Z() << endl;
403  break;
404  }
405 
406  fStartMC = new FairTrackParP(PosMC, MomMC, PosMCerr, MomMCerr, fCharge,
407  ocMC, djMC, dkMC);
408  for (int sj = nstep - 1; sj > -1; sj--) {
409  // cout<<"MC forward: to z="<<zbend[sj]<<endl;
410  bool isPropMC;
411  FairTrackParP *fResMC = PropToPlane(fStartMC, zbend[sj], +1,
412  isPropMC); // forward propagation
413  if (isPropMC) {
414  delete fStartMC;
415  fStartMC = fResMC;
416  vxmc[sj] = fResMC->GetX();
417  vxmc_err[sj] = fResMC->GetDX();
418  vymc_err[sj] = fResMC->GetDY();
419  vymc[sj] = fResMC->GetY();
420  vzmc[sj] = fResMC->GetZ();
421  TVector3 finDirMC(fResMC->GetPx(), fResMC->GetPy(),
422  fResMC->GetPz());
423  vpmc[sj] = finDirMC.Mag();
424  vthetamc[sj] = finDirMC.Theta();
425  vphimc[sj] = finDirMC.Phi();
426 
427  vvmc[sj] = fResMC->GetV();
428  vwmc[sj] = fResMC->GetW();
429  vtvmc[sj] = fResMC->GetTV();
430  vtwmc[sj] = fResMC->GetTW();
431  vvmc_err[sj] = fResMC->GetDV();
432  vwmc_err[sj] = fResMC->GetDW();
433  vtvmc_err[sj] = fResMC->GetDTV();
434  vtwmc_err[sj] = fResMC->GetDTW();
435  // cout<<"Next step;)"<<endl;
436  } else
437  break;
438  }
439  }
440  }
441  if (abs(fStartMC->GetZ() - fFittedTrkP.GetZ()) > 10 && fVerbose > 5)
442  break; // MC particle wasn't propagated to LMD plane ????
443  TClonesArray &clref1 = *fTrackParIni;
444  Int_t size1 = clref1.GetEntriesFast();
445 
446  FairTrackParP *fStartPst = new FairTrackParP(fFittedTrkP); // REC
447  FairTrackParP *fStartPstMCLMD =
448  new FairTrackParP(*fStartMCLMD); // MC in LMD, should be treated as REC
449  for (int js = 1; js < nstep; js++) {
450  // cout<<"step#"<<js<<endl;
451  bool isProp;
452  // cout<<"REC backward "<<endl;
453  FairTrackParP *fResPst =
454  PropToPlane(fStartPst, zbend[js], -1, isProp); // back propagation
455  if (isProp) {
456  delete fStartPst;
457  fStartPst = fResPst;
458  if (fVerbose > 5) {
459  bool isPropMClmd;
460  FairTrackParP *fResPstMCLMD = PropToPlane(
461  fStartPstMCLMD, zbend[js], -1, isPropMClmd); // back propagation
462  if (isPropMClmd) {
463  TVector3 finPosMCLMD(fResPstMCLMD->GetX(), fResPstMCLMD->GetY(),
464  fResPstMCLMD->GetZ());
465  TVector3 finDirMCLMD(fResPstMCLMD->GetPx(), fResPstMCLMD->GetPy(),
466  fResPstMCLMD->GetPz());
467 
468  TVector3 finPos(fResPst->GetX(), fResPst->GetY(), fResPst->GetZ());
469  TVector3 finDir(fResPst->GetPx(), fResPst->GetPy(),
470  fResPst->GetPz());
471 
472  // double pnt[3]={finPosMC.X(),finPosMC.Y(),finPosMC.Z()};
473  // //Position where to get field strength
474  double pnt[3] = {vxmc[js], vymc[js], vzmc[js]};
475  double Bf[3]; // result goes here
476  // retrieve the field from the framework
477  pndField->Field(pnt, Bf); //[kGs]
478  fbx = Bf[0];
479  fby = Bf[1];
480  fbz = Bf[2];
481  fxrec = finPos.X();
482  fyrec = finPos.Y();
483  fxrec_err = fResPst->GetDX();
484  fyrec_err = fResPst->GetDY();
485  fzrec = finPos.Z();
486  fprec = finDir.Mag();
487  fthetarec = finDir.Theta();
488  fphirec = finDir.Phi();
489  fxmc = vxmc[js];
490  fymc = vymc[js];
491  fzmc = vzmc[js];
492  fxmc_err = vxmc_err[js];
493  fymc_err = vymc_err[js];
494  fpmc = vpmc[js];
495  fthetamc = vthetamc[js];
496  fphimc = vphimc[js];
497 
498  fxmclmd = finPosMCLMD.X();
499  fymclmd = finPosMCLMD.Y();
500  fxmclmd_err = fResPstMCLMD->GetDX();
501  fymclmd_err = fResPstMCLMD->GetDY();
502  fzmclmd = finPosMCLMD.Z();
503  fpmclmd = finDirMCLMD.Mag();
504  fthetamclmd = finDirMCLMD.Theta();
505  fphimclmd = finDirMCLMD.Phi();
506 
507  fvmc = vvmc[js];
508  fwmc = vwmc[js];
509  ftvmc = vtvmc[js];
510  ftwmc = vtwmc[js];
511  fvmc_err = vvmc_err[js];
512  fwmc_err = vwmc_err[js];
513  ftvmc_err = vtvmc_err[js];
514  ftwmc_err = vtwmc_err[js];
515 
516  fvrec = fResPst->GetV();
517  fwrec = fResPst->GetW();
518  ftvrec = fResPst->GetTV();
519  ftwrec = fResPst->GetTW();
520  fvrec_err = fResPst->GetDV();
521  fwrec_err = fResPst->GetDW();
522  ftvrec_err = fResPst->GetDTV();
523  ftwrec_err = fResPst->GetDTW();
524 
525  fvmclmd = fResPstMCLMD->GetV();
526  fwmclmd = fResPstMCLMD->GetW();
527  ftvmclmd = fResPstMCLMD->GetTV();
528  ftwmclmd = fResPstMCLMD->GetTW();
529  fvmclmd_err = fResPstMCLMD->GetDV();
530  fwmclmd_err = fResPstMCLMD->GetDW();
531  ftvmclmd_err = fResPstMCLMD->GetDTV();
532  ftwmclmd_err = fResPstMCLMD->GetDTW();
533 
534  tprop->Fill();
535  delete fStartPstMCLMD;
536  fStartPstMCLMD = fResPstMCLMD;
537  } else
538  break;
539  }
540  } else
541  break;
542  }
543 
546  int ierr = 0;
547  FairTrackParH *fStart = new (clref1[size1]) FairTrackParH(fStartPst, ierr);
548  TClonesArray &clref = *fTrackParGeane;
549  Int_t size = clref.GetEntriesFast();
550  FairTrackParH *fRes = new (clref[size]) FairTrackParH();
551  fPro->SetPoint(vtx);
552  fPro->PropagateToPCA(1, -1); // back-propagate to point
553  Bool_t isProp = fPro->Propagate(fStart, fRes, PDGCode);
554 
555  // ///----------------------------------------------------------------------
556  delete fStartPst;
557 
559 
560  //------- TEST of calculation errors -------
561  TVector3 gPos(fRes->GetX(), fRes->GetY(), fRes->GetZ());
562  TVector3 gMom(fRes->GetPx(), fRes->GetPy(), fRes->GetPz());
563 
564  TVector3 gErrPos(fRes->GetDX(), fRes->GetDY(), fRes->GetDZ());
565  TVector3 gErrMom(fRes->GetDPx(), fRes->GetDPy(), fRes->GetDPz());
566  // cout<<" "<<endl;
567  if (fVerbose > 2) {
568  // cout<<"================= %%%% ===================="<<endl;
569 
570  cout << "gPos:" << endl;
571  gPos.Print();
572  // cout<<"difference between final and initial point =
573  // "<<(-StartPos.Mag()+gPos.Mag())<<endl;
574  cout << "gErrPos:" << endl;
575  gErrPos.Print();
576  // cout<<"difference between final and initial point error =
577  // "<<(-StartPosErr.Mag()+gErrPos.Mag())<<endl;
578 
579  cout << "gMom: " << gMom.Mag() << endl;
580  gMom.Print();
581  // cout<<"difference between initial and final momentum =
582  // "<<(StartMom.Mag()-gMom.Mag())<<endl;
583  cout << "gErrMom: " << gErrMom.Mag() << endl;
584  // cout<<"difference between initial and final momentum error=
585  // "<<(StartMomErr.Mag()-gErrMom.Mag())<<endl;
586  gErrMom.Print();
587  if (fVerbose > 5 && (fTracks->GetEntries() == fMCTracks->GetEntries())) {
588  PndMCTrack *mctrk = (PndMCTrack *)(fMCTracks->At(i));
589  int motherID = mctrk->GetMotherID();
590  if (motherID >= 0) break;
591  TVector3 MomMC = mctrk->GetMomentum();
592  TVector3 PosMC = mctrk->GetStartVertex();
593  // cout<<"mcPos:"<<endl;
594  // PosMC.Print();
595  // cout<<"mcMom: "<<MomMC.Mag()<<endl;
596  // MomMC.Print();
597 
598  double pnt[3] = {PosMC.X(), PosMC.Y(),
599  PosMC.Z()}; // Position where to get field strength
600  double Bf[3]; // result goes here
601  // retrieve the field from the framework
602 
603  pndField->Field(pnt, Bf); //[kGs]
604  // cout<<"^^^^ Bf = ("<<Bf[0]<<", "<<Bf[1]<<", "<<Bf[2]<<")"<<endl;
605  fbx = Bf[0];
606  fby = Bf[1];
607  fbz = Bf[2];
608  fxrec = gPos.X();
609  fyrec = gPos.Y();
610  fzrec = gPos.Z();
611  fprec = gMom.Mag();
612  fthetarec = gMom.Theta();
613  fphirec = gMom.Phi();
614  fxmc = PosMC.X();
615  fymc = PosMC.Y();
616  fzmc = PosMC.Z();
617  fpmc = MomMC.Mag();
618  fthetamc = MomMC.Theta();
619  fphimc = MomMC.Phi();
620 
621  FairTrackParH *fStartMClmd = new FairTrackParH(fStartPstMCLMD, ierr);
622  fPro->SetPoint(vtx);
623  fPro->PropagateToPCA(1, -1); // back-propagate to point
624  FairTrackParH *fResMCLMD = new FairTrackParH();
625  fPro->Propagate(
626  fStartMClmd, fResMCLMD,
627  PDGCode); // Bool_t isPropMC = //[R.K.03/2017] unused variable
628  TVector3 gPosMC(fResMCLMD->GetX(), fResMCLMD->GetY(),
629  fResMCLMD->GetZ());
630  TVector3 gMomMC(fResMCLMD->GetPx(), fResMCLMD->GetPy(),
631  fResMCLMD->GetPz());
632  fxmclmd = gPosMC.X();
633  fymclmd = gPosMC.Y();
634  fzmclmd = gPosMC.Z();
635  fpmclmd = gMomMC.Mag();
636  fthetamclmd = gMomMC.Theta();
637  fphimclmd = gMomMC.Phi();
638  tprop->Fill();
639  delete fStartPstMCLMD;
640  }
641  cout << "================= %%%% ====================" << endl;
642  }
643  //------------------------------------------
644 
645  if (isProp == kTRUE) {
646  new ((*fTrackParFinal)[counterGeaneTrk])
647  FairTrackParH(*(fRes)); // save Track
648  // new((*fTrackParFinal)[counterGeaneTrk]) FairTrackParP(*(fRes));
650  counterGeaneTrk++;
651  if (fVerbose > 2) cout << "***** isProp TRUE *****" << endl;
652  } else {
653  if (fVerbose > 2) {
654  cout << "!!! Back-propagation with GEANE didn't return result !!!"
655  << endl;
656  cout << "StartPos:" << endl;
657  StartPos.Print();
658  cout << "StartPosErr:" << endl;
659  StartPosErr.Print();
660 
661  cout << "StartMom: " << StartMom.Mag() << endl;
662  StartMom.Print();
663  cout << "StartMomErr: " << StartMomErr.Mag() << endl;
664  StartMomErr.Print();
665  }
666  new ((*fTrackParFinal)[counterGeaneTrk]) FairTrackParH(); // save NULL
667  // new((*fTrackParFinal)[counterGeaneTrk]) FairTrackParP(); //save
668  //NULL
669  counterGeaneTrk++;
670  }
671  }
672  // fMCTracks->Delete();
673  // fMCHits->Delete();
674  if (fVerbose > 2) cout << "PndLmdBPtestTask::Exec END!" << endl;
675 }
676 
677 void PndLmdBPtestTask::Finish() {
678  if (fVerbose > 5) {
679  TTree *nout1 = tprop->CloneTree();
680  nout1->Write();
681  }
682 }
683 
684 FairTrackParP *PndLmdBPtestTask::PropToPlane(FairTrackParP *fStartPst,
685  double zpos, int dir,
686  bool &isProp) {
687  TVector3 stStartPos(fStartPst->GetX(), fStartPst->GetY(), fStartPst->GetZ());
688  if (zpos > 1e3) { // external assumption about mom magnitude [use it while in
689  // BOX and we know mom should be diff]
690  TVector3 MomStartPos(fStartPst->GetPx(), fStartPst->GetPy(),
691  fStartPst->GetPz());
692  MomStartPos *= fPbeam / MomStartPos.Mag(); //<- commented as TEST
693  fStartPst->SetPx(MomStartPos.X());
694  fStartPst->SetPy(MomStartPos.Y());
695  fStartPst->SetPz(MomStartPos.Z()); // correct mom.magnitude
696  }
697  // propagate plane-to-plane
698  TVector3 ist(fStartPst->GetIVer());
699  TVector3 jst(fStartPst->GetJVer());
700  TVector3 kst(fStartPst->GetKVer());
701  TVector3 oc(0, 0, zpos);
702  TVector3 dj(1., 0., 0.);
703  TVector3 dk(0., 1., 0.);
704  dj.SetMag(1);
705  dk.SetMag(1);
706  fPro->PropagateFromPlane(jst, kst); // 1st detector plane
707  fPro->PropagateToPlane(oc, dj, dk); // virtual plane at fixed z
708  if (dir < 0) fPro->setBackProp();
709  FairTrackParP *fResPst = new FairTrackParP();
710  isProp = fPro->Propagate(fStartPst, fResPst, fPDGid);
711  return fResPst;
712 }
713 
714 std::map<int, std::vector<int> > PndLmdBPtestTask::AssignHitsToTracks() {
715  std::map<int, std::vector<int> > result;
716  for (int i = 0; i < fMCHits->GetEntriesFast(); i++) { // get all MC Hits
717  PndSdsMCPoint *myPoint =
718  (PndSdsMCPoint *)(fMCHits->At(i)); // sort MCHits with Tracks
719  // PndMCTrack* myTrack =
720  // (PndMCTrack*)(fMCTracks->At(myPoint->GetTrackID()));
721  result[myPoint->GetTrackID()].push_back(i);
722  }
723  return result;
724 }
725 
726 ClassImp(PndLmdBPtestTask);
PndLmdBPtestTask::fGeoH
PndGeoHandling * fGeoH
Definition: PndLmdBPtestTask.h:64
z0
Double_t z0
Definition: checkhelixhit.C:62
fVerbose
int fVerbose
Definition: poormantracks.C:24
PndLmdBPtestTask::fTrackParFinal
TClonesArray * fTrackParFinal
Definition: PndLmdBPtestTask.h:60
i
Int_t i
Definition: run_full.C:25
PndLmdBPtestTask::Finish
virtual void Finish()
Definition: PndLmdBPtestTask.cxx:677
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndTrackCand::GetSortedHit
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
pnt
TClonesArray * pnt
Definition: checkdedx_helixhit.C:12
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndLmdBPtestTask::AssignHitsToTracks
std::map< int, std::vector< int > > AssignHitsToTracks()
Definition: PndLmdBPtestTask.cxx:714
PndLmdBPtestTask::fDetName
TClonesArray * fDetName
Definition: PndLmdBPtestTask.h:61
PndLmdBPtestTask::Exec
virtual void Exec(Option_t *opt)
Definition: PndLmdBPtestTask.cxx:204
PndLmdBPtestTask::fTracks
TClonesArray * fTracks
Definition: PndLmdBPtestTask.h:78
PndLmdBPtestTask::Init
virtual InitStatus Init()
Definition: PndLmdBPtestTask.cxx:121
fParticle
Int_t * fParticle
Definition: run_full.C:24
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
zpos
const Double_t zpos
Definition: createRootFscGeometryFile.C:11
PndLmdBPtestTask::fMCHits
TClonesArray * fMCHits
Definition: PndLmdBPtestTask.h:75
PndLmdBPtestTask::SetParContainers
virtual void SetParContainers()
Definition: PndLmdBPtestTask.cxx:197
is
static int is
Definition: ranlxd.cxx:374
PndGeoHandling::Instance
static PndGeoHandling * Instance()
Definition: PndGeoHandling.cxx:36
PndLmdBPtestTask::pndField
FairField * pndField
Definition: PndLmdBPtestTask.h:70
PndLmdBPtestTask::fPro
FairGeanePro * fPro
Definition: PndLmdBPtestTask.h:63
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
PndLmdBPtestTask::fMCTracks
TClonesArray * fMCTracks
Definition: PndLmdBPtestTask.h:76
PndSdsMergedHit::GetSecondMCHit
Int_t GetSecondMCHit() const
Definition: PndSdsMergedHit.h:78
PndLmdBPtestTask::fRecHits
TClonesArray * fRecHits
Definition: PndLmdBPtestTask.h:81
ClassImp
ClassImp(PndLmdBPtestTask)
PndLmdBPtestTask::fTrackParIni
TClonesArray * fTrackParIni
Definition: PndLmdBPtestTask.h:59
PndLmdBPtestTask::fTrackParGeane
TClonesArray * fTrackParGeane
Definition: PndLmdBPtestTask.h:58
PndTrack::GetRefIndex
Int_t GetRefIndex() const
Definition: PndTrack.h:36
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
PndTrackCandHit::GetHitId
Int_t GetHitId() const
Definition: PndTrackCandHit.h:55
PndLmdBPtestTask::fRecCandTracks
TClonesArray * fRecCandTracks
Definition: PndLmdBPtestTask.h:82
PndTrack::GetParamFirst
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
PndLmdBPtestTask::PropToPlane
FairTrackParP * PropToPlane(FairTrackParP *fStartPst, double zpos, int dir, bool &isProp)
Definition: PndLmdBPtestTask.cxx:684
Generated on Sun Aug 18 2019 03:01:35 for FairRoot/PandaRoot by   doxygen 1.8.8