FairRoot/PandaRoot
Functions
geo/Check_particle_path.C File Reference
#include <TString.h>
#include <TH1.h>
#include <TH2.h>
#include <sstream>
#include <iostream>
#include <vector>
#include <map>
#include <TVectorT.h>
#include <TCanvas.h>
#include <TApplication.h>
#include <TChain.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TClonesArray.h>
#include <TRotation.h>
#include <TFile.h>
#include <TGaxis.h>
#include <PndMCTrack.h>
#include <PndSdsMCPoint.h>
#include <FairRunAna.h>
#include <FairGeane.h>
#include <FairRtdbRun.h>
#include <FairRuntimeDb.h>
#include <FairParRootFileIo.h>
#include <TGraph.h>
#include <TF1.h>
#include <algorithm>
#include <fstream>
#include <TStyle.h>
#include <TColor.h>

Go to the source code of this file.

Functions

const int nplanes (120)
 
double last_percent (-1.)
 
void DrawProgressBar (int len, double percent)
 
Double_t function_beampipe (Double_t *x, Double_t *par)
 
int Check_particle_path ()
 
int main ()
 

Function Documentation

int Check_particle_path ( )

Definition at line 49 of file geo/Check_particle_path.C.

References DigiFile, Double_t, DrawProgressBar(), fRun, function_beampipe(), PndMCTrack::GetMomentum(), PndMCTrack::GetPdgCode(), PndSdsMCPoint::GetPosition(), PndSdsMCPoint::GetSensorID(), PndMCTrack::GetStartVertex(), if(), PndMCTrack::IsGeneratorCreated(), nEvents, nHits, nlines, nplanes(), out, parFile, parInput1, rtdb, startEvent, storePath, TString, x, and z.

Referenced by main().

