ROOTPWA
primaryVertexGen.cc
Go to the documentation of this file.
1 /*
2  * author: Prometeusz (Promme) jasinski
3  * jasinski@kph.uni-mainz.de Promme@web.de
4  *
5  * This file contains some methods simulating the Primary Vertex distribution
6  * in the target cell as well as the incoming beam properties
7  *
8  * or use the methods here to adapt it in your event generator
9  *
10  * (2010.03.03)
11  * - implementing 2008 h- beam simulation
12  *
13  * (2010.06.16)
14  * - moved to rootpwa svn repository and changed to class structure
15  * - changed to a class
16  *
17  */
18 
19 #include <iostream>
20 
21 #include "TH1.h"
22 #include "TH2.h"
23 #include "TCanvas.h"
24 #include "TFile.h"
25 #include "TStyle.h"
26 #include "TColor.h"
27 #include "TROOT.h"
28 #include "TVector3.h"
29 #include "TLorentzVector.h"
30 #include "TRandom.h"
31 
32 #include "primaryVertexGen.h"
33 
34 
35 using namespace std;
36 using namespace rpwa;
37 
38 primaryVertexGen::primaryVertexGen(string histfilename,
39  double beam_part_mass,
40  double mean_beam_energy,
41  double mean_beam_energy_spread)
42  : _histogramfile(NULL),
43  //static const string histogramfilename("properties_2008/primary_vertex_properties.root");
44  _hist_angles_vert_mean(NULL),
45  _hist_angles_horiz_mean(NULL),
46  _hist_angles_vert_sigma(NULL),
47  _hist_angles_horiz_sigma(NULL),
48  _hist_vertex_distr_xy(NULL),
49  _hist_vertex_distr_z(NULL),
50  _beam_part_mass(beam_part_mass),
51  _beam_energy_mean(mean_beam_energy),
52  _beam_energy_sigma(mean_beam_energy_spread)
53 {
54  _histograms_loaded = loadHistograms(histfilename, true);
55 }
56 
57 
58 primaryVertexGen::~primaryVertexGen() {
59  if(_histogramfile) {
60  _histogramfile->Close();
61  }
62  delete _histogramfile;
63  delete _hist_angles_vert_mean;
64  delete _hist_angles_horiz_mean;
65  delete _hist_angles_vert_sigma;
66  delete _hist_angles_horiz_sigma;
67  delete _hist_vertex_distr_xy;
68  delete _hist_vertex_distr_z;
69 }
70 
71 
72 bool primaryVertexGen::loadHistograms(string filename, bool plot) {
73  bool result(false);
74  if(plot) {
75  gROOT->SetStyle("Plain");
76  gStyle->SetPalette(1);
77  gesPalette();
78  }
79  TFile* histogramfile = new TFile(filename.c_str());
80  if(histogramfile->IsZombie()) {
81  cout << " Error: Could not read the given file containing histograms! " << endl;
82  return result;
83  }
84  _hist_angles_vert_mean = (TH2*)gDirectory->Get("hist_angles_vert_mean");
85  _hist_angles_horiz_mean = (TH2*)gDirectory->Get("hist_angles_horiz_mean");
86  _hist_angles_vert_sigma = (TH2*)gDirectory->Get("hist_angles_vert_sigma");
87  _hist_angles_horiz_sigma = (TH2*)gDirectory->Get("hist_angles_horiz_sigma");
88  _hist_vertex_distr_xy = (TH2*)gDirectory->Get("hist_vertex_distr_xy");
89  _hist_vertex_distr_z = (TH1*)gDirectory->Get("hist_vertex_distr_z");
90  // check if it worked
91  if(!_hist_angles_vert_mean ||
92  !_hist_angles_horiz_mean ||
93  !_hist_angles_vert_sigma ||
94  !_hist_angles_horiz_sigma ||
95  !_hist_vertex_distr_xy ||
96  !_hist_vertex_distr_z)
97  {
98  cout << " Error: histograms not found!" << endl;
99  return result;
100  }
101  /*
102  for (int x = 0; x < n_steps_x; x++){
103  for (int y = 0; y < n_steps_y; y++){
104  hist_angles_vert_mean->SetBinContent(x+1, y+1, mean_vert[x][y]);
105  hist_angles_vert_sigma->SetBinContent(x+1, y+1, sigma_vert[x][y]);
106  hist_angles_horiz_mean->SetBinContent(x+1, y+1, mean_horiz[x][y]);
107  hist_angles_horiz_sigma->SetBinContent(x+1, y+1, sigma_horiz[x][y]);
108  }
109  }
110  */
111  // draw the output
112  if(plot) {
113  TCanvas *analyze_beam_properties_output = new TCanvas("analyze_beam_properties_output", "analyze_beam_properties_output method output", 600, 900);
114  analyze_beam_properties_output->Divide(2,3);
115  analyze_beam_properties_output->cd(1);
116  _hist_angles_vert_mean->Draw("COLZ");
117  gPad->Update();
118  analyze_beam_properties_output->cd(2);
119  _hist_angles_vert_sigma->Draw("COLZ");
120  gPad->Update();
121  analyze_beam_properties_output->cd(3);
122  _hist_angles_horiz_mean->Draw("COLZ");
123  gPad->Update();
124  analyze_beam_properties_output->cd(4);
125  _hist_angles_horiz_sigma->Draw("COLZ");
126  gPad->Update();
127  analyze_beam_properties_output->cd(5);
128  gPad->SetLogz();
129  _hist_vertex_distr_xy->Draw("COLZ");
130  gPad->Update();
131  analyze_beam_properties_output->cd(6);
132  _hist_vertex_distr_z->Draw("");
133  gPad->Update();
134  analyze_beam_properties_output->Print("Fitted_beam_property_distributions.pdf");
135  }
136  result = true;
137  return result;
138 }
139 
140 
141 bool primaryVertexGen::check() {
142  return _histograms_loaded;
143 }
144 
145 
146 TVector3& primaryVertexGen::getVertex(const float cutR,
147  const float cutZ_low,
148  const float cutZ_high)
149 {
150  if(!_histograms_loaded) {
151  TVector3* result = new TVector3();
152  return *result;
153  }
154  double x;
155  double y;
156  double z;
157  bool tryagain(true);
158  while(tryagain) {
159  _hist_vertex_distr_xy->GetRandom2(x, y);
160  z = _hist_vertex_distr_z->GetRandom();
161  if((sqrt(x*x+y*y) < cutR) && (z > cutZ_low) && (z < cutZ_high)) {
162  tryagain = false;
163  }
164  }
165  TVector3* result = new TVector3(x,y,z);
166  return *result;
167 }
168 
169 
170 TVector3& primaryVertexGen::getBeamDir(const TVector3 vertex) {
171  TVector3* result = new TVector3(0.,0.,1.);
172  if(!_histograms_loaded) {
173  return *result;
174  }
175  // check whever we are in the valid ranges of the histograms
176  double sigma_horiz= _hist_angles_horiz_sigma->GetBinContent(
177  _hist_angles_horiz_sigma->FindBin(vertex.X(), vertex.Y()));
178  double sigma_vert = _hist_angles_vert_sigma->GetBinContent(
179  _hist_angles_vert_sigma->FindBin(vertex.X(), vertex.Y()));
180  if((sigma_horiz == 0.) || (sigma_vert == 0.)) {
181  //cout << " Warning in Get_Beam: out of histogram range " << endl;
182  result->SetZ(0.);
183  return *result;
184  }
185  // if check passed retrieve interpolated values
186  double mean_horiz = _hist_angles_horiz_mean->Interpolate(vertex.X(), vertex.Y());
187  double mean_vert = _hist_angles_vert_mean->Interpolate(vertex.X(), vertex.Y());
188  sigma_horiz = _hist_angles_horiz_sigma->Interpolate(vertex.X(), vertex.Y());
189  sigma_vert = _hist_angles_vert_sigma->Interpolate(vertex.X(), vertex.Y());
190  // now we can retrieve the direction randomly
191  double angle_horiz = gRandom->Gaus(mean_horiz, sigma_horiz);
192  double angle_vert = gRandom->Gaus(mean_vert , sigma_vert);
193  result->RotateX(-angle_vert);
194  result->RotateY(angle_horiz);
195  return *result;
196 }
197 
198 TLorentzVector& primaryVertexGen::getBeamPart(const TVector3 beam_dir) {
199  double energy = gRandom->Gaus(_beam_energy_mean, _beam_energy_sigma);
200  TVector3 _beam_dir(beam_dir);
201  // m² = E² - p² -> |p| = sqrt(E²-m²)
202  _beam_dir.SetMag(sqrt(energy*energy-_beam_part_mass*_beam_part_mass));
203  TLorentzVector* result = new TLorentzVector(_beam_dir.X(), _beam_dir.Y(), _beam_dir.Z(), energy);
204  return *result;
205 }
206 
207 void primaryVertexGen::gesPalette(int i) {
208  TPad foo; // never remove this line :-)))
209  if(i == 0) {
210  const Int_t NRGBs = 5;
211  const Int_t NCont = 255;
212  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
213  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
214  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
215  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
216  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
217  gStyle->SetNumberContours(NCont);
218  }
219 
220  if(i == 1) {
221  const UInt_t Number = 3;
222  Double_t Red[Number] = { 1.00, 0.00, 0.00};
223  Double_t Green[Number] = { 0.00, 1.00, 0.00};
224  Double_t Blue[Number] = { 1.00, 0.00, 1.00};
225  Double_t Length[Number] = { 0.00, 0.50, 1.00 };
226  Int_t nb=255;
227  TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
228  gStyle->SetNumberContours(nb);
229  }
230 }
231 
232 void Beam_simulation() {
233  primaryVertexGen primaryVertexGen;
234  if(!primaryVertexGen.check()) {
235  return;
236  }
237  // some cross check histograms
238  //const int n_steps_x = 21;
239  //const int n_steps_y = 21;
240  //const float cell_size_x = 4.;
241  //const float cell_size_y = 4.;
242 
243  TH2F *hist_vertex_distr_xy_sim = new TH2F("hist_vertex_distr_xy_sim", "simulated vertex X Y distribution", 200, -2., 2., 200, -2., 2.);
244  TH1F *hist_vertex_distr_z_sim = new TH1F("hist_vertex_distr_z_sim","simulated vertex Z distribution",1000, -100, 0);
245  TH1F *hist_angle_distr_horiz = new TH1F("hist_angle_distr_horiz","horizontal angle distribution",1000, -0.001, 0.001);
246  TH1F *hist_angle_distr_vert = new TH1F("hist_angle_distr_vert ","vertical angle distribution",1000, -0.001, 0.001);
247  TH1F *hist_energy_distr_sim = new TH1F("hist_energy_distr_sim", "simulated energy distribution",1000, 170, 210);
248  TH2F* hist_beamtrack_horiz_plane = new TH2F("hist_beamtrack_horiz_plane", "beamtrack in horizontal plane", 1000, 0, 6000, 100, -5, 5);
249  TH2F* hist_beamtrack_vert_plane = new TH2F("hist_beamtrack_vert_plane", "beamtrack in vertical plane", 1000, 0, 6000, 100, -5, 5);
250  TH1F* hist_horiz_beamdivergence_at_CEDAR = new TH1F("hist_horiz_beamdivergence_at_CEDAR", "horizontal beam divergence at CEDAR region", 1000, -400e-6, 400e-6);
251  TH1F* hist_vert_beamdivergence_at_CEDAR = new TH1F("hist_vert_beamdivergence_at_CEDAR", "vertical beam divergence at CEDAR region", 1000, -400e-6, 400e-6);
252 
253  int counter = 0;
254  TCanvas* canvas = new TCanvas("canvas", "simulated distributions", 900, 600);
255  canvas->Divide(3,2);
256  TCanvas* canvas2 = new TCanvas("canvas2", "simulated tracks", 600, 600);
257  canvas2->Divide(1,2);
258  TCanvas* canvas3 = new TCanvas("canvas3", "simulated distributions at CEDAR region", 600, 400);
259  canvas3->Divide(2,1);
260  for(int i = 0; i < 1000000; i++) {
261  TVector3 vertex = primaryVertexGen.getVertex();
262  TVector3 beam_dir = primaryVertexGen.getBeamDir(vertex);
263  if(beam_dir.Mag() == 0) {
264  //cout << " skipping " << endl;
265  continue;
266  }
267  hist_vertex_distr_xy_sim->Fill(vertex.X(), vertex.Y());
268  hist_vertex_distr_z_sim->Fill(vertex.Z());
269  TLorentzVector beam_part = primaryVertexGen.getBeamPart(beam_dir);
270  //cout << beam_part.M() << endl;
271  double azi = beam_part.Px()/beam_part.Pz();
272  double dip = beam_part.Py()/beam_part.Pz();
273  hist_angle_distr_horiz->Fill(azi);
274  hist_angle_distr_vert->Fill(dip);
275  hist_energy_distr_sim->Fill(beam_part.E());
276 
277  // check the beam spot 30m downstream by extrapolation from 0 to 6000 cm
278  for(int i = 0; i < 1000; i++) {
279  // get the distance between pos z of the vertex to the point
280  // to be extrapolated
281  float pos_z = i*6;
282  float extrapol_dist_z = pos_z - vertex.Z();
283  // now we can calculate the factor for the pointing vector for extrapolation
284  float extrapol_fac = extrapol_dist_z/beam_dir.Z();
285  // this factor has to be applied in all 3 dimensions
286  hist_beamtrack_horiz_plane->Fill(pos_z, vertex.X()+beam_dir.X()*extrapol_fac);
287  hist_beamtrack_vert_plane ->Fill(pos_z, vertex.Y()+beam_dir.Y()*extrapol_fac);
288  }
289  // have a look on the properties at the CEDAR region
290  // transport the vertex position 30m downstream
291  float pos_z = 3000;
292  float extrapol_dist_z = pos_z - vertex.Z();
293  float extrapol_fac = extrapol_dist_z/beam_dir.Z();
294  vertex.SetXYZ(vertex.X()+beam_dir.X()*extrapol_fac,
295  vertex.Y()+beam_dir.Y()*extrapol_fac,
296  pos_z);
297 
298  // transport the particle track back to the CEDAR postion upstream
299  // by using the transportation matrix output given by Lau 26.01.2009
300  /*
301  * 1COMPASS HADRON OPTICS 200 GEV V.2007 TRANSPORT RUN26/01/09
302  0POSITION TYPE STRENGTH * H O R I Z O N T A L * V E R T I C A L * D I S P E R S I O N
303  METERS T*M,T/M*M * R11 R12 R21 R22 * R33 R34 R43 R44 * R16 R26 R36 R46
304  T/M**2*M * MM/MM MM/MR MR/MM MR/MR * MM/MM MM/MR MR/MM MR/MR * MM/PC MR/PC MM/PC MR/PC
305  ************************************************************************************************************************************
306  30.000 3 * 1.000 30.000 0.000 1.000 * 1.000 30.000 0.000 1.000 * 0.000 0.000 0.000 0.000
307  69.396 3 * 0.670 47.160 -0.021 0.000 * 0.924 73.374 -0.014 0.000 * 0.000 0.000 0.772 0.022
308  75.713 3 CED2 * 0.536 47.160 -0.021 0.000 * 0.838 73.374 -0.014 0.000 * 0.000 0.000 0.908 0.022
309  75.913 3 * 0.531 47.160 -0.021 0.000 * 0.835 73.374 -0.014 0.000 * 0.000 0.000 0.912 0.022
310  82.230 3 CED1 * 0.398 47.160 -0.021 0.000 * 0.749 73.374 -0.014 0.000 * 0.000 0.000 1.048 0.022
311  ^very parallel beam ^very parallel beam :)
312  */
313 
314  // compass is measuring in cm!
315  // x is horizontal plane
316  // y is vertical plane
317  // the azimuth lies in x
318  // the dip in y
319  // at 82.230 m upstream the CEDAR in Lau's table
320  float x_displacement_upstream = vertex.X()*(0.398) + azi*(47.160*100); // mm/mrad -> cm/rad = * 1000/10 = * 100
321  float y_displacement_upstream = vertex.Y()*(0.749) + dip*(73.374*100);
322  // at 69.396 m downstream the CEDAR in Lau's table
323  float x_displacement_downstream = vertex.X()*(0.670) + azi*(47.160*100);
324  float y_displacement_downstream = vertex.Y()*(0.924) + dip*(73.374*100);
325  // and calculate the beam divergence
326  float divergence_horiz = (x_displacement_downstream-x_displacement_upstream) / ((82.230-69.396) * 100);
327  float divergence_vert = (y_displacement_downstream-y_displacement_upstream) / ((82.230-69.396) * 100);
328  hist_horiz_beamdivergence_at_CEDAR->Fill(divergence_horiz);
329  hist_vert_beamdivergence_at_CEDAR->Fill(divergence_vert);
330 
331  ++counter;
332  if((counter % 10000) == 0) {
333  canvas->cd(1);
334  gPad->Clear();
335  gPad->SetLogz();
336  hist_vertex_distr_xy_sim->Draw("COLZ");
337  gPad->Update();
338  canvas->cd(2);
339  gPad->Clear();
340  hist_vertex_distr_z_sim->Draw("");
341  gPad->Update();
342  canvas->cd(3);
343  gPad->Clear();
344  hist_angle_distr_horiz->Draw("");
345  gPad->Update();
346  canvas->cd(4);
347  gPad->Clear();
348  hist_angle_distr_vert->Draw("");
349  gPad->Update();
350  canvas->cd(5);
351  gPad->Clear();
352  hist_energy_distr_sim->Draw("");
353  gPad->Update();
354  canvas2->cd(1);
355  gPad->Clear();
356  hist_beamtrack_horiz_plane->Draw("COLZ");
357  gPad->Update();
358  canvas2->cd(2);
359  gPad->Clear();
360  hist_beamtrack_vert_plane->Draw("COLZ");
361  gPad->Update();
362  canvas3->cd(1);
363  gPad->Clear();
364  hist_horiz_beamdivergence_at_CEDAR->Draw();
365  gPad->Update();
366  canvas3->cd(2);
367  gPad->Clear();
368  hist_vert_beamdivergence_at_CEDAR->Draw();
369  gPad->Update();
370  }
371  }
372  canvas2->Print("Simulated_tracks.pdf");
373  canvas->Print("Simulated_distributions.pdf");
374  canvas3->Print("Simulated_distributions_at_CEDAR.pdf");
375 }
376 
Generated by   doxygen 1.8.1.2