Pluto v6.3.0
Loading...
Searching...
No Matches
PBeamLineSimulation.cc
Go to the documentation of this file.
1
2//
3// Beam line transport
4//
5// Author: T. Hennino (core), J. Biernat
6// Pluto implementation by IF
7// Written: 14.05.2013
8//
10
11
12using namespace std;
13#include <sstream>
14#include <iostream>
15#include <iomanip>
16#include <fstream>
17
18#include "PBeamLineSimulation.h"
19#include "PDefines.h"
20
21
26
27PBeamLineSimulation::~PBeamLineSimulation() {
28 if (mybeam) delete (mybeam);
29 if (mytarget) delete (mytarget);
30};
31
32
33PBeamLineSimulation::PBeamLineSimulation(Char_t *id, Char_t *de) :
34 PBeamSmearing(id, de) {
35 beam = NULL;
36 mybeam = NULL;
37 target = NULL;
38 mytarget = NULL;
39 parent = NULL;
40
41 beam_counter = 0;
42 detector_counter = 0;
43 num_target = -1;
44 global_p = 0.;
45 beamtube_size_x = beamtube_size_y = -1.;
46
47 ii = 0;
48 ind1 = 0;
49 ind23 = 0;
50};
51
52PDistribution* PBeamLineSimulation::Clone(const char *) const {
53 return new PBeamLineSimulation((const PBeamLineSimulation &)* this);
54};
55
56void PBeamLineSimulation::MakeVars() {
57 if (!projector) projector = new PProjector();
58
59 in_xi = makeStaticData()->GetBatchValue("_beam_x");
60 in_yi = makeStaticData()->GetBatchValue("_beam_y");
61 in_dp = makeStaticData()->GetBatchValue("_beam_dp");
62 in_xpi = makeStaticData()->GetBatchValue("_beam_px");
63 in_ypi = makeStaticData()->GetBatchValue("_beam_py");
64
65 filter = makeStaticData()->GetBatchValue("#_beamline_filter");
66};
67
68Bool_t PBeamLineSimulation::AddEquation(char *command) {
69 MakeVars();
70 return projector->AddCommand(command);
71};
72
73Bool_t PBeamLineSimulation::AddHistogram(TH2 *histo, const char *command) {
74 MakeVars();
75 return projector->AddHistogram(histo, command, 0);
76};
77
78Int_t PBeamLineSimulation::AddBeamParticle(char *name) {
79
80 //create a clone of the string
81 char * dummy = new char[strlen(name)+1];
82 strcpy(dummy, name);
83
84 //create key in data base:
85 Int_t key1 = makeStaticData()->
86 MakeDirectoryEntry("batch_objects", NBATCH_NAME, LBATCH_NAME, dummy);
87
88 //check if we use this key already
89 int used_id = -1;
90 if (makeGlobalBulk()->GetSizeBranches() && makeGlobalBulk()->GetKeysBranches()) {
91 for (int br=0; br<*(makeGlobalBulk()->GetSizeBranches()); br++) {
92 if ((makeGlobalBulk()->GetKeysBranches())[br] == key1) used_id = br;
93 }
94 } else {
95 Error("AddBeamParticle","No access to branches (no PReaction?)");
96 }
97 if (used_id == -1) {
98 if (*(makeGlobalBulk()->GetSizeBranches()) != MAX_NUM_BRANCHES) {
99 (makeGlobalBulk()->GetKeysBranches()[*(makeGlobalBulk()->GetSizeBranches())]) = key1;
100 used_id = *(makeGlobalBulk()->GetSizeBranches());
101 (*(makeGlobalBulk()->GetSizeBranches()))++;
102 } else {
103 Error("AddDetector","MAX_NUM_BRANCHES reached");
104 return -1;
105 }
106 }
107
108 //generate beam particles and put them to the stack
109
110 if (beam_counter == MAX_NUM_BEAM_PARTICLES) {
111 Error("AddBeamParticle", "MAX_NUM_BEAM_PARTICLES reached");
112 return -1;
113 }
114
115 beam_particles[beam_counter] = new PParticle();
116 beam_id[beam_counter] = used_id;
117 beam_counter++;
118
119 return beam_counter;
120};
121
122Int_t PBeamLineSimulation::AddDetector(char *name, Double_t distance) {
123 //adds a detector "name" at place "distance" [in mm]
124 //(distance must be negative if placed before the target)
125
126 if (detector_counter == MAX_NUM_BEAM_PARTICLES) {
127 Error("AddDetector","MAX_NUM_BEAM_PARTICLES reached");
128 return -1;
129 }
130
131 detector_name[detector_counter] = name;
132 detector_distance[detector_counter] = distance;
133 detector_counter++;
134
135 return detector_counter - 1;
136}
137
138Bool_t PBeamLineSimulation::TargetIsElement(Int_t n) {
139 if (n<1 || n>=nelem) return kFALSE;
140 num_target = n;
141 for (int i=1; i<nelem; i++)
142 along[i] -= along[n];
143
144 return kTRUE;
145}
146
147Bool_t PBeamLineSimulation::InitBeamLine(char *filename) {
148 //opens the ASCII description of the beam line
149
150 nelem = 1;
151
152 std::ifstream plik1;
153 plik1.open(filename);
154 string str;
155
156 Bool_t readfile = kTRUE;
157
158 while (readfile) {
159
160 str.clear();
161 int i = nelem;
162
163 plik1>>along[i]>>txt2>>aa>>sig11>>txt5;
164 plik1>>aa>>sig22>>txt5>>r12;
165 plik1>>aa>>sig33>>txt5>>r13>>r23;
166 plik1>>aa>>sig44>>txt5>>r14>>r24>>r34;
167 plik1>>aa>>sig55>>txt5>>r15>>r25>>r35>>r45;
168 plik1>>aa>>sig66>>txt5>>r16>>r26>>r36>>r46>>r56;
169
170 getline(plik1,str,'\n');
171 getline(plik1,str,'\n');
172
173 plik1>>T11[i]>>T12[i]>>T13[i]>>T14[i]>>aa>>T16[i];
174 plik1>>T21[i]>>T22[i]>>T23[i]>>T24[i]>>aa>>T26[i];
175 plik1>>T31[i]>>T32[i]>>T33[i]>>T34[i]>>aa>>T36[i]>>std::fixed;
176 // std::cout<<T31[i]<<endl;
177 plik1>>T41[i]>>T42[i]>>T43[i]>>T44[i]>>aa>>T46[i];
178 plik1>>aa>>aa>>aa>>aa>>aa>>aa;
179 plik1>>aa>>aa>>aa>>aa>>aa>>aa;
180 getline(plik1,str);
181
182 //C
183 //C lecture second ordre T1ij
184 //C
185 getline(plik1,str);
186 plik1>>ind1>>ind23>>T111[i];
187 plik1>>ind1>>ind23>>T112[i]>>ind1>>ind23>>T122[i];
188 plik1>>ind1>>ind23>>T113[i]>>ind1>>ind23>>T123[i]>>ind1>>ind23>>T133[i];
189 plik1>>ind1>>ind23>>T114[i]>>ind1>>ind23>>T124[i]>>ind1>>ind23>>T134[i]>>ind1>>ind23>>T144[i];
190 plik1>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa;
191 plik1>>ind1>>ind23>>T116[i]>>ind1>>ind23>>T126[i]>>ind1>>ind23>>T136[i]>>ind1>>ind23>>T146[i]>>ind1>>ind23>>aa>>ind1>>ind23>>T166[i];
192 getline(plik1,str);
193 //plik2<<".."<<i<<".."<<T166[i]<<endl;
194
195 //C
196 //C lecture second ordre T2ij
197 //C
198 plik1>>ii>>ii>>T211[i];
199 plik1>>ii>>ii>>T212[i]>>ii>>ii>>T222[i];
200 plik1>>ii>>ii>>T213[i]>>ii>>ii>>T223[i]>>ii>>ii>>T233[i];
201 plik1>>ii>>ii>>T214[i]>>ii>>ii>>T224[i]>>ii>>ii>>T234[i]>>ii>>ii>>T244[i];
202 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
203 plik1>>ii>>ii>>T216[i]>>ii>>ii>>T226[i]>>ii>>ii>>T236[i]>>ii>>ii>>T246[i]>>ii>>ii>>aa>>ii>>ii>>T266[i];
204 getline(plik1,str);
205 //plik2<<".."<<i<<".."<<T266[i]<<endl;
206 //rT266 = T266[i];
207
208 //C
209 //C lecture second ordre T3ij
210 //C
211 plik1>>ind1>>ind23>>T311[i];
212 plik1>>ind1>>ind23>>T312[i]>>ind1>>ind23>>T322[i];
213 plik1>>ind1>>ind23>>T313[i]>>ind1>>ind23>>T323[i]>>ind1>>ind23>>T333[i];
214 plik1>>ind1>>ind23>>T314[i]>>ind1>>ind23>>T324[i]>>ind1>>ind23>>T334[i]>>ind1>>ind23>>T344[i];
215 plik1>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa>>ind1>>ind23>>aa;
216 plik1>>ind1>>ind23>>T316[i]>>ind1>>ind23>>T326[i]>>ind1>>ind23>>T336[i]>>ind1>>ind23>>T346[i]>>ind1>>ind23>>aa>>ind1>>ind23>>T366[i];
217 getline(plik1,str);
218 //plik2<<".."<<i<<".."<<T366[i]<<endl;
219 //rT366 = T366[i];
220
221 //C
222 //C lecture second ordre T4ij
223 //C
224
225 plik1>>ii>>ii>>T411[i];
226 plik1>>ii>>ii>>T412[i]>>ii>>ii>>T422[i];
227 plik1>>ii>>ii>>T413[i]>>ii>>ii>>T423[i]>>ii>>ii>>T433[i];
228 plik1>>ii>>ii>>T414[i]>>ii>>ii>>T424[i]>>ii>>ii>>T434[i]>>ii>>ii>>T444[i];
229 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
230 plik1>>ii>>ii>>T416[i]>>ii>>ii>>T426[i]>>ii>>ii>>T436[i]>>ii>>ii>>T446[i]>>ii>>ii>>aa>>ii>>ii>>T466[i];
231 getline(plik1,str);
232 //plik2<<".."<<i<<".."<<T466[i]<<endl;
233 //rT466 = T466[i];
234 //C
235 //C lecture second ordre T5ij
236 //C
237
238 plik1>>ii>>ii>>aa;
239 plik1>>ii>>ii>>aa>>ii>>ii>>aa;
240 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
241 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
242 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
243 plik1>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa>>ii>>ii>>aa;
244 getline(plik1,str);
245 getline(plik1,str,'\n');
246
247 if (!plik1.eof())
248 nelem++;
249 else
250 readfile = kFALSE;
251
252 } //end nelem
253
254 plik1.close();
255
256 Info("InitBeamLine", "Read %i elements from file %s", nelem-1, filename);
257
258 return kTRUE;
259}
260
261
262Bool_t PBeamLineSimulation::Init(void) {
263 //Read beam, target and parent
264
265 beam = GetParticle("beam");
266 target = GetParticle("target");
267 parent = GetParticle("parent");
268
269 if (!beam || !target) {
270 Error("Init", "beam or target not found");
271 return kFALSE;
272 }
273
274 if (!mybeam)
275 mybeam = new PParticle(beam);
276 else
277 *mybeam = *beam;
278 if (!mytarget)
279 mytarget= new PParticle(target);
280 else
281 *mytarget = target;
282
283 if (!parent){
284 Error("Init","Parent not found");
285 return kFALSE;
286 }
287
288 for (int i=0; i<detector_counter; i++) {
289 AddBeamParticle(detector_name[i]);
290 }
291
292 if (!projector) {
293 Error("Init", "PProjector not found, you have not specified the initial conditions");
294 return kFALSE;
295 }
296
297 return kTRUE;
298}
299
300Bool_t PBeamLineSimulation::Prepare(void) {
301 //We use the prepare function for smearing since it might affect
302 //the mass and momentum sampling done in the next steps
303
304 *filter = 1;
305
306 //restore saved particles
307 *beam = *mybeam;
308 *target = *mytarget;
309
310 PParticle *real_beam = target;
311
312 if (beam->P() > target->P()) {
313 //change "beam"
314 real_beam = beam;
315 }
316
317 //the global event vertex cannot be set at this point.
318 //this would require a major change in PReaction
319
320 //change momentum of "real_beam" here.....
321
322 //beam offset (in cm) at prod target in x
323 //beam offset (in cm) at prod target in y
324 Int_t num = 0;
325 projector->Modify(NULL, NULL, &num, 0);
326
327 xi = *in_xi/10.; //convert from Pluto standard (mm) to cm
328 yi = *in_xi/10.; //convert from Pluto standard (mm) to cm
329 xpi = *in_xpi*1000.; //convert from rad to mrad
330 ypi = *in_ypi*1000.; //convert from rad to mrad
331 dp = *in_dp * 100.; //convert from dp/p to percent (as discussed with T.H. at CM Prague)
332 //dp = 0.08;
333
334 //now loop over the number of elements
335 for(Int_t k = 1; k < nelem; k++){
336
337 z1 = T11[k]*xi;
338 z2 = T12[k]*xpi;
339 z3 = T13[k]*yi;
340 z4 = T14[k]*ypi;
341 z5 = T16[k]*dp;
342
343 xf = z1+z2+z3+z4+z5;
344 yf = T31[k]*xi+T32[k]*xpi+T33[k]*yi+T34[k]*ypi+T36[k]*dp; //yf=T31(k)*xi+T32(k)*xpi+T33(k)*yi+T34(k)*ypi+T36(k)*dp
345 // std::cout<<"T31 = .. "<<T31[k]<<"..T32 = .. "<<T32[k]<<"..T33 = .. "<<T33[k]<<"..T34 = .. "<<T34[k]<<"..T36 = .. "<<T36[k]<<std::fixed<<endl;
346 xpf = T21[k]*xi+T22[k]*xpi+T23[k]*yi+T24[k]*ypi+T26[k]*dp;
347 ypf = T41[k]*xi+T42[k]*xpi+T43[k]*yi+T44[k]*ypi+T46[k]*dp;
348 xf1 = xf;
349 yf1 = yf;
350 xpf1 = xpf;
351 ypf1 = ypf;
352
353 xf = xf1+T111[k]*xi*xi+T112[k]*xi*xpi+T113[k]*xi*yi;
354 xf = xf+T114[k]*xi*ypi+T116[k]*xi*dp;
355 xf = xf+T122[k]*xpi*xpi+T123[k]*xpi*yi+T124[k]*xpi*ypi;
356 xf = xf+T126[k]*xpi*dp;
357 xf = xf+T133[k]*yi*yi+T134[k]*yi*ypi+T136[k]*yi*dp;
358 xf = xf+T144[k]*ypi*ypi+T146[k]*ypi*dp;
359 xf = xf+T166[k]*dp*dp;
360
361 // xf=xf1+T126[k]*xpi*dp+T166[k]*dp*dp;
362 // xf=xf1+T126[k]*xpi*dp;
363
364 xpf = xpf1+T211[k]*xi*xi+T212[k]*xi*xpi+T213[k]*xi*yi;
365 xpf = xpf+T214[k]*xi*ypi+T216[k]*xi*dp;
366 xpf = xpf+T222[k]*xpi*xpi+T223[k]*xpi*yi+T224[k]*xpi*ypi;
367 xpf = xpf+T226[k]*xpi*dp;
368 xpf = xpf+T233[k]*yi*yi+T234[k]*yi*ypi+T236[k]*yi*dp;
369 xpf = xpf+T244[k]*ypi*ypi+T246[k]*ypi*dp;
370 xpf = xpf+T266[k]*dp*dp;
371
372 yf = yf1+T311[k]*xi*xi+T312[k]*xi*xpi+T313[k]*xi*yi;
373 yf = yf+T314[k]*xi*ypi+T316[k]*xi*dp;
374 yf = yf+T322[k]*xpi*xpi+T323[k]*xpi*yi+T324[k]*xpi*ypi;
375 yf = yf+T326[k]*xpi*dp;
376 yf = yf+T333[k]*yi*yi+T334[k]*yi*ypi+T336[k]*yi*dp;
377 yf = yf+T344[k]*ypi*ypi+T346[k]*ypi*dp;
378 yf = yf+T366[k]*dp*dp;
379
380 ypf = ypf1+T411[k]*xi*xi+T412[k]*xi*xpi+T413[k]*xi*yi;
381 ypf = ypf+T414[k]*xi*ypi+T416[k]*xi*dp;
382 ypf = ypf+T422[k]*xpi*xpi+T423[k]*xpi*yi+T424[k]*xpi*ypi;
383 ypf = ypf+T426[k]*xpi*dp;
384 ypf = ypf+T433[k]*yi*yi+T434[k]*yi*ypi+T436[k]*yi*dp;
385 ypf = ypf+T444[k]*ypi*ypi+T446[k]*ypi*dp;
386 ypf = ypf+T466[k]*dp*dp;
387 //C
388 //C yf=yf1+T346[k]*ypi*dp+T366[k]*dp*dp
389 //C yf=yf1+T346[k]*ypi*dp
390 // rxf =
391 xf2 = xf;
392 xpf2 = xpf;
393 yf2 = yf;
394 ypf2 = ypf;
395
396 out_xi[k] = xf2*10.; //convert back from cm. to mm
397 out_yi[k] = yf2*10.; //convert back from cm. to mm
398 out_dp[k] = dp/100.; //convert back from % to fraction
399 out_xpi[k] = xpf2/1000.; //convert back from mrad to rad
400 out_ypi[k] = ypf2/1000.; //convert back from mrad to rad
401
402 //cout << out_xi[k] << ":" << out_yi[k] << ":" << out_xpi[k] << ":" << out_ypi[k] << endl;
403
404 }
405
406 Double_t p = global_p * (1. + out_dp[num_target]);
407 Double_t p_z = p / sqrt(1.+ tan(out_xpi[num_target])*tan(out_xpi[num_target]) +
408 tan(out_ypi[num_target])*tan(out_ypi[num_target]));
409 Double_t p_x = p_z * tan(out_xpi[num_target]);
410 Double_t p_y = p_z * tan(out_ypi[num_target]);
411
412 real_beam->SetXYZM(p_x, p_y, p_z, makeStaticData()->GetParticleMass(real_beam->ID()));
413
414 //real_beam->Print();
415
416 parent->Reconstruct();
417
418 //parent->Print();
419
420 //change particles at detector places
421 for (int i=0; i<beam_counter; i++) {
422
423 for (int k=1; k<nelem-1; k++) {
424 if (along[k]*1000. <= detector_distance[i] && detector_distance[i] < along[k+1]*1000.) {
425
426 Double_t frac = (detector_distance[i] - along[k]*1000.) / (along[k+1]*1000. - along[k]*1000.);
427
428 Double_t my_xi = out_xi[k] + frac*(out_xi[k+1] - out_xi[k]);
429 Double_t my_yi = out_yi[k] + frac*(out_yi[k+1] - out_yi[k]);
430 Double_t my_dp = out_dp[k] + frac*(out_dp[k+1] - out_dp[k]);
431 Double_t my_xpi = out_xpi[k] + frac*(out_xpi[k+1] - out_xpi[k]);
432 Double_t my_ypi = out_ypi[k] + frac*(out_ypi[k+1] - out_ypi[k]);
433 p = global_p * (1. + my_dp);
434
435 //cout << p << ":" << my_xpi << ":" << my_ypi << endl;
436
437 p_z = p / sqrt(1.+ tan(my_xpi)*tan(my_xpi) + tan(my_ypi)*tan(my_ypi));
438 p_x = p_z * tan(my_xpi);
439 p_y = p_z * tan(my_ypi);
440 *(beam_particles[i]) = *(real_beam);
441 beam_particles[i]->SetVertex(my_xi, my_yi, detector_distance[i]);
442 beam_particles[i]->SetPx(p_x);
443 beam_particles[i]->SetPy(p_y);
444 beam_particles[i]->SetPz(p_z);
445
446 //cout << p_x << ":" << p_y << ":" << p_z << endl;
447
448 }
449
450 }
451 }
452
453 for (int k=1; k<nelem-1; k++) {
454 if (beamtube_size_x > 0 && out_xi[k] > beamtube_size_x) {
455 *filter = 0;
456 }
457 if (beamtube_size_y > 0 && out_yi[k] > beamtube_size_y) {
458 *filter = 0;
459 }
460 }
461
462 for (int i=0; i<beam_counter; i++)
463 (*(makeGlobalBulk()->GetBranchNum(beam_id[i]))) = 0;
464
465
466 for (int i=0; i<beam_counter; i++) {
467 *((makeGlobalBulk()->GetBranchArray(beam_id[i]))
468 [*(makeGlobalBulk()->GetBranchNum(beam_id[i]))] )
469 = *(beam_particles[i]);
470 (*(makeGlobalBulk()->GetBranchNum(beam_id[i])))++;
471 }
472
473 return kTRUE;
474}
475
ClassImp
ClassImp(PBeamLineSimulation) PBeamLineSimulation
Definition PBeamLineSimulation.cc:22
MAX_NUM_BEAM_PARTICLES
#define MAX_NUM_BEAM_PARTICLES
Definition PBeamLineSimulation.h:9
makeGlobalBulk
PBulkInterface * makeGlobalBulk()
MAX_NUM_BRANCHES
#define MAX_NUM_BRANCHES
Definition PDefines.h:19
LBATCH_NAME
#define LBATCH_NAME
Definition PDefines.h:14
NBATCH_NAME
#define NBATCH_NAME
Definition PDefines.h:13
makeStaticData
PStaticData * makeStaticData()
PBeamLineSimulation::AddDetector
Int_t AddDetector(char *name, Double_t distance)
Definition PBeamLineSimulation.cc:122
PBeamLineSimulation::AddEquation
Bool_t AddEquation(char *command)
Definition PBeamLineSimulation.cc:68
PBeamLineSimulation::TargetIsElement
Bool_t TargetIsElement(Int_t n)
Definition PBeamLineSimulation.cc:138
PBeamLineSimulation::AddHistogram
Bool_t AddHistogram(TH2 *histo, const char *command="")
Definition PBeamLineSimulation.cc:73
PBeamLineSimulation::Clone
PDistribution * Clone(const char *delme=NULL) const
Definition PBeamLineSimulation.cc:52
PBeamLineSimulation::InitBeamLine
Bool_t InitBeamLine(char *filename)
Definition PBeamLineSimulation.cc:147
PBulkInterface::GetSizeBranches
Int_t * GetSizeBranches()
Definition PBulkInterface.h:60
PBulkInterface::GetBranchArray
PParticle ** GetBranchArray(Int_t i)
Definition PBulkInterface.h:63
PBulkInterface::GetKeysBranches
Int_t * GetKeysBranches()
Definition PBulkInterface.h:61
PBulkInterface::GetBranchNum
Int_t * GetBranchNum(Int_t i)
Definition PBulkInterface.h:62
PDistribution::GetParticle
PParticle * GetParticle(const Char_t *pflag=NULL)
PParticle::ID
Int_t ID() const
Definition PParticle.h:80
PParticle::Reconstruct
void Reconstruct(void)
PStaticData::GetBatchValue
Double_t * GetBatchValue(const char *name, Int_t make_val=1)
p
PParticle * p
Definition dp_mult.C:24
filename
TString filename("eta_dalitz_a")
n
Int_t n
Definition etap+omeg.C:21
histo
TH1F * histo
Definition test_phi2170_ingo.C:16