49  {
50  cout << " analysis tool vers. 1.3 " << endl;
51 
52  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
53  gSystem->Load("libSds");
54  gSystem->Load("libSdsReco");
55  gSystem->Load("libLmd");
56  gSystem->Load("libLmdReco");
57  gSystem->Load("libLmdTrk");
58 
59  TPad foo; // never remove this line :-)))
60  if(1){
61  gROOT->SetStyle("Plain");
62  const Int_t NRGBs = 5;
63  const Int_t NCont = 255;
64  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
65  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
66  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
67  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
68  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
69  gStyle->SetNumberContours(NCont);
70  gStyle->SetTitleFont(10*13+2,"xyz");
71  gStyle->SetTitleSize(0.06, "xyz");
72  gStyle->SetTitleOffset(0.8,"y");
73  gStyle->SetTitleOffset(1.3,"z");
74  gStyle->SetLabelFont(10*13+2,"xyz");
75  gStyle->SetLabelSize(0.06,"xyz");
76  gStyle->SetLabelOffset(0.009,"xyz");
77  gStyle->SetPadBottomMargin(0.16);
78  gStyle->SetPadTopMargin(0.16);
79  gStyle->SetPadLeftMargin(0.10);
80  gStyle->SetPadRightMargin(0.10);
81  gStyle->SetOptTitle(1);
82  gStyle->SetOptStat(1);
83  gROOT->ForceStyle();
84  gStyle->SetFrameFillColor(0);
85  gStyle->SetFrameFillStyle(0);
86  TGaxis::SetMaxDigits(3);
87  }
88 
89  // ------------------------------------------------------------------------
90 
91  TString storePath = "./";
92  TString startEvent = "0";
93  // ---- Input files --------------------------------------------------------
94  TString simMC = storePath + "/Lumi_MC_";
95  simMC += startEvent;
96  simMC += ".root";
97  TChain tMC("pndsim");
98  tMC.Add(simMC);
99 
100  TString DigiFile = storePath + "/Lumi_digi_";
101  DigiFile += startEvent;
102  DigiFile += ".root";
103  TChain tdigiHits("pndsim");
104  tdigiHits.Add(DigiFile);
105 
106  // ---- Output file ----------------------------------------------------------------
107  TString fakefile = storePath + "/fake";
108  fakefile += startEvent;
109  fakefile += ".root";
110 
111  // ---- Output file ----------------------------------------------------------------
112  TString out = storePath + "/beampipe_results";
113  out += startEvent;
114  out += ".root";
115  TFile *output_histfile = new TFile(out, "RECREATE");
116  // ---------------------------------------------------------------------------------
117 
118  // Parameter file << needed for geane back tracking
119  TString parFile = storePath+"/Lumi_Params_";
120  parFile += startEvent;
121  parFile += ".root";
122 
123  // ----- Reconstruction run -------------------------------------------
124  FairRunAna *fRun = new FairRunAna();
125  fRun->SetInputFile(simMC);
126  fRun->AddFriend(DigiFile);
127  //fRun->AddFriend(RecoFile);
128  //fRun->AddFriend(CandFile);
129  //fRun->AddFriend(TrkFile);
130  fRun->SetOutputFile(fakefile);
131  // ------------------------------------------------------------------------
132 
133  // ----- Parameter database --------------------------------------------
134  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
135  FairParRootFileIo* parInput1 = new FairParRootFileIo(kTRUE);
136  parInput1->open(parFile.Data());
137  rtdb->setFirstInput(parInput1);
138 
139  //--- MC info -----------------------------------------------------------------
140  TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
141  tMC.SetBranchAddress("MCTrack", &true_tracks); //True Track to compare
142 
143  TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
144  tMC.SetBranchAddress("LMDPoint", &true_points); //True Points to compare
145 
146  //----------------------------------------------------------------------------------
147 
148  // histograms per plane
149  //TH2* hist_xy [nplanes];
150  //TH1* hist_theta_in [nplanes];
151 
152  TH2* hist_xz = new TH2F("hist_xz", "hist_xz", 150, 0, 150, 1000, -40, 40);
153  TH2* hist_yz = new TH2F("hist_yz", "hist_yz", 150, 0, 150, 1000, -40, 40);
154  TH2* hist_dxdz_z = new TH2F("hist_dxdz_z", "hist_dxdz_z", 150, 0, 150,
155  60000, -1e-3, 51.e-3);
156  TH2* hist_dydz_z = new TH2F("hist_dydz_z", "hist_dydz_z", 150, 0, 150,
157  1000, -5e-3, 5e-3);
158  TH2* hist_mom_z = new TH2F("hist_mom_z", "hist_mom_z", 150, 0, 150,
159  3000, -1.e-5, 1.e-5);
160  TH1* hist_n_events = new TH1F("hist_n_events", "hist_n_events", 150, 0, 150);
161 
162  //TCanvas temp_canvas("temp_canvas", "canvas for initialization", 100, 100);
163  //temp_canvas.cd();
164  // loop over planes
165  //cout << " constructing histograms for " << nplanes << " planes with " << nsensors_per_plane << " sensors per plane " << endl;
166 
167  /*
168  for (int iplane = 0; iplane < nplanes; iplane++){
169  stringstream hist_name;
170  stringstream hist_title;
171 
172  hist_name << "hist_xy_plane_" << iplane;
173  hist_title << "xy hit distribution plane " << iplane;
174  hist_xy[iplane] = new TH2F(hist_name.str().c_str(), hist_title.str().c_str(),
175  100, -40, 40, 100, -40, 40);
176  hist_xy[iplane]->Draw();
177  hist_xy[iplane]->GetXaxis()->SetTitle("X [cm]");
178  hist_xy[iplane]->GetYaxis()->SetTitle("Y [cm]");
179 
180  hist_name .str("");
181  hist_title.str("");
182 
183  hist_name << "hist_theta_init_plane_" << iplane;
184  hist_title << "initial #Theta distribution plane " << iplane;
185  hist_theta_in[iplane] = new TH1F(hist_name.str().c_str(), hist_title.str().c_str(),
186  150, 0, 15e-3);
187  hist_theta_in[iplane]->Draw();
188  hist_theta_in[iplane]->GetXaxis()->SetTitle("#Theta [rad]");
189  hist_theta_in[iplane]->GetYaxis()->SetTitle("entries");
190 
191  }*/
192 
193  int nEvents = tMC.GetEntries();
194  cout << " reading " << nEvents << " Events " << endl;
195  // store the x and z measurements for a later fit of a graph
196  vector<double> pos_x;
197  vector<double> pos_z;
198  map<int,int> hits_z; // count the hits per plane
199  map<int,TH2*> hists_ang_acceptance; // angular acceptance as a function of z
200  map<int,TH2*>::iterator hists_ang_acceptance_it, hists_ang_acceptance_itnorm;
201  // create the histogram to normalize to.
202  hists_ang_acceptance[-1] = new TH2F("hists_ang_acceptance", "generated angular distribution",200, -10, 10, 200, -10, 10);//, 200, 2., 10., 150., -3.141, 3.141);
203  hists_ang_acceptance_itnorm = hists_ang_acceptance.find(-1);
204  // optionally read existing values in order to create a mean graph
205  // for several scenarios. If you want to combine values, put the
206  // x_z_values_out.dat into the next folder to analyze and name it
207  // there as x_z_values_in.dat
208  ifstream x_z_values_in("x_z_values_in.dat");
209  if (nEvents <= 1000 && x_z_values_in.is_open()){
210  cout << " ************ reading x_z_values_in.dat ************ " << endl;
211  double x, z;
212  int nlines(0);
213  while(x_z_values_in.good()){
214  x_z_values_in >> x >> z;
215  if (x_z_values_in.good()){
216  pos_x.push_back(x);
217  pos_z.push_back(z);
218  nlines++;
219  } else break;
220  }
221  x_z_values_in.close();
222  cout << " read " << nlines << " lines "<< endl;
223  }
224 
225  for (Int_t j = 0; j < nEvents ; j++) {
226  DrawProgressBar(50, (j + 1) / ((double) nEvents));
227  tMC.GetEntry(j);
228  tdigiHits.GetEntry(j);
229  const int nParticles = true_tracks->GetEntriesFast();
230 
231  int npbar = 0;
232  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
233  // if(nParticles!=4) continue;
234  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(iParticle);
235  Int_t mcID = mctrk->GetPdgCode();
236  if (mctrk->IsGeneratorCreated()) {//mcID == -2212) {
237  npbar++;
238  // if(nParticles!=4) cout<<"For event #"<<j<<" mcID="<<mcID<<endl;
239  TVector3 momMC = mctrk->GetMomentum();
240  hists_ang_acceptance[-1]->Fill(momMC.X()/momMC.Z()*1e3, momMC.Y()/momMC.Z()*1e3);
241  //hists_ang_acceptance[-1]->Fill(MomMC.Theta()*1e3, MomMC.Phi());
242  }
243  }
244  if (0 && npbar != 1)
245  cout << "found " << npbar << " anti-protons" << endl;
246 
247  int nHits = true_points->GetEntriesFast();
248  int nSdsHits = 0;
249  for (Int_t iHit = 0; iHit < nHits; iHit++) {
250  PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
251  if (mcpoint->GetTrackID() < 0) continue;
252  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(
253  mcpoint->GetTrackID());
254  if (mctrk->IsGeneratorCreated()) { // take only anti-protons from the generator
255  if (mcpoint) {
256  TVector3 _mcpoint((mcpoint->GetPosition())); // position at the sensor entrance
257  TVector3 _mctrack(mcpoint->GetPx(), mcpoint->GetPy(),
258  mcpoint->GetPz()); // momentum of the track at the entrance
259 
260  TVector3 momMC = mctrk->GetMomentum(); // momentum in the primary vertex
261  TVector3 posMC = mctrk->GetStartVertex(); // position of the primary vertex
262 
263  //det_id = mcpoint->GetDetectorID(); // store the detector id
264  //sens_id = mcpoint->GetSensorID(); // store the sensor id
265  // calculate the plane and sensor on this plane
266  int plane = mcpoint->GetSensorID();
267 
268  if (plane < 0 || plane >= nplanes) {
269  cout
270  << " warning: calculated plane or sensor is wrong: plane "
271  << plane << endl;
272  } else {
273  // is it an unknown plane so far?
274  if (hists_ang_acceptance.find(plane) == hists_ang_acceptance.end()){
275  stringstream hist_name;
276  stringstream hist_title;
277  hist_name << "hists_ang_acceptance_plane_" << plane;
278  hist_title << "angular acceptance at z = " << plane*10. << " cm " << endl;
279  hists_ang_acceptance[plane] = new TH2F(hist_name.str().c_str(), hist_title.str().c_str(),
280  200, -10, 10, 200, -10, 10);//200, 2., 10., 150., -3.141, 3.141);
281  }
282  hists_ang_acceptance[plane]->Fill(momMC.X()/momMC.Z()*1e3, momMC.Y()/momMC.Z()*1e3);//(momMC.Theta()*1e3, momMC.Phi());
283  if (j < 1000){
284  pos_x.push_back(_mcpoint.X());
285  pos_z.push_back(_mcpoint.Z());
286  }
287  hist_xz->Fill(plane, _mcpoint.X());
288  hist_yz->Fill(plane, _mcpoint.Y());
289  hist_dxdz_z->Fill(plane, _mctrack.X() / _mctrack.Z());
290  hist_dydz_z->Fill(plane, _mctrack.Y() / _mctrack.Z());
291  hist_mom_z->Fill(plane, (_mctrack.Mag()-momMC.Mag())/momMC.Mag());
292  //hist_xy [plane] ->Fill(x_in,y_in);
293  //hist_theta_in [plane] ->Fill(ptheta_in);
294  hits_z[plane]++;
295  }
296  nSdsHits++;
297  }
298  }
299  }
300  if (0 && nHits > 0)
301  cout << "found " << nHits << " hit(s) with " << nSdsHits
302  << " sds hit(s)" << endl;
303  }
304 
305  ofstream x_z_values_out("x_z_values_out.dat");
306  if (nEvents <= 1000 && x_z_values_out.is_open()){
307  cout << " ************ writing x_z_values_out.dat ************ " << endl;
308  for (unsigned int iline = 0; iline < pos_x.size(); iline++){
309  x_z_values_out << pos_x[iline] << " " << pos_z[iline] << "\n";
310  }
311  x_z_values_out.close();
312  cout << " wrote " << pos_x.size() << " lines "<< endl;
313  }
314 
315  //std::sort(pos_x.begin(), pos_x.end());
316  //std::sort(pos_z.begin(), pos_z.end());
317  TGraph* graph_x_z = new TGraph( pos_x.size(), (double*)&pos_z[0], (double*)&pos_x[0]);
318  graph_x_z->SetNameTitle("graph_x_z","x over z hits");
319  //graph_x_z->SetMarkerStyle(4);
320  graph_x_z->SetLineWidth(0);
321  TF1* func_beampipe = new TF1("func_beampipe", function_beampipe, 0., pos_z[pos_z.size()-1], 3);
322  // par[0]: where bending starts (dipole region)
323  // par[1]: the bending radius
324  // par[2]: the bending angle
325  func_beampipe->SetParameters(355., 6000., 40.e-3);
326  func_beampipe->SetParNames(
327  "bending start [cm] ",
328  "bending radius [cm] ",
329  "bend by phi [rad] ");
330  //func_beampipe->FixParameter(2,40.068e-3);
331  gStyle->SetOptFit(111);
332 
333  TGaxis::SetMaxDigits(3);
334  TCanvas* canvas1 = new TCanvas("canvas1", "x over z hits", 1200, 800);
335  //canvas1->Divide(2,2);
336  canvas1->cd(0);
337  //hist_xz->Draw("COLZ");
338  graph_x_z->Draw("AP");
339  graph_x_z->Fit(func_beampipe);
340  graph_x_z->GetXaxis()->SetTitle("z [cm]");
341  graph_x_z->GetYaxis()->SetTitle("x [cm]");
342  func_beampipe->Draw("same");
343  cout << " x displacement at 11 m is " << func_beampipe->Eval(1100.) << endl;
344  gPad->Update();
345  //func_beampipe->Draw("same");
346  canvas1->Print("beampipe_results.png");
347  // create residual distribution
348  TCanvas* canvas6 = new TCanvas("canvas6", "residuals", 1200, 800);
349  TH2F* hist_res_x_z = new TH2F("hist_res_x_z", "residuals of x over z fit", 150, 0, 1500, 100, -0.1, 0.1);
350  for (unsigned int iz = 0; iz < pos_z.size(); iz++){
351  hist_res_x_z->Fill(pos_z[iz], pos_x[iz]-func_beampipe->Eval(pos_z[iz]));
352  }
353  hist_res_x_z->Draw("COLZ");
354  hist_res_x_z->GetXaxis()->SetTitle("x [cm]");
355  hist_res_x_z->GetYaxis()->SetTitle("x - x_{fit} [cm]");
356  canvas6->Print("beampipe_results.ps(");
357 
358  TCanvas* canvas2 = new TCanvas("canvas2", "canvas2", 1200, 800);
359  //canvas2->Divide(2,2);
360  canvas2->cd(0);
361  hist_yz->Draw("COLZ");
362  canvas2->Print("beampipe_results.ps(");
363  TCanvas* canvas3 = new TCanvas("canvas3", "canvas3", 1200, 800);
364  //canvas3->Divide(2,2);
365  canvas3->cd(0);
366  hist_dxdz_z->Draw("COLZ");
367  canvas3->Print("beampipe_results.ps(");
368  TCanvas* canvas4 = new TCanvas("canvas4", "canvas4", 1200, 800);
369  //canvas4->Divide(2,2);
370  canvas4->cd(0);
371  hist_dydz_z->Draw("COLZ");
372  canvas4->Print("beampipe_results.ps(");
373  TCanvas* canvas5 = new TCanvas("canvas5", "canvas5", 1200, 800);
374  //canvas5->Divide(2,2);
375  canvas5->cd(0);
376  hist_mom_z->Draw("COLZ");
377  canvas5->Print("beampipe_results.ps(");
378  TCanvas* canvas7 = new TCanvas("canvas6", "canvas6", 1200, 400);
379  gPad->SetLeftMargin(0.1);
380  gPad->SetRightMargin(0.05);
381  canvas7->cd(0);
382  // check the acceptance of the beam pipe as a function of z
383  // normalize to the plane before
384  map<int, int>::iterator hits_z_it, hits_z_it_next;
385  for (hits_z_it = hits_z.begin(); hits_z_it != hits_z.end(); hits_z_it++){
386  hits_z_it_next = hits_z_it;
387  hits_z_it_next++;
388  if (hits_z_it_next == hits_z.end()) break;
389  if (hits_z_it->second == 0) continue;
390  if (hits_z_it_next->first != hits_z_it->first+1){
391  cout << " sorry, entries are not sorted! " << hits_z_it_next->first << " != " << hits_z_it->first << " +1 "<< endl;
392  continue;
393  }
394  hist_n_events->Fill(hits_z_it_next->first, ((double)hits_z_it_next->second)/((double)hits_z_it->second));
395  }
396  hist_n_events->GetXaxis()->SetTitle("z [dm]");
397  hist_n_events->GetYaxis()->SetTitle("relative acceptance");
398  hist_n_events->Draw("hist text");
399  canvas7->Print("beampipe_results.ps(");
400 
401  TCanvas* canvas8 = new TCanvas("canvas8", "canvas8", 800, 400);
402  canvas8->Divide(2,1);
403  canvas8->cd(1);
404  gPad->SetLeftMargin(0.16);
405  gPad->SetRightMargin(0.16);
406  canvas8->cd(1);
407  hists_ang_acceptance_itnorm->second->GetXaxis()->SetTitle("dX/dZ x 10^{3}");
408  hists_ang_acceptance_itnorm->second->GetYaxis()->SetTitle("dY/dZ x 10^{3}");
409  hists_ang_acceptance_itnorm->second->Draw("COLZ");
410  canvas8->cd(2);
411  gPad->SetLeftMargin(0.16);
412  gPad->SetRightMargin(0.16);
413 
414  // check the acceptance of the beam pipe as a function of z
415  // normalize to the plane -1
416  for (hists_ang_acceptance_it = hists_ang_acceptance.begin();
417  hists_ang_acceptance_it != hists_ang_acceptance.end();
418  hists_ang_acceptance_it++){\
419  if (hists_ang_acceptance_it == hists_ang_acceptance_itnorm){
420  continue;
421  }
422  canvas8->cd(2);
423  hists_ang_acceptance_it->second->Divide(hists_ang_acceptance_itnorm->second);
424  hists_ang_acceptance_it->second->SetMaximum(1.);
425  hists_ang_acceptance_it->second->GetXaxis()->SetTitle("dX/dZ x 10^{3}");
426  hists_ang_acceptance_it->second->GetYaxis()->SetTitle("dY/dZ x 10^{3}");
427  hists_ang_acceptance_it->second->Draw("COLZ");
428  canvas8->Update();
429  canvas8->Print("beampipe_results.ps(");
430  canvas8->Print("beampipe_acceptance.gif+");
431  }
432  canvas8->Clear();
433  canvas8->Print("beampipe_results.ps)");
434 
435 
436  // convert to pdf
437  cout << "converting to pdf " << system("ps2pdf beampipe_results.ps")
438  << endl;
439  cout << "deleting ps file " << system("rm beampipe_results.ps") << endl;
440 
441  output_histfile->Write();
442 
443  output_histfile->Close();
444  return 0;
445 }
nplanes
const int nplanes(120)
function_beampipe
Double_t function_beampipe(Double_t *x, Double_t *par)
Definition: geo/Check_particle_path.C:467
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
nlines
Int_t nlines
Definition: crosstag.C:23
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
DrawProgressBar
void DrawProgressBar(int len, double percent)
Definition: geo/Check_particle_path.C:447
startEvent
Int_t startEvent
Definition: emc/run_linkCreation.C:7
storePath
TString storePath
Definition: runGeaneTestSimBox.cxx:50
PndSdsMCPoint::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsMCPoint.h:89
PndMCTrack::IsGeneratorCreated
Bool_t IsGeneratorCreated(void) const
Definition: PndMCTrack.h:84
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
DigiFile
TString DigiFile
Definition: emc/run_linkCreation.C:14
nHits
int nHits
Definition: RiemannTest.C:16
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
parFile
TString parFile
Definition: hit_dirc.C:14
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
z
Double_t z
Definition: createRootFscGeometryFile.C:146
out
TFile * out
Definition: reco_muo.C:20
parInput1
FairParRootFileIo * parInput1
Definition: hit_dirc.C:67
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
if
if(fWindowIsBox)
Definition: createRootGeoFile.C:134
void DrawProgressBar ( int  len,
double  percent 
)

