FairRoot/PandaRoot
Functions | Variables
Promme/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 <TGeoVolume.h>
#include <TGraph.h>
#include <TF1.h>
#include <algorithm>
#include <fstream>
#include <TStyle.h>
#include <TColor.h>
#include <FairGeoLoader.h>
#include <FairGeoInterface.h>
#include <FairGeoBuilder.h>
#include <FairGeoPcon.h>
#include <FairGeoMedia.h>
#include <TGeoManager.h>
#include <TImage.h>
#include <TLine.h>
#include <stdlib.h>
#include <TPolyLine.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)
 
TPolyLine * Get_PolyLine_circle (float x, float y, float r)
 
int Check_particle_path ()
 
int main ()
 

Variables

string geo_beampipe = "/home/jasinski/bin/pandaroot/macro/lmd/Promme/beampipe_201209.root"
 

Function Documentation

int Check_particle_path ( )

Definition at line 67 of file Promme/Check_particle_path.C.

References DigiFile, Double_t, DrawProgressBar(), fRun, function_beampipe(), Get_PolyLine_circle(), 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().

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

Definition at line 539 of file Promme/Check_particle_path.C.

References i, and last_percent().

Referenced by Check_particle_path().

539  {
540  if ((int) (last_percent * 100) == (int) (percent * 100))
541  return;
542  //cout << " drawing " << endl;
543  cout << "\x1B[2K"; // Erase the entire current line.
544  cout << "\x1B[0E"; // Move to the beginning of the current line.
545  string progress;
546  for (int i = 0; i < len; ++i) {
547  if (i < static_cast<int> (len * percent)) {
548  progress += "=";
549  } else {
550  progress += " ";
551  }
552  }
553  cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
554  << "%";
555  flush( cout); // Required.
556  last_percent = percent;
557 }
i
Int_t i
Definition: run_full.C:25
last_percent
double last_percent(-1.)
Double_t function_beampipe ( Double_t *  x,
Double_t *  par 
)

Definition at line 559 of file Promme/Check_particle_path.C.

References cos(), and sin().

Referenced by Check_particle_path().

560 {
561  double pos_z = x[0];
562  double result = 0.;
563  double end_seg_upstream = par[0]; // where bending starts with
564  double r_bend = par[1]; // a radius
565  double phi_bend = par[2]; // and the angle of the circle path
566  // the rest is fully determined
567  double end_seg_bend = end_seg_upstream+sin(phi_bend)*r_bend;
568  // the first straight part
569  if (pos_z < end_seg_upstream){
570  result = 0.;
571  }
572  // bending part is a part of a circle
573  if (end_seg_upstream <= pos_z && pos_z < end_seg_bend){
574  result = r_bend - cos(atan((pos_z-end_seg_upstream)/r_bend))*r_bend;
575  }
576  // straight part behind the dipole is a tangent to the circle
577  if (end_seg_bend <= pos_z){
578  result = (pos_z - end_seg_upstream - tan(phi_bend/2.)*r_bend)*tan(phi_bend);
579  }
580  return result;
581 }
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
TPolyLine * Get_PolyLine_circle ( float  x,
float  y,
float  r 
)

Definition at line 583 of file Promme/Check_particle_path.C.

References angle, cos(), Double_t, i, Pi, sin(), x, and y.

Referenced by Check_particle_path().

583  {
584  const int nsegments = 360;
585  Double_t _x[nsegments+1] = {0.,}; // +1 : Endpoint is Startpoint
586  Double_t _y[nsegments+1] = {0.,};
587  for (int i = 0; i <= nsegments; i++){
588  double angle = 2*TMath::Pi()*((double) i)/((double) nsegments);
589  _x[i] = r*sin(angle)+x;
590  _y[i] = r*cos(angle)+y;
591  }
592  TPolyLine *pline_circle = new TPolyLine(nsegments+1,_x,_y);
593  return pline_circle;
594 }
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
r
double r
Definition: RiemannTest.C:14
i
Int_t i
Definition: run_full.C:25
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
x
Double_t x
Definition: createRootFscGeometryFile.C:159
y
Double_t y
Definition: createRootFscGeometryFile.C:159
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
double last_percent ( -  1.)

Referenced by DrawProgressBar().

int main ( void  )

Definition at line 596 of file Promme/Check_particle_path.C.

References Check_particle_path().

596  {
597  //TApplication myapp("myapp", 0, 0);
598  Check_particle_path();
599  //myapp.Run();
600  return 0;
601 }
Check_particle_path
int Check_particle_path()
Definition: Promme/Check_particle_path.C:67
const int nplanes ( 120  )

Referenced by Check_particle_path().

Variable Documentation

string geo_beampipe = "/home/jasinski/bin/pandaroot/macro/lmd/Promme/beampipe_201209.root"

Definition at line 65 of file Promme/Check_particle_path.C.

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