FairRoot/PandaRoot
plotTrackCands.C
Go to the documentation of this file.
1 // Macros used to visualize the event by plotting
2 // B field
3 // MC points
4 // hits
5 // found track candidates
6 
7 
8 #ifndef HOUGH_TEST
9 #define HOUGH_TEST
10 
11 
12 #include <math.h>
13 
14 
15 #include "TString.h"
16 #include "TROOT.h"
17 #include "TSystem.h"
18 #include "TH1.h"
19 #include "TH2.h"
20 #include "TCanvas.h"
21 #include "TLine.h"
22 #include "TLorentzVector.h"
23 #include "TStyle.h"
24 #include "TRandom.h"
25 #include "TStopwatch.h"
26 #include "TGraph.h"
27 #include "TSpectrum2.h"
28 
29 
30 #include <stdio.h>
31 #include <iostream>
32 #include <iomanip>
33 using namespace std;
34 
35 // Panda Libs
36 //#include "RhoBase/TCandidate.h"
37 //#include "RhoBase/TCandList.h"
38 //#include "RhoBase/TFitParams.h"
39 //#include "RhoBase/VAbsMicroCandidate.h"
40 //#include "RhoSelector/TPidSelector.h"
41 
42 #include "FairLogger.h"
43 //#include "FairRunAna.h"
44 #include "FairRuntimeDb.h"
45 #include "FairParRootFileIo.h"
46 
47 #include "PndFtsHit.h"
48 #include "FairMCPoint.h"
49 #include "FairHit.h"
50 
51 
52 //#include "AnalysisTools/PndEventReader.h"
53 //#include "AnalysisTools/PndMcTruthMatch.h"
54 #include "AnalysisTools/PndAnalysis.h"
55 
56 #include "AnalysisTools/Fitter/Rho4CFitter.h"
57 #include "AnalysisTools/Fitter/RhoKinVtxFitter.h"
58 #include "AnalysisTools/Fitter/RhoKinFitter.h"
59 
60 #include "AnalysisTools/PndAnaPidCombiner.h"
61 #include "AnalysisTools/PndAnaPidSelector.h"
62 #include "FairRunSim.h"
63 #include "PndMultiField.h"
64 #include "PndMCTrack.h"
65 #include "PndFtsHoughTrackCand.h"
66 
67 
68 
69 
70 
71 
72 
73 
74 
75 int PlotMCTracks(int iEvent, TString branchName, TString tcaName, TGraph *mcPoints, TFile* fSim = 0)
76 {
77  // if no file is specified the newest open file is used
78 
79  // Uncomment if you want to test the routine without calling it in a macro (by copy pasting into the interpreter)
80  // TString branchName="FTSPoint"; TString tcaName="PndFtsPoint";
81  // TString branchName="STTPoint"; TString tcaName="PndSttPoint";
82 
83  int verbose = 0;
84 
85  // use last file if none is specified
86  if (0==fSim) {
87  fSim = gFile;
88  }
89 
90  TTree* tSim = (TTree*) (gFile->Get("pndsim"));
91  //tSim->StartViewer();
92  TClonesArray* Points = new TClonesArray(tcaName);
93  tSim->SetBranchAddress(branchName, &Points);
94  TBranch* br = tSim->FindBranch(branchName);
95 
96  if (br == 0 || br->GetEntries() == 0 || iEvent >= br->GetEntries())
97  return;
98 
99 
100  tSim->GetEntry(iEvent); // Look at event iEvent
101  const Int_t nPoints = Points->GetEntriesFast();
102  // store coordinates of points for plotting in TGraph
103  Double_t X[nPoints];
104  Double_t Y[nPoints]; // not used for plotting
105  Double_t Z[nPoints];
106 
107 
108  for (int iPoint = 0; iPoint < nPoints; iPoint++)
109  {
110  TObject* obj = Points->At(iPoint);
111  if (iPoint < 1 && verbose > 5)
112  {
113  std::cout << obj->ClassName() << std::endl;
114  std::cout << "obj->InheritsFrom(FairMCPoint) = "
115  << obj->InheritsFrom("FairMCPoint") << std::endl;
116  std::cout << "obj->InheritsFrom(FairHit) = "
117  << obj->InheritsFrom("FairHit") << std::endl;
118  }
119  TVector3 point;
120  if (obj->InheritsFrom("FairMCPoint"))
121  {
122  FairMCPoint* myPoint = (FairMCPoint*) obj;
123  myPoint->Position(point);
124  }
125  else if (obj->InheritsFrom("FairHit"))
126  {
127  FairHit* myHit = (FairHit*) obj;
128  myHit->Position(point);
129  }
130  else
131  {
132  std::cout
133  << "Error, the hit or point does neither inherit from FairHit nor from FairMCPoint \n";
134  }
135 
136  // Hit coordinates
137  X[iPoint] = point.X();
138  Y[iPoint] = point.Y();
139  Z[iPoint] = point.Z();
140 
141  //mcPoints->Fill(Z, X);
142 
143  if (verbose > 2)
144  {
145  std::cout << "Point " << iPoint << " : ";
146  //myPoint->Print();
147  std::cout << "Z=" << Z[iPoint] << "cm X=" << X[iPoint] << "cm Y=" << Y[iPoint] << "cm"
148  << std::endl;
149  }
150 
151  }
152 
153  if (0<nPoints)
154  {
155  mcPoints->DrawGraph(nPoints, Z, X, "P,SAME");
156  // DrawGraph(Int_t n, const Double_t *x=0, const Double_t *y=0, Option_t *option="");
157  }
158 
159 
160  tSim->ResetBranchAddresses();
161  Points->Delete(); // maybe Points->Clear(); would be enough...
162  return 0;
163 }
164 
165 
166 int PlotMCTracksPrintBField(PndMultiField *fField, int iEvent,
167  TString branchName, TString tcaName, TH2D* mcPoints)
168 {
169  // File that you want to access needs to be opened first (and not other file after it), because the file is accessed via gFile
170  // TODO I could also pass in the handle instead of using gFile
171 
172  // Uncomment if you want to test the routine without calling it in a macro (by copy pasting into the interpreter)
173  // TString branchName="FTSPoint"; TString tcaName="PndFtsPoint";
174  // TString branchName="STTPoint"; TString tcaName="PndSttPoint";
175 
176  int verbose = 6;
177 
178  TTree* tSim = (TTree*) (gFile->Get("pndsim"));
179  //tSim->StartViewer();
180  if (verbose > 5)
181  {
182  std::cout << "Got the tree" << std::endl;
183  }
184 
185  TClonesArray* Points = new TClonesArray(tcaName);
186  if (verbose > 5)
187  {
188  std::cout << "Created new TCA" << std::endl;
189  }
190  tSim->SetBranchAddress(branchName, &Points);
191  if (verbose > 5)
192  {
193  std::cout << "Set the branchName" << std::endl;
194  }
195  TBranch* br = tSim->FindBranch(branchName);
196  if (verbose > 5)
197  {
198  std::cout << "Found the branch" << std::endl;
199  }
200  if (br == 0 || br->GetEntries() == 0 || iEvent >= br->GetEntries())
201  return;
202  if (verbose > 3)
203  {
204  std::cout << "Passed the test" << std::endl;
205  }
206  tSim->GetEntry(iEvent); // Look at event iEvent
207  if (verbose > 2)
208  {
209  std::cout << "Number of points: " << Points->GetEntriesFast()
210  << std::endl;
211  }
212  for (int iPoint = 0; iPoint < Points->GetEntriesFast(); iPoint++)
213  {
214  TObject* obj = Points->At(iPoint);
215  if (iPoint < 1 && verbose > 5)
216  {
217  std::cout << obj->ClassName() << std::endl;
218  std::cout << "obj->InheritsFrom(FairMCPoint) = "
219  << obj->InheritsFrom("FairMCPoint") << std::endl;
220  std::cout << "obj->InheritsFrom(FairHit) = "
221  << obj->InheritsFrom("FairHit") << std::endl;
222  }
223  TVector3 point;
224  if (obj->InheritsFrom("FairMCPoint"))
225  {
226  FairMCPoint* myPoint = (FairMCPoint*) obj;
227  myPoint->Position(point);
228  }
229  else if (obj->InheritsFrom("FairHit"))
230  {
231  FairHit* myHit = (FairHit*) obj;
232  myHit->Position(point);
233  }
234  else
235  {
236  std::cout
237  << "Error, the hit or point does neither inherit from FairHit nor from FairMCPoint \n";
238  }
239 
240  // Hit coordinates
241  Double_t X = point.X();
242  Double_t Y = point.Y();
243  Double_t Z = point.Z();
244 
245  mcPoints->Fill(Z, X);
246 
247  // Get the magnetic field at the point location
248  // This draws the magnetic field in pointer fField into histogram Bfield
249 
250  Double_t po[3], BB[3];
251  Double_t Btot = 0;
252 
253  po[0] = X;
254  po[1] = Y;
255  po[2] = Z;
256 
257  fField->GetFieldValue(po, BB); //return value in KG (G3)
258 
259  if (verbose > 2)
260  {
261  std::cout << "Point " << iPoint << " : ";
262  //myPoint->Print();
263  std::cout << "Z=" << Z << "cm X=" << X << "cm Y=" << Y << "cm BX="
264  << BB[0] / 10. << "T BY=" << BB[1] / 10. << "T BZ="
265  << BB[2] / 10. << "T" << std::endl;
266  }
267 
268  }
269 
270  mcPoints->SetMarkerColor(2);
271  mcPoints->SetMarkerStyle(29);
272  mcPoints->SetMarkerSize(1);
273 
274  tSim->ResetBranchAddresses();
275  Points->Delete(); // maybe Points->Clear(); would be enough...
276  return 0;
277 
278 }
279 
280 
281 
282 
283 
284 
285 
286 
287 
288 int plotHoughTrackCand(PndFtsHoughTrackCand *trackCand, TGraph* returnTGraph)
289 {
290 
291  // calculate two points on lines, many points for parabola and return result in TGraph*
292 
293 
294  const Double_t fZLineParabola = trackCand->getZLineParabola();
295  const Double_t fZParabolaLine = trackCand->getZParabolaLine();
296  TVector3 pos;
297 
298 
299  const Int_t nPoints = 50; // first 2 for line before dipole, last 2 for line after dipole, rest for parabola within dipole
300  Double_t z[nPoints] = { fZLineParabola-500., fZLineParabola }; // all coordinates are in PANDA coordinate system
301  Double_t x[nPoints];
302  for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint)
303  {
304  switch (iPoint)
305  {
306  case 0:
307  z[iPoint] = fZLineParabola-100.;
308  break;
309  case 1:
310  z[iPoint] = fZLineParabola;
311  break;
312  case nPoints-2:
313  z[iPoint] = fZParabolaLine;
314  break;
315  case nPoints-1:
316  z[iPoint] = fZParabolaLine+100.;
317  break;
318  default:
319  z[iPoint] = (fZParabolaLine-fZLineParabola)/(nPoints-4);
320  }
321  pos = trackCand->getPos(z[iPoint]);
322  x[iPoint] = pos.X();
323  }
324 
325  returnTGraph->DrawGraph(nPoints,z,x,"LP,SAME") ;
326  return 0;
327 }
328 
329 
330 
331 
332 
333 
334 Double_t MakeHoughParabolaFitwithBfield(PndMultiField *fField, Double_t pzinv, Double_t theta,
335  Double_t weight, Double_t zOffset, Double_t hitshiftinx, Double_t zscaling, const Double_t c,
336  const Bool_t keepBConstant, TGraph* ParabolaFitInXZPlane)
337 {
338 
339  int verbose = 0;
340  // (zreal=zOffset, xreal=-hitshiftinx) = (zshifted = 0, xshifted = 0)
341  // zshifted = zreal - zOffset
342  // xshifted = xreal + hitshiftinx
343  // zreal = zshifted + zOffset
344  // xreal = xshifted - hitshiftinx
345 
346 
347 
348  // There is a bug for theta == 0
349  // calculate constants
350  // There is a bug for By == 0
351 
352  // 0 == theta is a special case which I did not want to implement. if statement is a workaround
353  if (0 == theta)
354  {
355  theta = 1. / 100.;
356  }
357 
358  const Int_t iterationsForBfield = 10; // 10 should be enough
359 
360  // for B field access
361  Double_t Byminus = 0., Byplus = 0.;
362  Double_t BMeanPlus = 0., BMeanMinus;
363  Double_t po[3], BB[3];
364  Double_t xplusshifted = 0., xminusshifted = 0., lastxminusshifted = 0., lastxplusshifted = 0., lastxminusreal= 0., lastxplusreal = 0., zshifted = 0.;
365  po[1] = 0.; // always use B-field for y == 0
366 
367  const Double_t meinpi = 3.14159265;
368  Double_t realtheta = theta / 360. * 2. * meinpi;
369  const Double_t tantheta = tan(realtheta);
370  const Double_t pzstuff = 1. / pzinv / 2. / c / tantheta / sin(realtheta);
371  const Int_t upperlimit = 500, lowerlimit = -200;
372  const Int_t nPoints = ceil((upperlimit-lowerlimit)/zscaling);
373  Double_t xplusreal[nPoints];
374  Double_t xminusreal[nPoints];
375  Double_t zreal[nPoints];
376  Int_t howManyPointsInBField = 0;
377 
378 
379 
380  if (2<verbose) { std::cout << "hitshiftinx is " << hitshiftinx << std::endl; }
381 
382  for (Int_t indexForArrays = 0, zshifted = lowerlimit; indexForArrays < nPoints; ++indexForArrays, zshifted = zshifted + 1. / zscaling)
383  {
384  // (zreal=zOffset, xreal=-hitshiftinx) = (zshifted = 0, xshifted = 0)
385  // zshifted = zreal - zOffset
386  // xshifted = xreal + hitshiftinx
387  // zreal = zshifted + zOffset
388  // xreal = xshifted - hitshiftinx
389 
390 
391  zreal[indexForArrays] = zshifted + zOffset;
392  // iterate to get the correct B field
393  for (Int_t iTry = 0; iTry<iterationsForBfield; ++iTry)
394  {
395 
396  if (kTRUE == keepBConstant)
397  {
398  // do not take B field into account for plotting
399  Byplus = 1.;
400  Byminus = 1.;
401  // do not iterate
402  iTry = iterationsForBfield+10;
403  }
404  else
405  {
406  // try to plot the track taking the B-field in y direction into account
407 
408  po[2] = zreal[indexForArrays];
409 
410  lastxminusreal = lastxminusshifted-hitshiftinx;
411  lastxminusreal = lastxminusshifted-hitshiftinx;
412 
413 
414  po[0] = lastxminusreal;
415  fField->GetFieldValue(po, BB); //return value in KG (G3)
416  Byminus = BB[1] / 10.; // y-component of magnetic field in Tesla
417 
418  po[0] = lastxplusreal;
419  fField->GetFieldValue(po, BB); //return value in KG (G3)
420  Byplus = BB[1] / 10.; // y-component of magnetic field in Tesla
421 
422  if (6<verbose)
423  {
424  std::cout << iTry << ". try: z+zOffset = " << zreal[indexForArrays] << " cm Byminus = "
425  << Byminus << " T Byplus = " << Byplus << " T"
426  << std::endl;
427  std::cout << "lastxplusreal = " << lastxplusreal << " cm lastxminusreal = " << lastxminusreal << " cm" << std::endl;
428  }
429 
430  // This is just a workaround as I have to divide by By (parabola becomes straight line if By is small)
431  if (0 == Byminus)
432  {
433  Byminus = 1. / 10000.;
434  }
435  if (0 == Byplus)
436  {
437  Byplus = 1. / 10000.;
438  }
439  }
440 
441  Double_t ztantheta = zshifted / tantheta;
442  Double_t a = -ztantheta + pzstuff / Byminus;
443  xminusshifted = a - sqrt(a * a - zshifted / pzinv / c / sin(realtheta) / Byminus - ztantheta * ztantheta);
444  a = -ztantheta + pzstuff / Byplus;
445  xplusshifted = a + sqrt(a * a - zshifted / pzinv / c / sin(realtheta) / Byplus - ztantheta * ztantheta);
446 
447 
448  lastxminusshifted = xminusshifted;
449  lastxplusshifted = xplusshifted;
450 
451  //std::cout << "\n\n(z+zOffset, xminus, xplus) = (" << z<< ", " << xminus << ", " << xplus << ")\n\n\n";
452  } // iterated to find correct B field
453 
454  xplusreal[indexForArrays] = xplusshifted - hitshiftinx;
455  xminusreal[indexForArrays] = xminusshifted - hitshiftinx;
456 
457 
458  if (zreal[indexForArrays] > 415. && zreal[indexForArrays] < 490.)
459  {
460  // TODO:: Try alternative: determine maximum of B field
461 
462  //the following lines average the B field (probably not what I want to do...)
463  ++howManyPointsInBField;
464  // integrate B field along the path of the particle (might be helpful for reconstructing pz
465  po[2] = zreal[indexForArrays];
466  po[0] = xplusreal[indexForArrays];
467  fField->GetFieldValue(po, BB); //return value in KG (G3)
468  BMeanPlus += BB[1]/10.;
469 
470  po[0] = xminusreal[indexForArrays];
471  fField->GetFieldValue(po, BB); //return value in KG (G3)
472  BMeanMinus += BB[1]/10.;
473  }
474 
475  } // for loop over z coordinate
476 
477  BMeanPlus = BMeanPlus/howManyPointsInBField;
478  BMeanMinus = BMeanMinus/howManyPointsInBField;
479 
480  if (1<verbose)
481  {
482  cout << "??? BMeanPlus = " << BMeanPlus << std::endl;
483  cout << "??? BMeanMinus = " << BMeanMinus << std::endl;
484  }
485 
486  ParabolaFitInXZPlane->DrawGraph(nPoints, zreal, xplusreal,"P,SAME");
487  ParabolaFitInXZPlane->DrawGraph(nPoints, zreal, xminusreal,"P,SAME");
488  return max(BMeanPlus,BMeanMinus);
489 }
490 
491 
492 // This procedure draws the magnetic field in pointer fField into histogram Bfield
493 int DrawField(PndMultiField *fField, TH2F* Bfield)
494 {
495  int verbose = 0;
496  // For magnetic field access
497  Double_t po[3], BB[3];
498  Double_t Btot = 0.;
499 
500  po[1] = 0.; // set y position to 0 (because FTS currently does not provide a y position
501 
502  for (Double_t z = -199.; z < 800.; z += 2.)
503  {
504  po[2] = z;
505  for (Double_t x = -249.; x < 250.; x += 2.)
506  {
507  // cout << "X = " << x << endl;
508  po[0] = x;
509  fField->GetFieldValue(po, BB); //return value in KG (G3)
510  if (7<verbose)
511  {
512  std::cout << "B_y = " << BB[1]/10. << " T for Z = " << po[2] << " cm X = " << po[0] << " cm" << std::endl;
513  }
514  Btot = TMath::Sqrt(BB[0] * BB[0] + BB[1] * BB[1] + BB[2] * BB[2]);
515  Bfield->Fill(z, x, Btot / 10.); // This plots B_total
516  //Bfield->Fill(z,x,BB[1]/ 10.); // This plots By = y-component of B field
517  }
518  }
519  return 0;
520 }
521 
522 
523 
524 
525 
526 
527 
528 
529 
530 
531 
532 
533 
534 
535 
536 
537 
538 
539 
540 
541 
542 
543 
544 
545 
546 
547 
548 
549 
550 
551 
552 
553 
554 
555 
556 
557 
558 
559 
560 
561 
562 
563 
564 // -----------------------------------------------------------------------
565 // -----------------------------------------------------------------------
566 // ------------------------ MACRO STARTS HERE --------------------
567 // -----------------------------------------------------------------------
568 // -----------------------------------------------------------------------
569 int plotTrackCands(int iEvent)
570 {
571 
572  UInt_t fVerbose = 0;
573 
574 
575  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
576  // basiclibs(); // this gives an error when compiling
577 
578  // Load the example libraries
579  gSystem->Load("libGeoBase");
580  gSystem->Load("libParBase");
581  gSystem->Load("libBase");
582  gSystem->Load("libField");
583  gSystem->Load("libFts");
584 
585 
586 
587  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
588  gROOT->SetStyle("Plain");
589  gStyle->SetMarkerStyle(20);
590  gStyle->SetOptStat(0);
591  gStyle->SetPalette(1);
592 
593 
594 
595  // -----------------------------------------------------------------------
596  // FileNames
597  TString RecoFile = "reco_complete.root";
598  TString DigiFile = "digi_complete.root";
599  TString SimFile = "sim_complete.root";
600  // -----------------------------------------------------------------------
601 
602 
603 
604 
605 
606 
607  // General stuff
608  const Double_t BeamMomentum = 15.; //6.991;
609  const Double_t meinpi = 3.14159265359;
610  const Int_t ResolutionX = 800, ResolutionY = 600; // for plotting
611 
612  // -----------------------------------------------------------------------
613 
614  // for B field access
615  Double_t By = 0.;
616  Double_t po[3], BB[3];
617  Double_t BMeanForParabolapz = 0;
618  Double_t BMeanForParabola = 0;
619 
620  // -----------------------------------------------------------------------
621 
622 
623 
624 
625 
626 
627  // Histograms
628  // TH2D* hits = new TH2D("hits", "hits", 276, 250, 800, 331, -330, 330);
629  // TString mcPointsHistoname = "mcPoints event ";
630  // mcPointsHistoname += iEvent;
631  // TH2D* mcPoints = new TH2D(mcPointsHistoname, mcPointsHistoname, 1200, -200, 1000, 500, -250,
632  // 250);
633 
634  TH2F* fFieldHist = NULL;
635 
636  fFieldHist = new TH2F("Bfield", "", 650, -200, 1100, 250, -250, 250);
637  fFieldHist->GetXaxis()->SetTitle("z [cm]");
638  fFieldHist->GetYaxis()->SetTitle("x [cm]");
639 
640 
641 
642  // TH2D* Fitinxzplane= new TH2D("Fitinxzplane", "Fitinxzplane", 800, 0, 800, 500, -250, 250);
643  TGraph* FitInXZPlane = new TGraph();
644  FitInXZPlane->SetMarkerColor(4);
645  FitInXZPlane->SetMarkerStyle(1);
646  FitInXZPlane->SetMarkerSize(1.1);
647  FitInXZPlane->SetLineColor(4);
648  FitInXZPlane->SetLineWidth(1.5);
649 
650 
651 
652 
653 
654  // Initialize the B-field // new method
655  PndMultiField *fField= new PndMultiField("FULL", BeamMomentum); //beam momentum is BeamMomentum in GeV/c (for scaling of dipole field)
656  fField->Init();
657 
658 
659 
660 
661  // // That is what I used before Mohammad did some changes
662  //
663  // // Create the Simulation run manager--------------------------------
664  // FairRunSim *fRun = new FairRunSim();
665  // //fRun->SetName(SimEngine.Data() );
666  // //fRun->SetOutputFile(OutputFile.Data());
667  // fRun->SetBeamMom(BeamMomentum);
668  // //fRun->SetMaterials(MediaFile.Data());
669  // //FairRuntimeDb *rtdb=fRun->GetRuntimeDb();
670  //
671  // PndMultiField *fField= new PndMultiField("FULL");
672  // fField->Init();
673 
674 
675  TCanvas *cEventZx = NULL;
676 
677  cEventZx = new TCanvas("zx plane", "zx plane", ResolutionX, ResolutionY);
678  // Plot the B-field
679  DrawField(fField, fFieldHist);
680  fFieldHist->Draw("cont1");
681 
682 
683 
684 
685 
686 
687 
688 
689 
690 
691 
692 
693 
694 
695 
696 
697 
698 
699 
700 
701 
702 
703 
704 
705 
706 
707 
708 
709  // // This code can be used to test the fit drawing method with B field
710  // for (Double_t thetatest=0.1; thetatest <30.; thetatest+=14.)
711  // {
712  // void MakeHoughParabolaFitwithBfield(PndMultiField *fField, Double_t pzinv, Double_t theta,
713  // Double_t weight, Double_t zOffset, Double_t hitshiftinx, Double_t zscaling, const Double_t c,
714  // const Bool_t keepBConstant, TH2D* Fitinxzplane)
715  // // make some test plots for 1/pz
716  // //MakeHoughParabolaFitwithBfield(fField,-1./300./thetatest,0.01,thetatest*thetatest,zOffset,0.,1.,c,keepBConstant,Fitinxzplane);
717  // // rotate in positive theta direction (counterclockwise)
718  // //MakeHoughParabolaFitwithBfield(fField,-1./1000.,thetatest,thetatest*thetatest,zOffset,0.,1.,c,keepBConstant,Fitinxzplane);
719  // // rotate in negative theta direction (clockwise)
720  // MakeHoughParabolaFitwithBfield(fField,-1./3000.,-thetatest,thetatest*thetatest,zOffset,0.,1.,c,keepBConstant,Fitinxzplane);
721  // }
722  //
723  //
724  // // This code can be used to test the fit drawing method without B field
725  // for (Double_t thetatest=0.1; thetatest <30.; thetatest+=14.)
726  // {
727  // // make some test plots for 1/pz
728  // //MakeHoughFit(-1./300./thetatest,0.01,thetatest*thetatest,zOffset,0.,1.,c,keepBConstant,Fitinxzplane);
729  // // rotate in positive theta direction (counterclockwise)
730  // //MakeHoughFit(-1./1000.,thetatest,thetatest*thetatest,zOffset,0.,1.,c,keepBConstant,Fitinxzplane);
731  // // rotate in negative theta direction (clockwise)
732  // MakeHoughFit(-1./3000.,-thetatest,thetatest*thetatest,zOffset,0.,1.,c,keepBConstant,Fitinxzplane);
733  // }
734 
735  // plot the fit
736  //const int ResolutionX = 800, ResolutionY = 600;
737  //TCanvas *cHoughFit=new TCanvas("cHoughFit","cHoughFit",ResolutionX,ResolutionY);
738 
739 
740  //
741 
742 
743 
744 
745 
746 
747 
748 
749 
750 
751 
752 
753 
754 
755 
756 
757 
758 
759 
760 
761 
762 
763 
764 
765 
766 
767 
768 
769 
770 
771  cEventZx->cd(0);
772  // Plot the MC Points
773  const int subdetectorCount = 11; //
774  TString branchNames[subdetectorCount] =
775  { "FTSPoint", "STTPoint", "MVDPoint", "GEMPoint", "DrcBarPoint",
776  "DrcPDPoint", "DskParticle", "PndDskFLGHit", "MdtPoint",
777  "FtofPoint", "EmcHit" };
778  TString tcaNames[subdetectorCount] =
779  { "PndFtsPoint", "PndSttPoint", "PndSdsMCPoint", "PndGemMCPoint",
780  "PndDrcBarPoint", "PndDrcPDPoint", "PndDskParticle", "PndDskFLGHit",
781  "PndMdtPoint", "PndFtofPoint", "PndEmcHit" };
782  std::cout << "Event number " << iEvent << std::endl;
783 
784 
785  TGraph *mcPoints = new TGraph();
786  mcPoints->SetMarkerColor(2);
787  mcPoints->SetMarkerStyle(29);
788  mcPoints->SetMarkerSize(1.4);
789  mcPoints->SetLineColor(2);
790  TFile* fSim = new TFile(SimFile.Data());
791  for (int iSubdetector = 0; iSubdetector < subdetectorCount; ++iSubdetector)
792  {
793  cout << iSubdetector
794  << ". subdetector: branchNames[iSubdetector], tcaNames[iSubdetector] = "
795  << branchNames[iSubdetector] << " " << tcaNames[iSubdetector]
796  << std::endl;
797  PlotMCTracks(iEvent, branchNames[iSubdetector], tcaNames[iSubdetector], mcPoints, fSim);
798  }
799  // do not close the file fSim, otherwise, you won't see any plots
800 
801 
802 
803  // get FTS Hits
804  TFile* fDigi = new TFile(DigiFile.Data());
805  TTree* tDigi = (TTree*) (fDigi->Get("pndsim"));
806  //tDigi->StartViewer();
807 
808  TClonesArray* fFtsHitArray = new TClonesArray("PndFtsHit");
809  tDigi->SetBranchAddress("FTSHit", &fFtsHitArray);
810  tDigi->GetEntry(iEvent); // Look at event iEvent
811 
812  // Plot FTS hits
813  TGraph* hits = new TGraph();
814  hits->SetMarkerStyle(2);
815  hits->SetMarkerColor(1);
816  hits->SetMarkerSize(2.0);
817  const Int_t nHits = fFtsHitArray->GetEntriesFast();
818  Double_t hitX[nHits];
819  Double_t hitY[nHits];
820  Double_t hitZ[nHits];
821  for (int iHit = 0; iHit < nHits; iHit++)
822  {
823  PndFtsHit* myHit = (PndFtsHit*) fFtsHitArray->At(iHit);
824  // if (kTRUE == onlyUseHitsFromNonSkewedStraws)
825  // {
826  if (0 != myHit->GetSkewed())
827  {
828  if (0<fVerbose) {std::cout << "Not drawing hit with index " << iHit << " , because it comes from a skewed straw! LayerID = " << myHit->GetLayerID() << std::endl;}
829  continue;
830  }
831 
832  // }
833  TVector3 hit;
834  myHit->Position(hit);
835 
836  hitX[iHit] = hit.X();
837  hitY[iHit] = hit.Y();
838  hitZ[iHit] = hit.Z();
839 
840 
841  if (0<fVerbose)
842  {
843  std::cout << "Hit " << iHit << " : ";
844  //myHit->Print();
845  std::cout << "X: " << hitX[iHit] << " Y: " << hitY[iHit] << " Z: " << hitZ[iHit]
846  << std::endl;
847  }
848  }
849  hits->DrawGraph(nHits, hitZ, hitX, "P,SAME");
850 
851 
852 
853  // TTree *tree=(TTree *) fSim->Get("pndsim") ;
854  // TClonesArray* mc_array=new TClonesArray("PndMCTrack");
855  // tree->SetBranchAddress("MCTrack", &mc_array);
856  // tree->GetEntry(iEvent);
857  // for (Int_t mcIndex = 0; mcIndex < mc_array->GetEntriesFast(); ++mcIndex)
858  // {
859  // PndMCTrack *mctrack = (PndMCTrack*)mc_array->At(mcIndex);
860  // // only look at primary particles
861  // if (mctrack->GetMotherID()!=-1) continue;
862  //
863  // Int_t mc_pid = mctrack->GetPdgCode();
864  // Double_t mc_px = mctrack->GetMomentum().X();
865  // Double_t mc_py = mctrack->GetMomentum().Y();
866  // Double_t mc_pz = mctrack->GetMomentum().Z();
867  // // Double_t mc_theta = mctrack->GetMomentum().Theta()*TMath::RadToDeg();
868  // // Double_t mc_phi = mctrack->GetMomentum().Phi()*TMath::RadToDeg();
869  // //
870  // // Double_t mc_theta_xz = atan2(mc_px,mc_pz)*TMath::RadToDeg();
871  //
872  // std::cout << "Primary particle with pid " << mc_pid << " has momentum = (" << mc_px << ", " << mc_py << ", " << mc_pz << ") GeV/c (at interaction point)" << std::endl;
873  //
874  // // std::cout<<setprecision(4);
875  // TString histoname = "Single #mu^{-} event with MC p_{z} = ";
876  // // manually round mc_pz to a.bc, for example: 4.71
877  // TString numMCTruth="", numReco="";
878  // numMCTruth += mc_pz+0.0005;
879  // numMCTruth.Resize(5);
880  // histoname += numMCTruth;
881  // histoname += " GeV/c <-> FTS PR p_{z} = ";
882  // numReco +=pzinvpeak+0.0005;
883  // numReco.Resize(5);
884  // histoname += numReco;
885  // histoname += " GeV/c"; //, red = points, black = FTS hits";
886  //
887  // fFieldHist->SetTitle(histoname);
888 
889 
890 
891  // get PndFtsHoughTrackCand
892  TFile* fReco = new TFile(RecoFile.Data());
893  TTree* tReco = (TTree*) (fReco->Get("pndsim"));
894  //tReco->StartViewer();
895 
896  TClonesArray* fTrackCands = new TClonesArray("PndFtsHoughTrackCand");
897  tReco->SetBranchAddress("PndFtsHoughTrackCand", &fTrackCands);
898  tReco->GetEntry(iEvent); // Look at event iEvent
899  cout << fTrackCands->GetEntriesFast()<< " Hough track cands. in event\n";
900 
901  // Plotting of HoughTrackCands
902  cEventZx->cd(0);
903  for (UInt_t iTrackCand=0; iTrackCand < fTrackCands->GetEntriesFast(); ++iTrackCand){
904  PndFtsHoughTrackCand *myTrackCand = (PndFtsHoughTrackCand*) fTrackCands->At(iTrackCand);
905  plotHoughTrackCand(myTrackCand, FitInXZPlane);
906  }
907  return 0;
908 }
909 
910 
911 #endif
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
fField
PndMultiField * fField
Definition: sim_emc_apd.C:97
fVerbose
int fVerbose
Definition: poormantracks.C:24
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
RecoFile
TString RecoFile
Definition: emc/run_linkCreation.C:15
PlotMCTracksPrintBField
int PlotMCTracksPrintBField(PndMultiField *fField, int iEvent, TString branchName, TString tcaName, TH2D *mcPoints)
Definition: plotTrackCands.C:166
DrawField
int DrawField(PndMultiField *fField, TH2F *Bfield)
Definition: plotTrackCands.C:493
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndMultiField::GetFieldValue
void GetFieldValue(const Double_t point[3], Double_t *bField)
Definition: PndMultiField.cxx:257
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
BeamMomentum
double BeamMomentum
Definition: sim_ftof_stof.C:17
Y
double Y
Definition: anaLmdDigi.C:68
PndFtsHoughTrackCand::getPos
TVector3 getPos(const Double_t zLabSys) const
Definition: PndFtsHoughTrackCand.cxx:331
PndFtsHit::GetLayerID
Int_t GetLayerID() const
Definition: PndFtsHit.h:74
MakeHoughParabolaFitwithBfield
Double_t MakeHoughParabolaFitwithBfield(PndMultiField *fField, Double_t pzinv, Double_t theta, Double_t weight, Double_t zOffset, Double_t hitshiftinx, Double_t zscaling, const Double_t c, const Bool_t keepBConstant, TGraph *ParabolaFitInXZPlane)
Definition: plotTrackCands.C:334
DigiFile
TString DigiFile
Definition: emc/run_linkCreation.C:14
plotTrackCands
int plotTrackCands()
Definition: plotMomRes.C:55
a
Int_t a
Definition: anaLmdDigi.C:126
nHits
int nHits
Definition: RiemannTest.C:16
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
z
Double_t z
Definition: createRootFscGeometryFile.C:146
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
X
double X
Definition: anaLmdDigi.C:68
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Z
double Z
Definition: anaLmdDigi.C:68
PndFtsHit::GetSkewed
Int_t GetSkewed() const
Definition: PndFtsHit.h:75
PndFtsHoughTrackCand
Class for saving a FTS track cand. for Hough transform based FTS PR.
Definition: PndFtsHoughTrackCand.h:47
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
plotHoughTrackCand
int plotHoughTrackCand(PndFtsHoughTrackCand *trackCand, TGraph *returnTGraph)
Definition: plotTrackCands.C:288
PlotMCTracks
int PlotMCTracks(int iEvent, TString branchName, TString tcaName, TGraph *mcPoints, TFile *fSim=0)
Definition: plotTrackCands.C:75
point
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
Generated on Wed Apr 3 2019 10:02:48 for FairRoot/PandaRoot by   doxygen 1.8.5