Definition at line 447 of file geo/Check_particle_path.C.

References i, and last_percent().

Referenced by Check_particle_path().

447  {
448  if ((int) (last_percent * 100) == (int) (percent * 100))
449  return;
450  //cout << " drawing " << endl;
451  cout << "\x1B[2K"; // Erase the entire current line.
452  cout << "\x1B[0E"; // Move to the beginning of the current line.
453  string progress;
454  for (int i = 0; i < len; ++i) {
455  if (i < static_cast<int> (len * percent)) {
456  progress += "=";
457  } else {
458  progress += " ";
459  }
460  }
461  cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
462  << "%";
463  flush( cout); // Required.
464  last_percent = percent;
465 }
last_percent
double last_percent(-1.)
i
Int_t i
Definition: run_full.C:25
Double_t function_beampipe ( Double_t *  x,
Double_t *  par 
)

Definition at line 467 of file geo/Check_particle_path.C.

References cos(), and sin().

Referenced by Check_particle_path().

468 {
469  double pos_z = x[0];
470  double result = 0.;
471  double end_seg_upstream = par[0]; // where bending starts with
472  double r_bend = par[1]; // a radius
473  double phi_bend = par[2]; // and the angle of the circle path
474  // the rest is fully determined
475  double end_seg_bend = end_seg_upstream+sin(phi_bend)*r_bend;
476  // the first straight part
477  if (pos_z < end_seg_upstream){
478  result = 0.;
479  }
480  // bending part is a part of a circle
481  if (end_seg_upstream <= pos_z && pos_z < end_seg_bend){
482  result = r_bend - cos(atan((pos_z-end_seg_upstream)/r_bend))*r_bend;
483  }
484  // straight part behind the dipole is a tangent to the circle
485  if (end_seg_bend <= pos_z){
486  result = (pos_z - end_seg_upstream - tan(phi_bend/2.)*r_bend)*tan(phi_bend);
487  }
488  return result;
489 }
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
par
Double_t par[3]
Definition: anaLmdCluster.C:146
x
Double_t x
Definition: createRootFscGeometryFile.C:159
double last_percent ( -  1.)

Referenced by DrawProgressBar().

int main ( void  )

Definition at line 491 of file geo/Check_particle_path.C.

References Check_particle_path().

491  {
492  //TApplication myapp("myapp", 0, 0);
493  Check_particle_path();
494  //myapp.Run();
495  return 0;
496 }
Check_particle_path
int Check_particle_path()
Definition: geo/Check_particle_path.C:49
const int nplanes ( 120  )

Referenced by Check_particle_path(), PndLmdPerformanceTask::Exec(), PndLmdPerformanceTask::Finish(), PndLmdPerformanceTask::Init(), reco_LMD(), and scattered_particles().

Generated on Sun Aug 18 2019 03:02:23 for FairRoot/PandaRoot by   doxygen 1.8.8