ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
circlegen.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file circlegen.cpp
1 // a. dion - author of code
2 // j. nagle - added more documentation (07/29/2011)
3 //
4 // routine generates a ROOT output with M hits in cylindrical detector layers
5 // (with radii and number of layers defined below in "vector<double> radii")
6 // from N (as user input) tracks following a helical pattern from a uniform solenoidal
7 // magnetic field.
8 //
9 // now with autogen and make, puts executable in install/bin
10 //
11 // recent addition for smearing of position resolution of hits (as test)
12 // see boolean "resolution_smear" below
13 //
14 // j.nagle - added ability to include noise hits
15 //
16 
17 #include <math.h>
18 #include "TFile.h"
19 #include "TTree.h"
20 #include "TRandom3.h"
21 #include <sstream>
22 #include <iostream>
23 #include <vector>
24 
25 using namespace std;
26 
27 bool intersect_circles(bool hel, double startx, double starty, double rad_det, double rad_trk, double cx, double cy, double& x, double& y)
28 {
29  double d2 = (cx*cx + cy*cy);
30  double d = sqrt(d2);
31  if(d > (rad_det + rad_trk))
32  {
33  return false;
34  }
35  if(d < fabs(rad_det - rad_trk))
36  {
37  return false;
38  }
39 
40  double r2 = rad_trk*rad_trk;
41 
42  double d_inv = 1./d;
43  double R2 = rad_det*rad_det;
44  double a = 0.5*(R2 - r2 + d2)*d_inv;
45  double h = a*d_inv;
46  double P2x = cx*h;
47  double P2y = cy*h;
48  h = sqrt(R2 - a*a);
49 
50  double ux = -cy*d_inv;
51  double uy = cx*d_inv;
52  double P3x1 = P2x + ux*h;
53  double P3y1 = P2y + uy*h;
54  ux = -ux;
55  uy = -uy;
56  double P3x2 = P2x + ux*h;
57  double P3y2 = P2y + uy*h;
58 
59  double d1_2 = (startx - P3x1)*(startx - P3x1) + (starty - P3y1)*(starty - P3y1);
60  double d2_2 = (startx - P3x2)*(startx - P3x2) + (starty - P3y2)*(starty - P3y2);
61 
62  if(d1_2 < d2_2)
63  {
64  x = P3x1;
65  y = P3y1;
66  }
67  else
68  {
69  x = P3x2;
70  y = P3y2;
71  }
72 
73  return true;
74 }
75 
76 
77 void print_usage()
78 {
79  cout<<"usage: circle_gen <# of events> <max # of tracks per event> <output file name>"<<endl;
80 }
81 
82 
83 int main(int argc, char** argv)
84 {
85  stringstream ss;
86  unsigned int nmaxtracks=0;
87  unsigned int nevents=0;
88 
89  bool all_same_track_number = true;
90 
91  // inclusion of noise hits by these fake tracks
92  bool faketracks = true;
93  int nfaketracks = 250;
94 
95  bool resolution_smear = true;
96  double smear_x = (70.0e-4/sqrt(12.))/sqrt(2.0); // 70 microns as example
97  double smear_y = (70.0e-4/sqrt(12.))/sqrt(2.0); // 70 microns as example
98  double smear_z = 425.0e-4/sqrt(12.); // 425 microns as example
99 
100 
101 
102  // how many layers to include....
103  //===============================
104 // sPHENIX style with 6 layers
105  int nlayers = 6;
106  vector<double> radii;
107  radii.assign(nlayers,0.);
108  radii[0]=2.5;
109  radii[1]=5.0;
110  radii[2]=10.0;
111  radii[3]=14.0;
112  radii[4]=40.0;
113  radii[5]=60.0;
114 
115  vector<double> smear_x_layer;smear_x_layer.assign(nlayers,0);
116  vector<double> smear_y_layer;smear_y_layer.assign(nlayers,0);
117  vector<double> smear_z_layer;smear_z_layer.assign(nlayers,0);
118  smear_x_layer[0] = (50.0e-4/sqrt(12.))/sqrt(2.0);
119  smear_y_layer[0] = (50.0e-4/sqrt(12.))/sqrt(2.0);
120  smear_z_layer[0] = (425.0e-4/sqrt(12.));
121  smear_x_layer[1] = (50.0e-4/sqrt(12.))/sqrt(2.0);
122  smear_y_layer[1] = (50.0e-4/sqrt(12.))/sqrt(2.0);
123  smear_z_layer[1] = (425.0e-4/sqrt(12.));
124  smear_x_layer[2] = (80.0e-4/sqrt(12.))/sqrt(2.0);
125  smear_y_layer[2] = (80.0e-4/sqrt(12.))/sqrt(2.0);
126  smear_z_layer[2] = (1000.0e-4/sqrt(12.));
127  smear_x_layer[3] = (80.0e-4/sqrt(12.))/sqrt(2.0);
128  smear_y_layer[3] = (80.0e-4/sqrt(12.))/sqrt(2.0);
129  smear_z_layer[3] = (1000.0e-4/sqrt(12.));
130  smear_x_layer[4] = (100.0e-4/sqrt(12.))/sqrt(2.0);
131  smear_y_layer[4] = (100.0e-4/sqrt(12.))/sqrt(2.0);
132  smear_z_layer[4] = (10000.0e-4/sqrt(12.));
133  smear_x_layer[5] = (100.0e-4/sqrt(12.))/sqrt(2.0);
134  smear_y_layer[5] = (100.0e-4/sqrt(12.))/sqrt(2.0);
135  smear_z_layer[5] = (10000.0e-4/sqrt(12.));
136 
137 
138  // current SVX
139 // int nlayers = 4;
140 // vector<double> radii;
141 // radii.assign(nlayers,0.);
142 // radii[0]=2.5;
143 // radii[1]=5.0;
144 // radii[2]=10.0;
145 // radii[3]=14.0;
146  cout << "Circle Gen Code : running with layers = " << nlayers << " and resolution smearing = " << resolution_smear << endl;
147 
148  double kappa_min=0.001;
149  double kappa_max=0.03;
150  double d_min=-0.2;
151  double d_max=0.2;
152  double phi_min=0.;
153  double phi_max=6.28318530717958623;
154  double dzdl_min = -0.5;
155  double dzdl_max = 0.5;
156  double z0_min=-5.;
157  double z0_max=5.;
158 
159  ss.clear();ss.str("");
160  ss<<argv[1];
161  ss>>nevents;
162  if(nevents==0){print_usage();return -1;}
163  ss.clear();ss.str("");
164  ss<<argv[2];
165  ss>>nmaxtracks;
166  if(nmaxtracks==0){print_usage();return -1;}
167  ss.clear();ss.str("");
168  ss<<argv[3];
169  if(ss.str() == ""){print_usage();return -1;}
170  string filename = ss.str();
171 
172  unsigned int track=0;
173  unsigned int nhits=0;
174  double x_hits[12];
175  double y_hits[12];
176  double z_hits[12];
177  int layer[12];
178  double trk_kappa=0.;
179  double trk_d=0.;
180  double trk_phi=0.;
181  double trk_dzdl=0.;
182  double trk_z0=0.;
183  unsigned char charge=0;
184  TFile* ofile = new TFile(filename.c_str(), "recreate");
185  TTree* etree = new TTree("events", "event list");
186  TTree* otree = new TTree("tracks", "tracks in a single event");
187  otree->Branch("track", &track, "track/i");
188  otree->Branch("nhits", &nhits, "nhits/i");
189  otree->Branch("x_hits", x_hits, "x_hits[nhits]/D");
190  otree->Branch("y_hits", y_hits, "y_hits[nhits]/D");
191  otree->Branch("z_hits", z_hits, "z_hits[nhits]/D");
192  otree->Branch("layer", layer, "layer[nhits]/i");
193  otree->Branch("trk_kappa", &trk_kappa, "trk_kappa/D");
194  otree->Branch("trk_d", &trk_d, "trk_d/D");
195  otree->Branch("trk_phi", &trk_phi, "trk_phi/D");
196  otree->Branch("trk_dzdl", &trk_dzdl, "trk_dzdl/D");
197  otree->Branch("trk_z0", &trk_z0, "trk_z0/D");
198  otree->Branch("charge", &charge, "charge/b");
199  etree->Branch("tracklist", &otree);
200 
201  double secondary_proportion=0.1;
202 
203 
204  TRandom3 rand;
205  double kappa=0.;
206  double d=0.;
207  double phi=0.;
208  double dzdl=0.;
209  double z0=0.;
210 
211  for(unsigned int ev=0;ev<nevents;ev++)
212  {
213  otree->Reset();
214  unsigned int ntracks = (unsigned int)(rand.Uniform(1., nmaxtracks));
215  if (all_same_track_number) ntracks = (unsigned int) nmaxtracks;
216 
217  for(unsigned int i=0;i<ntracks;i++)
218  {
219  kappa = rand.Uniform(kappa_min, kappa_max);
220  d=0.;
221  phi = rand.Uniform(phi_min, phi_max);
222  dzdl = rand.Uniform(dzdl_min, dzdl_max);
223  z0 = 0.;
224  bool hel=false;
225  double temp = rand.Uniform(-1., 1.);
226  if(temp > 0.){hel = true;}
227 
228  double ux=0.;
229  double uy=0.;
230  double tanphi=0.;
231  if((fabs(phi) < 7.85398163397448279e-01) || (fabs(phi - 3.1415926535897931) < 7.85398163397448279e-01))
232  {
233  tanphi = tan(phi);
234  ux = sqrt(1./(1. + tanphi*tanphi));
235  if(phi > 1.57079632679489656 && phi < 4.71238898038468967){ux=-ux;}
236  uy = ux*tanphi;
237  }
238  else
239  {
240  tanphi = tan(1.57079632679489656 - phi);//cotangent
241  uy = sqrt(1./(tanphi*tanphi + 1.));
242  if(phi > 3.1415926535897931){uy=-uy;}
243  ux = uy*tanphi;
244  }
245  double dx = ux*d;
246  double dy = uy*d;
247  double dispx = ux/kappa;
248  double dispy = uy/kappa;
249  double cx = dx + dispx;
250  double cy = dy + dispy;
251 
252  double x=0.;
253  double y=0.;
254  double z=0.;
255 
256  double startx=dx;
257  double starty=dy;
258  double startz=z0;
259  nhits=0;
260  for(unsigned int l=0;l<radii.size();l++)
261  {
262  if(intersect_circles(hel, startx, starty, radii[l], 1./kappa, cx, cy, x, y) == true)
263  {
264  layer[nhits]=l;
265  // could smear these by some resolution factor here...
266  x_hits[nhits]=x;
267  y_hits[nhits]=y;
268  if (resolution_smear) {
269  x_hits[nhits] += rand.Gaus(0.0,smear_x_layer[l]); // in cm units
270  y_hits[nhits] += rand.Gaus(0.0,smear_y_layer[l]); // in cm units
271  }
272 
273  double k = kappa;
274  double D = sqrt((startx-x)*(startx-x) + (starty-y)*(starty-y));
275 
276  double s=0.;
277  if(0.5*k*D > 0.1)
278  {
279  double v = 0.5*k*D;
280  if(v >= 0.999999){v=0.999999;}
281  s = 2.*asin(v)/k;
282  }
283  else
284  {
285  double temp1 = k*D*0.5;temp1*=temp1;
286  double temp2 = D*0.5;
287  s += 2.*temp2;
288  temp2*=temp1;
289  s += temp2/3.;
290  temp2*=temp1;
291  s += (3./20.)*temp2;
292  temp2*=temp1;
293  s += (5./56.)*temp2;
294  }
295  double dz = sqrt(s*s*dzdl*dzdl/(1. - dzdl*dzdl));
296  if(dzdl>0.){z = startz + dz;}
297  else{z = startz - dz;}
298 
299  z_hits[nhits]=z;
300  if (resolution_smear) {
301  z_hits[nhits] += rand.Gaus(0.0,smear_z_layer[l]); // in cm units
302  }
303 
304 
305  if((x_hits[nhits] == x_hits[nhits]) && (y_hits[nhits] == y_hits[nhits]) && (z_hits[nhits] == z_hits[nhits]))
306  {
307  nhits++;
308  startx = x;
309  starty = y;
310  startz = z;
311  }
312  else
313  {
314  break;
315  }
316  }
317  }
318 
319  if(nhits>0)
320  {
321  double v1x = x_hits[0];
322  double v1y = y_hits[0];
323  double v2x = x_hits[1];
324  double v2y = y_hits[1];
325 
326  double diffx = v2x-v1x;
327  double diffy = v2y-v1y;
328 
329  double cross = diffx*cy - diffy*cx;
330 
331  int cross_sign = 1;
332  if(cross < 0.){cross_sign=-1;}
333 // cout<<hel<<" "<<cross_sign<<endl;
334 
335 
336  trk_kappa=kappa;
337  trk_d=d;
338  trk_phi=phi;
339  trk_dzdl=dzdl;
340  trk_z0=z0;
341  charge=0;
342  if(hel==true){charge=1;}
343  track=i;
344 // for(int nnh=0;nnh<nhits;++nnh)
345 // {
346 // cout<<x_hits[nnh]<<" "<<y_hits[nnh]<<" "<<z_hits[nnh]<<endl;
347 // }
348  otree->Fill();
349  }
350 
351  } // end loop over tracks
352 
353 
354  // how about noise tracks (i.e. made up of random hits from nlayers different tracks (so same spatial distribution)
355  // j.nagle - add switch for this features....
356  if (faketracks) {
357 
358  for(unsigned int j=0;j<nfaketracks;j++) {
359 
360 
361 
362  for(unsigned int i=0;i<nlayers;i++)
363  {
364  nhits=0;
365  kappa = rand.Uniform(kappa_min, kappa_max);
366  d=0.;
367  phi = rand.Uniform(phi_min, phi_max);
368  dzdl = rand.Uniform(dzdl_min, dzdl_max);
369  z0 = 0.;
370  bool hel=false;
371  double temp = rand.Uniform(-1., 1.);
372  if(temp > 0.){hel = true;}
373 
374  double ux=0.;
375  double uy=0.;
376  double tanphi=0.;
377  if((fabs(phi) < 7.85398163397448279e-01) || (fabs(phi - 3.1415926535897931) < 7.85398163397448279e-01))
378  {
379  tanphi = tan(phi);
380  ux = sqrt(1./(1. + tanphi*tanphi));
381  if(phi > 1.57079632679489656 && phi < 4.71238898038468967){ux=-ux;}
382  uy = ux*tanphi;
383  }
384  else
385  {
386  tanphi = tan(1.57079632679489656 - phi);//cotangent
387  uy = sqrt(1./(tanphi*tanphi + 1.));
388  if(phi > 3.1415926535897931){uy=-uy;}
389  ux = uy*tanphi;
390  }
391  double dx = ux*d;
392  double dy = uy*d;
393  double dispx = ux/kappa;
394  double dispy = uy/kappa;
395  double cx = dx + dispx;
396  double cy = dy + dispy;
397 
398  double x=0.;
399  double y=0.;
400  double z=0.;
401 
402  double startx=dx;
403  double starty=dy;
404  double startz=z0;
405 
406  // only check one layer for each of these tracks
407  for(unsigned int l=i;l<i+1;l++)
408  {
409  if(intersect_circles(hel, startx, starty, radii[l], 1./kappa, cx, cy, x, y) == true)
410  {
411  layer[nhits]=l;
412  // could smear these by some resolution factor here...
413  x_hits[nhits]=x;
414  y_hits[nhits]=y;
415  if (resolution_smear) {
416  x_hits[nhits] += rand.Gaus(0.0,smear_x); // in cm units
417  y_hits[nhits] += rand.Gaus(0.0,smear_y); // in cm units
418  }
419 
420  double k = kappa;
421  double D = sqrt((startx-x)*(startx-x) + (starty-y)*(starty-y));
422 
423  double s=0.;
424  if(0.5*k*D > 0.1)
425  {
426  s = 2.*asin(0.5*k*D)/k;
427  }
428  else
429  {
430  double temp1 = k*D*0.5;temp1*=temp1;
431  double temp2 = D*0.5;
432  s += 2.*temp2;
433  temp2*=temp1;
434  s += temp2/3.;
435  temp2*=temp1;
436  s += (3./20.)*temp2;
437  temp2*=temp1;
438  s += (5./56.)*temp2;
439  }
440  double dz = sqrt(s*s*dzdl*dzdl/(1. - dzdl*dzdl));
441  if(dzdl>0.){z = startz + dz;}
442  else{z = startz - dz;}
443 
444  z_hits[nhits]=z;
445  if (resolution_smear) {
446  z_hits[nhits] += rand.Gaus(0.0,smear_z); // in cm units
447  }
448 
449  nhits++;
450 
451  startx = x;
452  starty = y;
453  startz = z;
454  }
455  }
456 
457  if(nhits>0)
458  {
459  trk_kappa=kappa;
460  trk_d=d;
461  trk_phi=phi;
462  trk_dzdl=dzdl;
463  trk_z0=z0;
464  charge=0;
465  //if(hel==true){charge=1;}
466  charge=1; // fake track anyway so no real charge
467  track=otree->GetEntries();
468  otree->Fill();
469  }
470 
471 
472 
473  } // end loop over tracks leaving hits in each layer
474 
475 
476 
477  } // end loop over fake tracks
478 
479  } // end if faketracks
480 
481  /*
482  ntracks = (unsigned int)((secondary_proportion)*((double)ntracks));
483 
484  //now the secondaries
485  for(unsigned int i=0;i<ntracks;i++)
486  {
487  kappa = rand.Uniform(kappa_min, kappa_max);
488  d = rand.Uniform(d_min, d_max);
489  phi = rand.Uniform(phi_min, phi_max);
490  dzdl = rand.Uniform(dzdl_min, dzdl_max);
491  z0 = rand.Uniform(z0_min, z0_max);
492  bool hel=false;
493  double temp = rand.Uniform(-1., 1.);
494  if(temp > 0.){hel = true;}
495 
496  double ux=0.;
497  double uy=0.;
498  double tanphi=0.;
499  if((fabs(phi) < 7.85398163397448279e-01) || (fabs(phi - 3.1415926535897931) < 7.85398163397448279e-01))
500  {
501  tanphi = tan(phi);
502  ux = sqrt(1./(1. + tanphi*tanphi));
503  if(phi > 1.57079632679489656 && phi < 4.71238898038468967){ux=-ux;}
504  uy = ux*tanphi;
505  }
506  else
507  {
508  tanphi = tan(1.57079632679489656 - phi);//cotangent
509  uy = sqrt(1./(tanphi*tanphi + 1.));
510  if(phi > 3.1415926535897931){uy=-uy;}
511  ux = uy*tanphi;
512  }
513  double dx = ux*d;
514  double dy = uy*d;
515  double dispx = ux/kappa;
516  double dispy = uy/kappa;
517  double cx = dx + dispx;
518  double cy = dy + dispy;
519 
520  double x=0.;
521  double y=0.;
522  double z=0.;
523 
524  double startx=dx;
525  double starty=dy;
526  double startz=z0;
527  nhits=0;
528  for(unsigned int l=0;l<radii.size();l++)
529  {
530  if(intersect_circles(hel, startx, starty, radii[l], 1./kappa, cx, cy, x, y) == true)
531  {
532  layer[nhits]=l;
533  x_hits[nhits]=x;
534  y_hits[nhits]=y;
535 
536  double k = kappa;
537  double D = sqrt((startx-x)*(startx-x) + (starty-y)*(starty-y));
538 
539  double s=0.;
540  if(0.5*k*D > 0.1)
541  {
542  s = 2.*asin(0.5*k*D)/k;
543  }
544  else
545  {
546  double temp1 = k*D*0.5;temp1*=temp1;
547  double temp2 = D*0.5;
548  s += 2.*temp2;
549  temp2*=temp1;
550  s += temp2/3.;
551  temp2*=temp1;
552  s += (3./20.)*temp2;
553  temp2*=temp1;
554  s += (5./56.)*temp2;
555  }
556  double dz = sqrt(s*s*dzdl*dzdl/(1. - dzdl*dzdl));
557  if(dzdl>0.){z = startz + dz;}
558  else{z = startz - dz;}
559 
560  z_hits[nhits]=z;
561 
562  nhits++;
563 
564  startx = x;
565  starty = y;
566  startz = z;
567  }
568  }
569  if(nhits>0)
570  {
571  trk_kappa=kappa;
572  trk_d=d;
573  trk_phi=phi;
574  trk_dzdl=dzdl;
575  trk_z0=z0;
576  charge=0;
577  if(hel==true){charge=1;}
578  track=i;
579  otree->Fill();
580  }
581  }
582  */
583 
584 
585  etree->Fill();
586  }
587 
588 
589  etree->Write();
590  ofile->Close();
591  ofile->Delete();
592 
593 
594 
595  return 0;
596 }
597 
598