Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fasthadv4.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fasthadv4.C
1 #include "fasthadv4.h"
2 
3 /*
4  First shot at a fast hadron simulation tool.
5  The "detector" is rectangular, nbin*nbin granularity, ddepth deep.
6  Dimensions are arbitrary (make sense only in context with the
7  electromagnetic and hadronic shower width parameters dwidthem, dwidthhad.
8  Impact point and theta, phi angles are arbitrary.
9  Each hadron enters the detector as minimum ionizing; afetr each (fixed)
10  MIP deposit throws a random number to decide whether to start a shower or not.
11  Once it decided to form a shower, the remaining length in the detector is
12  calculated, and serves as a weight for the energy it will deposit (the deposit
13  in any single step is constant now, so one plays with the number of steps to change
14  the total deposit).
15  Also, the total energy it will actually deposited is weighted (inversely) with
16  a formula, essentially with the log of the original hadron energy.
17  The total energy it will deposit is "dshe".
18  Next, it will be decided how much of "dshe" will be electromagnetic and
19  how much is hadronic ("reh"). Since the electromagnetic part usually starts
20  later, a spatial delay is assigned to the electromagnetic shower.
21  Finally the two showers (em, had) are generated in three dimensions,
22  gradually widening with depth, but randomized with a Gaussian.
23  */
24 
25 // Impact point and angle
26 double dimpx,dimpy,dimpxtheta,dimpyphi;
27 double dpath,dfullpath,ddepth,dtote,dtotemip,dtoteem,dtotehad;
28 double x,y,z,xhad,yhad,zhad,reh;
29 double dshe,dhade,deme;
30 double ddzhad,ddzem; // stepsize for had, em part (mip is fixed)
31 double dfudge;
32 double dcostheta,dcosphi,dsintheta,dsinphi;
33 double demdelay;
34 
35 
36 int nmipstep,nshstep;
37 int i,j;
38 
39 
40 int ntrysh,ntryem,ntryhad;
41 
42 
43  TH2D *h2dhits = new TH2D("h2dhits","Hit pattern",nbin,xlo,xhi,nbin,ylo,yhi);
44  TH2D *h2dlive = new TH2D("h2dlive","Live display",nbin,xlo,xhi,nbin,ylo,yhi);
45  TH1D *htote = new TH1D("htote","Total energy",100,0.0,15.0);
46  TH1D *htotemip = new TH1D("htotemip","Total energy",25,0.0,1.0);
47  TH1D *htoteem = new TH1D("htoteem","Total energy",160,0.0,8.0);
48  TH1D *htotehad = new TH1D("htotehad","Total energy",160,0.0,8.0);
49  TH2D *h2dmipnonmip = new TH2D("h2dmipnonmip","MIP vs non MIP energy",
50  80,0.0,8.0,50,0.0,0.5);
51  TH3D *h3dlivemip = new TH3D("h3dlivemip","3D deposit",
52  nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
53  TH3D *h3dliveem = new TH3D("h3dliveem","3D em deposit",
54  nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
55  TH3D *h3dlivehad = new TH3D("h3dlivehad","3D had deposit",
56  nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
57  TH3D *h3dliveall = new TH3D("h3dliveall","3D all deposit",
58  nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
59 
60 
61 
62 void fasthadv4(double dimpx=0.5, double dimpy=0.5,
63  double dimpxtheta=0.0, double dimpyphi=0.0)
64 {
65 
66 
67  nmipstep = ndefmipstep;
68  nshstep = ndefshstep;
69 
70  dcostheta = TMath::Cos(dimpxtheta);
71  dcosphi = TMath::Cos(dimpyphi);
72  dsintheta = TMath::Sin(dimpxtheta);
73  dsinphi = TMath::Sin(dimpyphi);
74  if(dcostheta<=0.0 || dcosphi<=0.0)
75  {
76  std::cout << "Something dead wrong with impact angle " << std::endl;
77  return;
78  }
79 
80  nmipstep = int(float(ndefmipstep)/(dcostheta*dcosphi));
81  nshstep = int(float(ndefshstep)/(dcostheta*dcosphi));
82  dfullpath = dlength/(dcostheta*dcosphi);
83  std::cout << nmipstep << " " << nshstep << " " << dfullpath << std::endl;
84 
85 
86 
87  dhade = dtothade / float(nshstep);
88  deme = dtothade / float(nshstep);
89  ddzem = dlength / float(ndefshstep);
90  ddzhad = dlength / float(ndefshstep);
91 
92 TCanvas * c10 = new TCanvas("c10","Live display",1200,600);
93  c10->Divide(2,1);
94 
95 
96  for(i=1; i<=ntry;i++)
97  {
98 
99  dpath = ddepth = dtote = dtotemip = dtoteem = dtotehad = dshe = 0.0;
100  h3dlivemip->Reset();
101  h3dliveem->Reset();
102  h3dlivehad->Reset();
103  h2dlive->Reset();
104 
105  for(j=1; j<=nmipstep; j++)
106  {
107  dpath = j*ddzmip;
108  x = dimpx + dpath*dsintheta;
109  y = dimpy + dpath*dsintheta;
110  z = dpath * dcostheta*dcosphi;
111  dtotemip +=ddedx;
112  h2dhits->Fill(x,y,ddedx);
113  // i%nlive == 0 ? h2dlive->Fill(x,y,ddedx) : ;
114  if(i%nlive == 0)
115  {
116  h2dlive->Fill(x,y,ddedx);
117  h3dlivemip->Fill(x,y,z);
118  }
119 
120  if(gRandom->Uniform()>dmipstop) break;
121  }
122  dtotemip *= 1.0 + dresol * gRandom->Gaus();
123  htotemip->Fill(dtotemip);
124  xhad = x; yhad = y;
125  zhad = dpath * dcostheta*dcosphi;
126  ddepth = dpath * dcostheta*dcosphi;
127  if(ddepth>=dlength)
128  {
129  htote->Fill(dtotemip);
130  continue; // done with this, generate the next particle
131  }
132 
133  // Get the remaining energy
134  dshe = dtothade - dtotemip;
135  // Estimate how much of it will be deposited (function of dtothade)
136  dshe *= TMath::Sqrt(1.0/(1.0+0.5*TMath::Log(dshe+1.0)));
137 
138  // Right now for simplicity the deposits per step are identical,
139  // so you manipulate the number of steps to get the proper deposit
140  ntrysh = int(nshstep * dshe/dtothade);
141 
142  // Ratio of electromagnetic and hadronic part of the shower
143  reh = 0.1 + 0.15*gRandom->Uniform();
144  reh /=(1.0+0.6*TMath::Log(dshe+1.0));
145 
146  ntryem = int(ntrysh * reh);
147  ntryhad = int(ntrysh * (1.0 - reh));
148  ntryem = int(ntryem * ((dlength-ddepth) / dlength));
149  ntryhad = int(ntryhad * ((dlength-ddepth) / dlength));
150 
151  // The em part of the hadron shower doesn't start immediately
152  // it is
153  demdelay = 2.0 + gRandom->Gaus();
154  demdelay > 0.0 ? : demdelay = -demdelay;
155  if(ntryem>0) ddzem = (dlength-ddepth)/float(ntryem);
156 
157 
158  dpath = 0.0;
159  for(j=1; j<ntryem; j++)
160  {
161  dpath = j*ddzem + demdelay;
162  if((ddepth + dpath * dcostheta*dcosphi)>dlength)
163  {
164  htote->Fill(dtoteem);
165  break;
166  }
167  x = xhad + dpath*(dsintheta+dwidthem*gRandom->Gaus());
168  y = yhad + dpath*(dsinphi+dwidthem*gRandom->Gaus());
169  z = zhad+dpath * dcostheta*dcosphi;
170  dtoteem += deme;
171  h2dhits->Fill(x,y,deme);
172  if(i%nlive == 0)
173  {
174  h2dlive->Fill(x,y,deme);
175  h3dliveem->Fill(x,y,z);
176  }
177  }
178  dtoteem *= 1.0 + dresol * gRandom->Gaus();
179  htoteem->Fill(dtoteem);
180 
181  if(ntryhad>0) ddzhad = (dlength-ddepth)/float(ntryhad);
182 
183 
184  dpath = 0.0;
185  for(j=1; j<ntryhad; j++)
186  {
187  dpath = j*ddzhad;
188  if((ddepth + dpath * dcostheta*dcosphi)>dlength)
189  {
190  htote->Fill(dtotehad);
191  break;
192  }
193  x = xhad + dpath*(dsintheta+dwidthhad*gRandom->Gaus());
194  y = yhad + dpath*(dsinphi+dwidthhad*gRandom->Gaus());
195  z = zhad+dpath * dcostheta*dcosphi;
196  dtotehad += dhade;
197  h2dhits->Fill(x,y,dhade);
198  if(i%nlive == 0)
199  {
200  h2dlive->Fill(x,y,dhade);
201  h3dlivehad->Fill(x,y,z);
202  }
203  }
204 
205  dtotehad *= 1.0 + dresol * gRandom->Gaus();
206  htotehad->Fill(dtotehad);
207 
208  htote->Fill(dtotemip+dtoteem+dtotehad);
209  h2dmipnonmip->Fill(dtoteem+dtotehad,dtotemip);
210 
211  if(i%nlive==0)
212  {
213  c10->cd(1);
214  h2dlive->Draw("colz");
215  c10->cd(2);
216  h3dlivemip->SetMarkerStyle(20);
217  h3dlivemip->SetMarkerColor(1);
218  h3dlivemip->SetMarkerSize(0.2);
219  h3dlivemip->Draw();
220  h3dliveem->SetMarkerStyle(20);
221  h3dliveem->SetMarkerColor(2);
222  h3dliveem->SetMarkerSize(0.4);
223  h3dliveem->Draw("SAME");
224  h3dlivehad->SetMarkerStyle(20);
225  h3dlivehad->SetMarkerColor(3);
226  h3dlivehad->SetMarkerSize(0.5);
227  h3dlivehad->Draw("SAME");
228  c10->Update();
229  std::cout << i << std::endl;
230  // if(i==300) c10->Print("evt300.jpg","jpg");
231 
232  h3dlivemip->Reset();
233  h3dliveem->Reset();
234  h3dlivehad->Reset();
235  h2dlive->Reset();
236 
237 
238 
239  }
240 
241  }
242 
243 
244 TFile *fout = new TFile("fasthadv4_out.root","RECREATE");
245 
246  h2dhits->Write();
247  htote->Write();
248  htotemip->Write();
249  htoteem->Write();
250  htotehad->Write();
251  h2dmipnonmip->Write();
252  h2dlive->Write();
253  h3dlivemip->Write();
254  h3dliveem->Write();
255  h3dlivehad->Write();
256  h3dliveall->Write();
257 
258 
259 
260  fout->Close();
261 
262 }
263 
264 
265 
266 
267 
268 
269 
270