Pluto v6.3.0
Loading...
Searching...
No Matches
PBremsstrahlung.cc
Go to the documentation of this file.
1
2//
3// Parameterization of the virtual pn/pp Bremsstrahlung
4// Ref. 19 and Ref. 20
5//
6//
7// Author: Wuestenfeld/Dohrmann
9
10
11#include "PBremsstrahlung.h"
12#include <cmath>
13
14Double_t PBremsstrahlung::Pol1(Double_t z, Double_t par1, Double_t par2) const {
15 Double_t x = z;
16 return par1+par2*x;
17}
18
19Double_t PBremsstrahlung::Pol2(Double_t z, Double_t par1, Double_t par2, Double_t par3) const {
20 Double_t x = z;
21 return par1+par2*x+par3*x*x;
22}
23
24#if 0
25Double_t PBremsstrahlung::Pol3(Double_t z, Double_t par1, Double_t par2, Double_t par3, Double_t par4)
26{
27 Double_t x=z;
28 return par1+par2*x+par3*x*x+par4*x*x*x;
29}
30#endif
31
32Double_t PBremsstrahlung::DLines(Double_t x, Double_t par0, Double_t par1, Double_t par2, Double_t par3) const {
33 Double_t m,n;
34
35 if(x < 1.9) {
36 m = (par1-par0)/(0.65);
37 n = par0-m*1.25;
38 return m*x+n;
39 } else {
40 if(x < 2.2) {
41 m = (par2-par1)/(0.3);
42 n = par1-m*1.9;
43 return m*x+n;
44 } else {
45 m = (par3-par2)/(1.3);
46 n = par2-m*2.2;
47 return m*x+n;
48 }
49 }
50}
51
52Double_t PBremsstrahlung::AGauss(Double_t x, Double_t *par) const
53{
54 Double_t result;
55
56 if(par[0] == 0) return 0;
57
58 if(x-par[2] > 0) {
59 result = par[0]*exp(-0.5*pow((x-par[2])/(TMath::Max(1.e-10,((1.0 + par[3])*par[1]))),2));
60 } else {
61 result = par[0]*exp(-0.5*pow((x-par[2])/(TMath::Max(1.e-10,((1.0 - par[3])*par[1]))),2));
62 }
63 //cout << x << "," << par[0] << "," << par[1] << "," << par[2] << "," << par[3] << "," <<result << endl;
64 //cout << (x-par[2]) << "," << ((1.0 + par[3])*par[1]) << "," << exp(-0.5*pow((x-par[2])/(TMath::Max(1.e-10,((1.0 - par[3])*par[1]))),2)) << endl;
65 return result;
66}
67
68
69PDistribution *PBremsstrahlung::Clone(const char*) const {
70 //clone the object
71 return new PBremsstrahlung((const PBremsstrahlung &)* this);
72};
73
74PBremsstrahlung::PBremsstrahlung(const Char_t *id, const Char_t *de, Int_t key) :
75 PChannelModel(id, de, key) {
76 //Constructor
77
78 if (is_channel < 0)
79 Fatal("PBremsstrahlung", "The model (%s) should be bound to CHANNELS only", de);
80
81 //try to figure out if the model is pp or pn
82 //in addition make some consistency checks
83 Int_t tid[11];
84 tid[0] = 10;
85
86 makeStaticData()->GetDecayModeByKey(primary_key, tid); // retrieve current mode info
87
88 if (tid[0] != 3) {
89 Fatal("PBremsstrahlung", "Only 3 body decay");
90 }
91 dilepton_position = -1;
92 neutron_position = -1;
93 for (int i=1; i<=3; i++) {
94 if (makeStaticData()->IsParticle(tid[i],"dilepton")) dilepton_position = i;
95 if (makeStaticData()->IsParticle(tid[i],"n")) neutron_position = i;
96 }
97
98 if (dilepton_position == -1) {
99 Fatal("PBremsstrahlung", "No dilepton found");
100 }
101 threshold = 2*makeStaticData()->GetParticleMass("p");
102 if (neutron_position != -1) {
103 threshold = makeStaticData()->GetParticleMass("p") + makeStaticData()->GetParticleMass("n");
104 }
105
106 p2_energy = 2.2;
107 bin = NULL;
108
109 //...just in the case that they are needed
110 mp = makeStaticData()->GetParticleMass("p");
111 mn = makeStaticData()->GetParticleMass("n");
112 model = 0;
113 author = 0;
114 graph = NULL;
115 graph2d = NULL;
116}
117
118void PBremsstrahlung::SetMode(Char_t mode) {
119 // model [0..4] where 0 => bare gNN Bremsstrahlung [default]
120 // 1 => gDN Delta contribution included
121 // 2 => cSum coherent sum of bare an gDN Bremsstrahlung
122 // 3 => N1520 contribution (Shyam/Mosel only)
123 // The following modes are untested (IF, 26.9.2008)
124 // 4 => FSI final state interactions included
125 // 5 => VMD Vector Meson Dominance contribution
126 // 6 => VMD and FSI contribution included
127 model = mode;
128}
129
130
131Bool_t PBremsstrahlung::Init(void) {
132 //looking for the dilepton. This is mandatory
133 dilepton = GetParticle("dilepton");
134 if (!dilepton) {
135 Warning("Init", "Dilepton not found");
136 return kFALSE;
137 }
138
139 parent = GetParticle("parent");
140 if (!parent) {
141 Warning("Init", "Parent not found");
142 return kFALSE;
143 }
144
145 //get in/outgoing nucleons...
146 p1 = GetParticle("p1");
147 p2 = GetParticle("p2");
148 p3 = GetParticle("p3");
149 p4 = GetParticle("p4");
150
151 return kTRUE;
152}
153
154Double_t PBremsstrahlung::EvalPar(const double *x, const double *) {
155 return Eval(x[0]);
156}
157
158
159Double_t PBremsstrahlung::EvalSM(Double_t x, Double_t, Double_t, Double_t) const {
160 Double_t kinLim;
161 Double_t M = x;
162 Double_t b0=0, b1=0, b2=0, b3=0;
163
164 if(neutron_position > -1) {
165 // p + n case
166
167 kinLim = TMath::Sqrt(pow((mn + mp),2)+2*mp* p2_energy)-(mn+mp);
168 if(M>kinLim) { // FD this is somewhat clumsy, but ...
169 return 0;
170 }
171 if(model == gNN) {
172 b0 = Pol2(p2_energy, 24.4393, -7.58498, 0.819536);
173 b1 = Pol2(p2_energy, -33.4414, 14.8302, -1.75497);
174 b2 = Pol2(p2_energy, 4.13174, -5.99696, 1.95849);
175 b3 = Pol2(p2_energy, 0.942944, 1.5077, -0.207335);
176 } else {
177 if(model == gDN) {
178 b0 = Pol2(p2_energy, 5.40183, 4.9095, -1.06483);
179 b1 = Pol2(p2_energy, 27.9443, -25.8368, 4.04411);
180 b2 = Pol2(p2_energy, 0.859852, -1.2709, 0.437889);
181 b3 = Pol2(p2_energy, 1.3053, 0.788249, -0.109513);
182 } else {
183 if(model == gN1520) {
184 b0 = Pol2(p2_energy, 41.7807, -22.0668, 2.91033);
185 b1 = Pol2(p2_energy, -41.2204, 22.8689, -3.01134);
186 b2 = Pol2(p2_energy, 6.69874, -9.94155, 3.39074);
187 b3 = Pol2(p2_energy, 0.795684, 1.37388, -0.172581);
188 //cout << b0 << " " << b1 << " " << b2 << " " << b3 << endl;
189 } else {
190 if(model == cSum) {
191 b0 = Pol2(p2_energy,7.65485, 3.43015, -0.811125);
192 b1 = Pol2(p2_energy, 14.7932, -15.5442, 2.47197);
193 b2 = Pol2(p2_energy, 0.379128, -0.559714, 0.193903);
194 b3 = Pol2(p2_energy, 1.60404, 0.522244, -0.0678947);
195 } else {
196 Error("EvalSM","Invalid Bremsstrahlung model!");
197 return 0;
198 }
199 }
200 }
201 }
202 //cout << x << endl;
203 return (1./1000000.)*pow(fabs(M - kinLim),b3)/fabs((b2*exp(b0*M+b1*M*M)));
204 } else {
205 // p + p case
206 kinLim = TMath::Sqrt(pow((mp + mp),2)+2*mp* p2_energy)-(mp+mp);
207 if(M>kinLim) { // FD this is somewhat clumsy again, but ...
208 Error("EvalSM", "M>kinLim");
209 return 0;
210 }
211 if(model == gNN) {
212 b0 = Pol2(p2_energy, 27.4811, -10.6609, 1.25601);
213 b1 = Pol2(p2_energy, -38.7574, 19.8781, -2.52863);
214 b2 = Pol2(p2_energy, 9.29954, -13.6681, 4.60585);
215 b3 = Pol2(p2_energy, 1.85104, 0.518799, -0.0646308);
216 } else {
217 if(model == gDN) {
218 b0 = Pol2(p2_energy, 2.13375, 7.31236, -1.44497);
219 b1 = Pol2(p2_energy, 49.7263, -40.9453, 6.24139);
220 b2 = Pol2(p2_energy, 1.93576, -2.96657, 1.02646);
221 b3 = Pol2(p2_energy, 0.621623, 1.22859, -0.160213);
222 } else {
223 if(model == gN1520) {
224 b0 = Pol2(p2_energy, 43.3482, -24.036, 3.24966);
225 b1 = Pol2(p2_energy, -43.0502, 26.0555, -3.58266);
226 b2 = Pol2(p2_energy, 8.37618, -12.8982, 4.7091);
227 b3 = Pol2(p2_energy, 1.47027, 0.557303, -0.0442349);
228 } else {
229 if(model == cSum) {
230 b0 = Pol2(p2_energy, 6.2699, 3.50518, -0.808666);
231 b1 = Pol2(p2_energy, 30.9533, -25.3514, 3.81197);
232 b2 = Pol2(p2_energy, 2.07252, -3.04563, 1.02282);
233 b3 = Pol2(p2_energy, 2.69517, -0.512592, 0.119903);
234 } else {
235 Error("EvalSM","Invalid Bremsstrahlung model!");
236 return 0;
237 }
238 }
239 }
240 }
241 return (1./1000000.)*pow(fabs(M - kinLim),b3)/fabs((b2*exp(b0*M+b1*M*M)));
242 }
243}
244
245
246Double_t PBremsstrahlung::Eval(Double_t x, Double_t y , Double_t z , Double_t t) const {
247
248 if (graph) return graph->Eval(x);
249 if (graph2d) return graph2d->Interpolate(x,p2_energy);
250
251 if (author == BREMS_SHYAM_MOSEL) return EvalSM(x, y, z, t);
252
253 Double_t kinLim;
254 Double_t M = x;
255 Double_t b0=0, b1=0, b2=0, b3=0;
256 Double_t parA[4];
257 Double_t parB[4];
258
259 if(neutron_position > -1) {
260 // p + n case
261 //cout << p2_energy << endl;
262 kinLim = TMath::Sqrt(pow((mn + mp),2)+2*mp* p2_energy)-(mn+mp);
263 if(M>kinLim) { // FD this is somewhat clumsy, but ...
264 return 0;
265 }
266 if(model == gNN) {
267 b0 = Pol2(p2_energy, 16.9411, -5.93276, 0.7807);
268 b1 = Pol2(p2_energy, -33.19, 18.8141, -2.91396);
269 b2 = Pol2(p2_energy, 0.0654566, -0.102335, 0.0592914);
270 b3 = Pol2(p2_energy, 1.99472, 0.447733, -0.0727639);
271 } else {
272 if(model == gDN) {
273 b0 = Pol2(p2_energy, 20.240726, -7.9911127, 0.94649249);
274 b1 = Pol2(p2_energy, -27.103907, 14.986975, -2.3138011);
275 b2 = Pol2(p2_energy, 0.057434153, -0.12646526, 0.07317695);
276 b3 = Pol2( p2_energy, 4.2753172, -1.5160941, 0.33231404);
277 } else {
278 if(model == cSum) {
279 b0 = Pol2( p2_energy, 12.370295, -2.379364, 0.045619611);
280 b1 = Pol2( p2_energy, -10.650187, 2.7560029, -0.16295315);
281 b2 = Pol2( p2_energy, 0.041419424, -0.069788419, 0.038073331);
282 b3 = Pol2( p2_energy, 1.9655238, 0.43345186, -0.04576372);
283 } else {
284 if(model == FSI) {
285 //???
286 } else {
287 if(model == VMD) {
288 } else {
289 if(model == cFsiVmd) {
290 b0 = DLines(p2_energy, 11.73, 9.873, 6.253, 7.787);
291 b1 = DLines(p2_energy, -7.589, -8.468, -30.1, -11.75);
292 b2 = DLines(p2_energy, 0.01877, 0.03573, 0.007521, 0.2954);
293 b3 = DLines(p2_energy, 1.595, 2.317, 11.43, 7.124);
294
295 if(p2_energy > 1.9) {
296 parA[0] = Pol2(p2_energy, -0.65, 0.09, 0.116);
297 parA[1] = Pol1(p2_energy, -0.008866, 0.00908);
298 parA[2] = Pol1(p2_energy, 0.7693, 0.004265);
299 parA[3] = Pol1(p2_energy, 0.5262, -0.2392);
300
301 parB[0] = Pol2(p2_energy, -0.7279, 0.4763, -0.04907);
302 parB[1] = Pol1(p2_energy, -0.1483, 0.09447);
303 parB[2] = Pol1(p2_energy, 0.7329, 0.009148);
304 parB[3] = Pol1(p2_energy, 0.08683, -0.2189);
305 } else {
306 for(Int_t i=0; i<4; i++) {
307 parA[i] = 0;
308 parB[i] = 0;
309 }
310 }
311 return (1./1000000.)*pow(fabs(M - kinLim),b3)/(b2*exp(b0*M+b1*M*M)) + AGauss(p2_energy,parA) + AGauss(p2_energy,parB);
312 } else {
313 Error("Eval", "Invalid Bremsstrahlung model!");
314 return 0;
315 }
316 }
317 }
318 //return pow(fabs(M - kinLim),b3)/(b2*exp(b0*M+b1*M*M)) + AGauss(p2_energy,parA) + AGauss(p2_energy,parB);
319 }
320 }
321 }
322 return (1./1000000.)*pow(fabs(M - kinLim),b3)/(b2*exp(b0*M+b1*M*M));
323 //FD probably could omit 'fabs' function here
324 } else {
325 // p + p case
326
327 kinLim = TMath::Sqrt(pow((mp + mp),2)+2*mp* p2_energy)-(mp+mp);
328 if(M > kinLim) { // FD this is somewhat clumsy again, but ...
329 return 0;
330 }
331
332 //FD 24/01/08 Take either this block (note return function)
333 //FD do NOT delete this set of parameters
334 /* b0 = Pol2(Teff,33.6382,-15.3112,2.10402);
335 b1 = Pol2(Teff,-48.9918,28.9945,-4.43253);
336 b2 = Pol3(Teff,-0.700971,1.71958,-1.39966,0.389064);
337 b3 = Pol2(Teff,2.16083,0.835047,-0.216302);
338 return pow(fabs(M *M - kinLim * kinLim),b3)/(b2*exp(b0*M+b1*M*M));
339 //FD cannot omit 'fabs' function there
340 */
341
342 //FD 24/01/08 Or that block (note different return function)
343 if(model == gNN) {
344 b0 = Pol2(p2_energy, 28.572, -15.518, 2.41939);
345 b1 = Pol2(p2_energy, -61.4211, 43.4459, -7.57627);
346 b2 = Pol2(p2_energy, 0.226395, -0.510502, 0.35198);
347 b3 = Pol2(p2_energy, 4.2694, -1.5697, 0.338622);
348 } else {
349 if(model == gDN) {
350 b0 = Pol2(p2_energy, 17.5256, -5.47771, 0.389949);
351 b1 = Pol2(p2_energy, -31.254, 18.2722, -2.80298);
352 b2 = Pol2(p2_energy, 0.0928247, -0.188658, 0.102528);
353 b3 = Pol2(p2_energy, 4.63276, -1.19151, 0.210233);
354 } else {
355 if(model == cSum) {
356 b0 = Pol2(p2_energy, 16.1204, -6.72895, 0.927332);
357 b1 = Pol2(p2_energy, -13.0257, 4.86514, -0.614537);
358 b2 = Pol2(p2_energy, 0.0421482, -0.0873232, 0.0551647);
359 b3 = Pol2(p2_energy, 2.48228, 0.668239, -0.0493283);
360 } else {
361 if(model == FSI) {
362
363 } else {
364 if(model == VMD) {
365 } else {
366 if(model == cFsiVmd) {
367 } else {
368 Error("Eval","Invalid Bremsstrahlung model!");
369 return 0;
370 }
371 }
372 }
373 }
374 //return pow(fabs(M - kinLim),b3)/(b2*exp(b0*M+b1*M*M)) + AGauss(p2_energy,parA) + AGauss(p2_energy,parB);
375 }
376 }
377 return (1./1000000.)*pow(fabs(M - kinLim),b3)/(b2*exp(b0*M+b1*M*M));
378 //FD probably could omit 'fabs' function here
379 }
380}
381
382Double_t PBremsstrahlung::GetWeight(void) {
383 //boost particle 2 such to be in particle 1 rest frame
384
385 Double_t m1 = p1->M();
386 Double_t m2 = p2->M();
387
388 Double_t mq = parent->M();
389
390 //m2 should be the target:
391 p2_energy = (mq*mq-m1*m1-m2*m2-2*m1*m2)/(2*m2);
392
393 if (p1->Is("n")) //neutron at rest
394 p2_energy = (mq*mq-m1*m1-m2*m2-2*m1*m2)/(2*m1);
395
396 Double_t val = Eval(dilepton->M());
397 //cout << p2_energy << ":" << dilepton->M() << ":" << val << endl;
398 return val;
399}
400
401Bool_t PBremsstrahlung::SampleMass(void) {
402
403 //boost particle 2 such to be in particle 1 rest frame
404
405 PParticle tmp(p2); //particle under investigation. Make better a copy
406 tmp.Boost(-p1->BoostVector());
407 p2_energy = tmp.KE();
408
409 if (!bin) {
410 bin = new PAdaptiveMesh(0, 3, 0.5);
411
412 bin->SetTF1((TF1*)this);
413//bin->Divide(5,3);
414 bin->Divide(3, 3);
415 bin->Print();
416 }
417
418 //dilepton->SetM(this->GetRandom());
419
420 dilepton->SetM(bin->GetRandom());
421
422 return kTRUE;
423}
424
425
426ClassImp(PBremsstrahlung)
ClassImp
ClassImp(PBeamLineSimulation) PBeamLineSimulation
Definition PBeamLineSimulation.cc:22
BREMS_SHYAM_MOSEL
#define BREMS_SHYAM_MOSEL
Definition PBremsstrahlung.h:16
makeStaticData
PStaticData * makeStaticData()
par
PParticle * par[11]
Definition analyze_g12C.C:15
par2
PParticle * par2[3]
Definition analyze_pp_delta_dalitz.C:11
PBremsstrahlung::Clone
PDistribution * Clone(const char *delme=NULL) const
Definition PBremsstrahlung.cc:69
PBremsstrahlung::EvalPar
Double_t EvalPar(const Double_t *x, const Double_t *params)
Definition PBremsstrahlung.cc:154
PBremsstrahlung::GetWeight
Double_t GetWeight(void)
Definition PBremsstrahlung.cc:382
PBremsstrahlung::SampleMass
Bool_t SampleMass(void)
Definition PBremsstrahlung.cc:401
PBremsstrahlung::SetMode
void SetMode(Char_t mode)
Definition PBremsstrahlung.cc:118
PBremsstrahlung::PBremsstrahlung
PBremsstrahlung(const Char_t *id, const Char_t *de, Int_t key)
Definition PBremsstrahlung.cc:74
PBremsstrahlung::Eval
Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Definition PBremsstrahlung.cc:246
PDistribution::GetParticle
PParticle * GetParticle(const Char_t *pflag=NULL)
PParticle::KE
virtual Double_t KE() const
Definition PParticle.h:102
PStaticData::GetParticleMass
double GetParticleMass(const int &id)
PStaticData::GetDecayModeByKey
void GetDecayModeByKey(const int, int *n)
m1
PProjector * m1
Definition dp_eta_subthreshold.C:18
n
Int_t n
Definition etap+omeg.C:21