Pluto v6.3.0
Loading...
Searching...
No Matches
PBremsstrahlungPlugin.cc
Go to the documentation of this file.
1
2//Plugin for the Bremsstrahlung from Kaptari/Kaempfer [L1] and
3//Shyam/Mosel [L2]
4//
5//It can be activated with
6//makeDistributionManager()->Exec("brems: kaptari");
7//
8//Commands can be formwarded to the plugin via:
9//makeDistributionManager()->Exec("brems: command");
10//
11//The plugin supports the following commands:
12//elastic Elastic term only
13//delta Delta term only
14//sum Coherent sum
15//fsi Add FSI a la Jost function
16//
17//For the absolute normalization it is suggested to use the weighting option:
18//
19//kin_max = arg arg=0...1 maximal kinematic range for di-lepton sampling
20//weighting activate weighting
21//
22//Example:
23//makeDistributionManager()->Exec("brems: kaptari ; kin_max=0.9 ; weighting");
24//
25//Author: I. Froehlich
26//Written: 27.7.2008
27//
28//References:
29//[L1] L.P. Kaptari and B. Kaempfer, Nucl. Phys. A 764 (2006) 338, arXiv:nucl-th/0504072
30//[L2] R. Shyam, U. Mosel (Giessen U.) , Phys.Rev.C67:065202,2003, hep-ph/0303035, and priv. comm
32
33
34#include "PBremsstrahlungPlugin.h"
35
36
37PBremsstrahlungPlugin::PBremsstrahlungPlugin(const Char_t *id, const Char_t *de):
38 PDistributionCollection(id, de) {
39
40 pn_gen = NULL;
41 pp_gen = NULL;
42 kin_max = 1.;
43 pn_fsi = pp_fsi = np_fsi = NULL;
44 pn_kk = pp_kk = pn_sm_total = pp_sm_total = NULL;
45 pn_sm_elastic = pp_sm_elastic = pn_sm_delta = pp_sm_delta
46 = pn_sm_n1520 = pp_sm_n1520 = NULL;
47 pn = pn_kk;
48 pp = pp_kk;
49 array_pp_sm_125 = array_pn_sm_125 = NULL;
50
51 RequiresPlugin("elementary");
52}
53
54Bool_t PBremsstrahlungPlugin::Activate(void) {
55
56 PDistributionManagerUtil * pdmutil = makeDistributionManagerUtil();
57
58 //First, we check if all composite particles are available
59
60// if (!makeStaticData()->IsParticleValid("p + p")) {
61// makeStaticData()->AddParticle(14014, "p + p",2.);
62// }
63// if (!makeStaticData()->IsParticleValid("p + n")) {
64// makeStaticData()->AddParticle(13014, "p + n",2.);
65// }
66// if (!makeStaticData()->IsParticleValid("n + p")) {
67// makeStaticData()->AddParticle(14013, "n + p",2.);
68// }
69// if (!makeStaticData()->IsParticleValid("d + p")) {
70// makeStaticData()->AddParticle(14045, "d + p",2.);
71// }
72// if (!makeStaticData()->IsParticleValid("p + d")) {
73// makeStaticData()->AddParticle(45014, "p + d",2.);
74// }
75
76 Int_t ipid[4];
77 ipid[0] = makeStaticData()->GetParticleID("p + n");
78 ipid[1] = makeStaticData()->GetParticleID("p");
79 ipid[2] = makeStaticData()->GetParticleID("n");
80 ipid[3] = makeStaticData()->GetParticleID("dilepton");
81
82 //Add decay
83 if (makeStaticData()->GetDecayKey(ipid, 3) < 0)
84 makeStaticData()->AddDecay(-1, "p + n -> p + n + dilepton (Bremsstrahlung)",
85 "p + n", "p,n,dilepton", 1.0 );
86
87 ipid[0] = makeStaticData()->GetParticleID("n + p");
88 if (makeStaticData()->GetDecayKey(ipid, 3) < 0)
89 makeStaticData()->AddDecay(-1, "n + p -> p + n + dilepton (Bremsstrahlung)",
90 "n + p", "p,n,dilepton", 1.0 );
91
92 ipid[0] = makeStaticData()->GetParticleID("p + p");
93 ipid[2] = makeStaticData()->GetParticleID("p");
94 if (makeStaticData()->GetDecayKey(ipid, 3) < 0)
95 makeStaticData()->AddDecay(-1, "p + p -> p + p + dilepton (Bremsstrahlung)",
96 "p + p", "p,p,dilepton", 1.0 );
97 //deuteron beam
98 ipid[0] = makeStaticData()->GetParticleID("d + p");
99 ipid[2] = makeStaticData()->GetParticleID("n + p");
100 if (makeStaticData()->GetDecayKey(ipid, 2) < 0)
101 makeStaticData()->AddDecay(-1, "d + p -> p + (n + p) (Quasi-Free)",
102 "d + p", "p,n + p", 1.0 );
103
104 pdmutil->AddSubGroup("bremsstrahlung", "Bremsstrahlung", "root");
105 pdmutil->SetGroup("bremsstrahlung");
106
107 return kTRUE;
108}
109
112
113
114Bool_t PBremsstrahlungPlugin::ExecCommand(const char *command, Double_t value) {
115
116 if (strcmp (command,"weighting") == 0) {
117
118 //Monto-Carlo integration
119 pn->EnableWeighting();
120 pn->SetExpectedWeightMean(-1);
121 pn->SetVersionFlag(VERSION_GENERATOR_MC);
122
123 pp->EnableWeighting();
124 pp->SetExpectedWeightMean(-1);
125 pp->SetVersionFlag(VERSION_GENERATOR_MC);
126
127
128 if (pn_gen) return kTRUE;
129
130 //TF1 object representing the di-lepton statistics:
131 TF1 *flat = new TF1("flat", "1", 0, kin_max);
132
133 //The "PInclusiveModel" can be used as a generator:
134 pn_gen =
135 new PInclusiveModel("flat@p + n_brems_p_n_dilepton/generator", "Dilepton generator", -1);
136
137 //The distribution template:
138 pn_gen->Add("q,parent");
139 pn_gen->Add("p,daughter");
140 pn_gen->Add("n,daughter");
141 pn_gen->Add("dilepton,daughter,primary");
142 pn_gen->SetSampleFunction(flat);
143 //Enable distribution as a generator
144 pn_gen->EnableGenerator();
145 makeDistributionManagerUtil()->Add(pn_gen);
146
147 //repeat this for the pp case:
148 pp_gen =
149 new PInclusiveModel("flat2@p + p_brems_p_p_dilepton/generator", "Dilepton generator", -1);
150
151 //The distribution template:
152 pp_gen->Add("q,parent");
153 pp_gen->Add("p,daughter");
154 pp_gen->Add("p,daughter");
155 pp_gen->Add("dilepton,daughter,primary");
156 pp_gen->SetSampleFunction(flat);
157 //Enable distribution as a generator
158 pp_gen->EnableGenerator();
159 makeDistributionManagerUtil()->Add(pp_gen);
160
161 return kTRUE;
162 } else if (strcmp (command,"kin_max") == 0) {
163 kin_max = value;
164
165 return kTRUE;
166 } else if (strcmp (command,"elastic") == 0) {
167 if (current_author) { //sm
168 pp = pp_sm_elastic;
169 pn = pn_sm_elastic;
170 return makeDistributionManagerUtil()->Enable("bremsstrahlung_sm_elastic");
171 } else {
172 pn->SetMode(0);
173 pp->SetMode(0);
174 }
175 return kTRUE;
176 } else if (strcmp (command,"delta") == 0) {
177 if (current_author) { //sm
178 pp = pp_sm_delta;
179 pn = pn_sm_delta;
180 return makeDistributionManagerUtil()->Enable("bremsstrahlung_sm_delta");
181 } else {
182 pn->SetMode(1);
183 pp->SetMode(1);
184 }
185 return kTRUE;
186 } else if (strcmp (command,"n1520") == 0) {
187 if (current_author) { //sm
188 pp = pp_sm_n1520;
189 pn = pn_sm_n1520;
190 return makeDistributionManagerUtil()->Enable("bremsstrahlung_sm_n1520");
191 } else {
192 Warning("ExecCommand","N1520 option not available for KK");
193 }
194 return kTRUE;
195 } else if (strcmp (command,"sum") == 0) {
196 if (current_author) { //sm
197 pp = pp_sm_total;
198 pn = pn_sm_total;
199 return makeDistributionManagerUtil()->Enable("bremsstrahlung_sm_total");
200 } else {
201 pn->SetMode(2);
202 pp->SetMode(2);
203 }
204 return kTRUE;
205 } else if (strcmp (command,"fsi") == 0) {
206
207 if (pn_fsi) return kTRUE;
208 makeDistributionManagerUtil()->AddSubGroup("fsi", "Final state effects", "root");
209 makeDistributionManagerUtil()->SetGroup("fsi");
210
211 pn_fsi = new PNNFSI("pn_fsi@p + n_brems_p_n_dilepton/nnfsi","PN FSI effects",-1);
212 pn_fsi->EnableWeighting();
213 pn_fsi->SetExpectedWeightMean(-1);
214 makeDistributionManagerUtil()->Add(pn_fsi);
215
216 np_fsi = new PNNFSI("np_fsi@n + p_brems_p_n_dilepton/nnfsi","PN FSI effects",-1);
217 np_fsi->EnableWeighting();
218 np_fsi->SetExpectedWeightMean(-1);
219 makeDistributionManagerUtil()->Add(np_fsi);
220
221 pp_fsi = new PNNFSI("pp_fsi@p + p_brems_p_p_dilepton/nnfsi","PP FSI effects",-1);
222 pp_fsi->EnableWeighting();
223 pp_fsi->SetExpectedWeightMean(-1);
224 pp_fsi->SetA0R0(-7.8098,2.767);
225
226 makeDistributionManagerUtil()->Add(pp_fsi);
227
228 return kTRUE;
229 } else if (strcmp (command,"kaptari") == 0) {
230 current_author = 0;
231
232 if (pp_kk) {
233 pp = pp_kk;
234 pn = pn_kk;
235 makeDistributionManagerUtil()->Enable("bremsstrahlung_kk");
236 return kTRUE;
237 }
238
239 PDistributionManagerUtil *pdmutil = makeDistributionManagerUtil();
240 pdmutil->AddSubGroup("bremsstrahlung_kk", "Bremsstrahlung from KK", "bremsstrahlung");
241 pdmutil->SetGroup("bremsstrahlung_kk");
242
243 pn = new PBremsstrahlung("p + n_brems_p_n_dilepton", "pn Kaptari/Kaempfer bremsstrahlung", -1);
244
245 pn->Add("q,parent");
246 pn->Add("p,grandparent,p1");
247 pn->Add("n,grandparent,p2");
248 pn->Add("p,daughter,p3");
249 pn->Add("n,daughter,p4");
250 pn->Add("dilepton", "daughter");
251
252 PBremsstrahlung *newmodel3 = new PBremsstrahlung("n + p_brems_p_n_dilepton", "pn Kaptari/Kaempfer bremsstrahlung [clone]",-1);
253 newmodel3->Add("dummy,parent");
254
255 pp = new PBremsstrahlung("p + p_brems_p_p_dilepton", "pp Kaptari/Kaempfer bremsstrahlung", -1);
256
257 pp->Add("q,parent");
258 pp->Add("p,grandparent,p1");
259 pp->Add("p,grandparent,p2");
260 pp->Add("p,daughter,p3");
261 pp->Add("p,daughter,p4");
262 pp->Add("dilepton", "daughter");
263
264 //make it known to the Pluto world:
265 pdmutil->Add(pp);
266 pdmutil->Add(pn);
267
268 pdmutil->Add(newmodel3);
269 pp_kk = pp;
270 pn_kk = pn;
271 return kTRUE;
272
273 } else if (strcmp (command,"shyam") == 0) {
274
275 char *path = getenv("PLUTO_DIR");
276 if (!path) {
277 Error ("ExecCommand", "$PLUTO_DIR not set: dat files cannot be located");
278 return kFALSE;
279 }
280
281 current_author = 1;
282
283 if (pp_sm_total) {
284 pp = pp_sm_total;
285 pn = pn_sm_total;
286 makeDistributionManagerUtil()->Enable("bremsstrahlung_sm_total");
287 return kTRUE;
288 }
289
290 //******* Let us read the array from the files:
291 array_pp_sm_125 = new PArray(1); //Dim = 1
292 array_pp_sm_125->Scaling(0.001); //convert to barn
293 array_pp_sm_125->OpenFile(path + TString("/plugins/brems/shyam_pp_1.25GeV.dat"), 0, 5);
294
295
296 array_pn_sm_125 = new PArray(2); //Dim = 2, 3 files
297 array_pn_sm_125->Scaling(0.001); //convert to barn
298 array_pn_sm_125->OpenFile(path + TString("/plugins/brems/shyam_pn_1GeV.dat"), 0, 5, 1);
299 array_pn_sm_125->OpenFile(path + TString("/plugins/brems/shyam_pn_1.25GeV.dat"), 0, 5, 1.25);
300 array_pn_sm_125->OpenFile(path + TString("/plugins/brems/shyam_pn_1.5GeV.dat"), 0, 5, 1.5);
301
302 //******* Modify groups
303 PDistributionManagerUtil * pdmutil = makeDistributionManagerUtil();
304 pdmutil->AddSubGroup("bremsstrahlung_sm", "Bremsstrahlung from Shyam/Mosel", "bremsstrahlung");
305 pdmutil->AddSubGroup("bremsstrahlung_sm_total",
306 "Bremsstrahlung from Shyam/Mosel (total)", "bremsstrahlung_sm");
307 pdmutil->AddSubGroup("bremsstrahlung_sm_elastic",
308 "Bremsstrahlung from Shyam/Mosel (elastic)", "bremsstrahlung_sm");
309 pdmutil->AddSubGroup("bremsstrahlung_sm_delta",
310 "Bremsstrahlung from Shyam/Mosel (delta)", "bremsstrahlung_sm");
311 pdmutil->AddSubGroup("bremsstrahlung_sm_n1520",
312 "Bremsstrahlung from Shyam/Mosel (N*1520)", "bremsstrahlung_sm");
313
314 //******* pp case TOTAL
315 pdmutil->SetGroup("bremsstrahlung_sm_total");
316 pp_sm_total = new PBremsstrahlung("pp_sm_total@p + p_brems_p_p_dilepton",
317 "pp Shyam/Mosel bremsstrahlung (total)", -1);
318
319 pp_sm_total->Add("q,parent");
320 pp_sm_total->Add("p,grandparent,p1");
321 pp_sm_total->Add("p,grandparent,p2");
322 pp_sm_total->Add("p,daughter,p3");
323 pp_sm_total->Add("p,daughter,p4");
324 pp_sm_total->Add("dilepton", "daughter");
325 pp_sm_total->SetFunc(array_pp_sm_125->GetTGraph(0,1)); //1=total
326 pdmutil->Add(pp_sm_total);
327
328 //******* pn case TOTAL
329 pn_sm_total = new PBremsstrahlung("pn_sm_total@p + n_brems_p_n_dilepton",
330 "pn Shyam/Mosel bremsstrahlung (total)", -1);
331
332 pn_sm_total->Add("q,parent");
333 pn_sm_total->Add("p,grandparent,p1");
334 pn_sm_total->Add("n,grandparent,p2");
335 pn_sm_total->Add("p,daughter,p3");
336 pn_sm_total->Add("n,daughter,p4");
337 pn_sm_total->Add("dilepton", "daughter");
338 TGraph2D *gr = array_pn_sm_125->GetTGraph2D(0,1); //1=total
339 pn_sm_total->SetFunc(gr);
340 pdmutil->Add(pn_sm_total);
341
342 PBremsstrahlung *newmodel3 = new PBremsstrahlung("np_sm_total@n + p_brems_p_n_dilepton",
343 "pn Shyam/Mosel bremsstrahlung (total) [clone]", -1);
344 newmodel3->Add("dummy,parent");
345 newmodel3->SetFunc(gr);
346 pdmutil->Add(newmodel3);
347
348 //******* pp case ELASTIC
349 pdmutil->SetGroup("bremsstrahlung_sm_elastic");
350 pp_sm_elastic = new PBremsstrahlung("pp_sm_elastic@p + p_brems_p_p_dilepton",
351 "pp Shyam/Mosel bremsstrahlung (elastic)",-1);
352
353 pp_sm_elastic->Add("p,grandparent,p1");
354 pp_sm_elastic->Add("p,grandparent,p2");
355 pp_sm_elastic->Add("q,parent");
356 pp_sm_elastic->Add("p,daughter,p3");
357 pp_sm_elastic->Add("p,daughter,p4");
358 pp_sm_elastic->Add("dilepton", "daughter");
359 pp_sm_elastic->SetFunc(array_pp_sm_125->GetTGraph(0,2)); //2=elastic
360 pdmutil->Add(pp_sm_elastic);
361
362 //******* pn case ELASTIC
363 pn_sm_elastic = new PBremsstrahlung("pn_sm_elastic@p + n_brems_p_n_dilepton",
364 "pn Shyam/Mosel bremsstrahlung (elastic)", -1);
365
366 pn_sm_elastic->Add("p,grandparent,p1");
367 pn_sm_elastic->Add("n,grandparent,p2");
368 pn_sm_elastic->Add("q,parent");
369 pn_sm_elastic->Add("p,daughter,p3");
370 pn_sm_elastic->Add("n,daughter,p4");
371 pn_sm_elastic->Add("dilepton", "daughter");
372 gr = array_pn_sm_125->GetTGraph2D(0,2); //2=elastic
373 pn_sm_elastic->SetFunc(gr);
374 pdmutil->Add(pn_sm_elastic);
375
376 newmodel3 = new PBremsstrahlung("np_sm_elastic@n + p_brems_p_n_dilepton",
377 "pn Shyam/Mosel bremsstrahlung (elastic) [clone]", -1);
378 newmodel3->Add("dummy,parent");
379 newmodel3->SetFunc(gr);
380 pdmutil->Add(newmodel3);
381
382 //******* pp case DELTA
383 pdmutil->SetGroup("bremsstrahlung_sm_delta");
384 pp_sm_delta = new PBremsstrahlung("pp_sm_delta@p + p_brems_p_p_dilepton",
385 "pp Shyam/Mosel bremsstrahlung (delta)", -1);
386
387 pp_sm_delta->Add("p,grandparent,p1");
388 pp_sm_delta->Add("p,grandparent,p2");
389 pp_sm_delta->Add("q,parent");
390 pp_sm_delta->Add("p,daughter,p3");
391 pp_sm_delta->Add("p,daughter,p4");
392 pp_sm_delta->Add("dilepton", "daughter");
393 pp_sm_delta->SetFunc(array_pp_sm_125->GetTGraph(0,3)); //3=delta
394 pdmutil->Add(pp_sm_delta);
395
396 //******* pn case DELTA
397 pn_sm_delta = new PBremsstrahlung("pn_sm_delta@p + n_brems_p_n_dilepton",
398 "pn Shyam/Mosel bremsstrahlung (delta)", -1);
399
400 pn_sm_delta->Add("p,grandparent,p1");
401 pn_sm_delta->Add("n,grandparent,p2");
402 pn_sm_delta->Add("q,parent");
403 pn_sm_delta->Add("p,daughter,p3");
404 pn_sm_delta->Add("n,daughter,p4");
405 pn_sm_delta->Add("dilepton", "daughter");
406 gr = array_pn_sm_125->GetTGraph2D(0,3); //3=delta
407 pn_sm_delta->SetFunc(gr);
408 pdmutil->Add(pn_sm_delta);
409
410 newmodel3 = new PBremsstrahlung("np_sm_delta@n + p_brems_p_n_dilepton",
411 "pn Shyam/Mosel bremsstrahlung (delta) [clone]", -1);
412 newmodel3->Add("dummy,parent");
413 newmodel3->SetFunc(gr);
414 pdmutil->Add(newmodel3);
415
416 //******* pp case N1520
417 pdmutil->SetGroup("bremsstrahlung_sm_n1520");
418 pp_sm_n1520 = new PBremsstrahlung("pp_sm_n1520@p + p_brems_p_p_dilepton",
419 "pp Shyam/Mosel bremsstrahlung (n1520)", -1);
420
421 pp_sm_n1520->Add("p,grandparent,p1");
422 pp_sm_n1520->Add("p,grandparent,p2");
423 pp_sm_n1520->Add("q,parent");
424 pp_sm_n1520->Add("p,daughter,p3");
425 pp_sm_n1520->Add("p,daughter,p4");
426 pp_sm_n1520->Add("dilepton", "daughter");
427 pp_sm_n1520->SetFunc(array_pp_sm_125->GetTGraph(0,4)); //4=n1520
428 pdmutil->Add(pp_sm_n1520);
429
430 //******* pn case N1520
431 pn_sm_n1520 = new PBremsstrahlung("pn_sm_n1520@p + n_brems_p_n_dilepton",
432 "pn Shyam/Mosel bremsstrahlung (n1520)", -1);
433
434 pn_sm_n1520->Add("p,grandparent,p1");
435 pn_sm_n1520->Add("n,grandparent,p2");
436 pn_sm_n1520->Add("q,parent");
437 pn_sm_n1520->Add("p,daughter,p3");
438 pn_sm_n1520->Add("n,daughter,p4");
439 pn_sm_n1520->Add("dilepton", "daughter");
440 gr = array_pn_sm_125->GetTGraph2D(0,4); //4=n1520
441 pn_sm_n1520->SetFunc(gr);
442 pdmutil->Add(pn_sm_n1520);
443
444 newmodel3 = new PBremsstrahlung("np_sm_n1520@n + p_brems_p_n_dilepton",
445 "pn Shyam/Mosel bremsstrahlung (n1520) [clone]", -1);
446 newmodel3->Add("dummy,parent");
447 newmodel3->SetFunc(gr);
448 pdmutil->Add(newmodel3);
449
450 pp = pp_sm_total;
451 pn = pn_sm_total;
452
453 return kTRUE;
454 }
455
456 return kFALSE;
457}
458
459ClassImp(PBremsstrahlungPlugin)
460
461
462
ClassImp
ClassImp(PBeamLineSimulation) PBeamLineSimulation
Definition PBeamLineSimulation.cc:22
makeDistributionManagerUtil
PDistributionManagerUtil * makeDistributionManagerUtil()
VERSION_GENERATOR_MC
#define VERSION_GENERATOR_MC
Definition PDistribution.h:27
makeStaticData
PStaticData * makeStaticData()
PBremsstrahlungPlugin::ExecCommand
Bool_t ExecCommand(const char *command, Double_t value)
Definition PBremsstrahlungPlugin.cc:114
PBremsstrahlungPlugin::PBremsstrahlungPlugin
PBremsstrahlungPlugin(const Char_t *id, const Char_t *de)
Definition PBremsstrahlungPlugin.cc:37
PBremsstrahlung::SetFunc
void SetFunc(TGraph *gr)
Definition PBremsstrahlung.h:39
PDistributionCollection::RequiresPlugin
Bool_t RequiresPlugin(const char *name)
PDistributionManagerUtil::Add
int Add(PDistribution *dist)
PDistributionManagerUtil::AddSubGroup
Int_t AddSubGroup(const Char_t *id, const Char_t *de, const Char_t *sub)
PDistributionManagerUtil::SetGroup
void SetGroup(const Char_t *id)
PDistributionManagerUtil::Enable
Bool_t Enable(const Char_t *id)
PDistribution::Add
Int_t Add(const Char_t *opt)
PStaticData::AddDecay
int AddDecay(int didx, const char *name, const char *parent, const char *daughters, double br)
PStaticData::GetParticleID
int GetParticleID(const char *id, int warn=1)
flat
TF1 * flat
Definition use_eta_generator.C:10