base_sim.C
1 |
void base_sim(Int_t nEvents = 1000){ |
---|---|
2 |
//---------------------Files-----------------------------------------------
|
3 |
TString outFile= "sim.root";
|
4 |
TString parFile= "par.root";
|
5 |
// ------------------------------------------------------------------------
|
6 |
|
7 |
// ----- Timer --------------------------------------------------------
|
8 |
TStopwatch timer; |
9 |
timer.Start(); |
10 |
// ----
|
11 |
// ----- Create simulation run ----------------------------------------
|
12 |
FairRunSim* run = new FairRunSim();
|
13 |
/** Select transport engine
|
14 |
* TGeant3
|
15 |
* TGeant4
|
16 |
**/
|
17 |
run->SetName("TGeant4"); // Transport engine |
18 |
run->SetOutputFile(outFile.Data()); // Output file
|
19 |
// ------------------------------------------------------------------------
|
20 |
// ----- Runtime database ---------------------------------------------
|
21 |
FairRuntimeDb* rtdb = run->GetRuntimeDb(); |
22 |
// ------------------------------------------------------------------------
|
23 |
// ----- Create media -------------------------------------------------
|
24 |
run->SetMaterials("media_my.geo"); // Materials |
25 |
// ------------------------------------------------------------------------
|
26 |
|
27 |
//-------- Set MC event header --------------------------------------------
|
28 |
ERMCEventHeader* header = new ERMCEventHeader();
|
29 |
run->SetMCEventHeader(header); |
30 |
//---------------------------------
|
31 |
// ----- Create detectors ----------------------------------------------
|
32 |
FairModule* cave= new ERCave("CAVE"); |
33 |
cave->SetGeometryFileName("cave.geo");
|
34 |
run->AddModule(cave); |
35 |
|
36 |
ERCollimator* collimator = new ERCollimator();
|
37 |
collimator->SetGeometryFileName("collimator.mygeo.root");
|
38 |
run->AddModule(collimator); |
39 |
|
40 |
|
41 |
/*
|
42 |
FairModule* target = new ERTarget("Target", kTRUE,1);
|
43 |
target->SetGeometryFileName("target.mygeo.root");
|
44 |
run->AddModule(target);
|
45 |
*/
|
46 |
|
47 |
// Det definition
|
48 |
/* Select verbosity level
|
49 |
* 0 - only standard logs
|
50 |
* 1 - Print points after each event
|
51 |
*/
|
52 |
Int_t verbose = 0;
|
53 |
|
54 |
ERDetector* detector= new ERDetector("TestDetector", kTRUE,verbose); |
55 |
detector->SetGeometryFileName("base.geo.root");
|
56 |
//detector->SetGeometryFileName("base.geo_2.root");
|
57 |
detector->AddSensetive("gas");
|
58 |
detector->AddSensetive("plastic");
|
59 |
run->AddModule(detector); |
60 |
|
61 |
ERDetector* detector1= new ERDetector("TestDetector", kTRUE,verbose); |
62 |
detector1->SetGeometryFileName("base.geo_2.root");
|
63 |
detector1->AddSensetive("gas1");
|
64 |
detector1->AddSensetive("plastic1");
|
65 |
run->AddModule(detector1); |
66 |
|
67 |
ERDetector* detector2= new ERDetector("TestDetector", kTRUE,verbose); |
68 |
detector2->SetGeometryFileName("base.geo_4.root");
|
69 |
detector2->AddSensetive("gas2");
|
70 |
detector2->AddSensetive("plastic2");
|
71 |
run->AddModule(detector2); |
72 |
// FairModule* target = new ERTarget("BeamDetTarget", kTRUE, 1);
|
73 |
//target->SetGeometryFileName("target.h.geo.root");
|
74 |
//run->AddModule(target);
|
75 |
//
|
76 |
//--------------------------------------------------------
|
77 |
// ----- Create PrimaryGenerator --------------------------------------
|
78 |
FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
|
79 |
|
80 |
//Ion 15N
|
81 |
Int_t A = 15;
|
82 |
Int_t Z = 7;
|
83 |
Int_t Q = 3;
|
84 |
|
85 |
ERIonMixGenerator* generator = new ERIonMixGenerator("15N", Z, A, Q, 1); |
86 |
Double32_t kin_energy = 0.043; //GeV |
87 |
//generator->SetKinESigma(kin_energy, 0);
|
88 |
//generator->SetPSigma(6.7835, 6.7835*0.003);
|
89 |
generator->SetKinESigma(kin_energy, 0);
|
90 |
|
91 |
generator->SpreadingOnTarget(); |
92 |
|
93 |
Double32_t theta = 0;
|
94 |
Double32_t sigmaTheta = 0.004*TMath::RadToDeg(); |
95 |
generator->SetThetaSigma(theta, sigmaTheta); |
96 |
|
97 |
generator->SetPhiRange(0, 360); |
98 |
|
99 |
Double32_t distanceToTarget = 200;
|
100 |
Double32_t sigmaOnTarget = 0.5; |
101 |
generator->SetSigmaXYZ(0, 0, -distanceToTarget, sigmaOnTarget, sigmaOnTarget); |
102 |
generator->SetBoxXYZ(-0.4,-0.4,0.4,0.4, -distanceToTarget); |
103 |
|
104 |
//generator->AddBackgroundIon("26P", 15, 26, 15, 0.25);
|
105 |
//generator->AddBackgroundIon("26S", 16, 26, 16, 0.25);
|
106 |
//generator->AddBackgroundIon("24Si", 14, 24, 14, 0.25);
|
107 |
|
108 |
primGen->AddGenerator(generator); |
109 |
/*
|
110 |
Int_t pdgId = 2212; //
|
111 |
Double32_t theta1 = 0.; // polar angle distribution
|
112 |
Double32_t theta2 = 0.;
|
113 |
//Double32_t momentum = 0.05; //GeV
|
114 |
Double32_t kin_energy = .05; //GeV
|
115 |
Double_t mass = TDatabasePDG::Instance()->GetParticle(pdgId)->Mass();
|
116 |
Double32_t momentum = TMath::Sqrt(kin_energy*kin_energy + 2.*kin_energy*mass); //GeV
|
117 |
|
118 |
FairBoxGenerator* boxGen = new FairBoxGenerator(pdgId, 1);
|
119 |
boxGen->SetThetaRange(theta1, theta2);
|
120 |
boxGen->SetPRange(momentum, momentum);
|
121 |
boxGen->SetPhiRange(0,360);
|
122 |
boxGen->SetXYZ(0.,0., -100.);
|
123 |
|
124 |
primGen->AddGenerator(boxGen);
|
125 |
*/
|
126 |
run->SetGenerator(primGen); |
127 |
// ------------------------------------------------------------------------
|
128 |
//-------Set visualisation flag to true------------------------------------
|
129 |
run->SetStoreTraj(kTRUE); |
130 |
|
131 |
//-------Set LOG verbosity -----------------------------------------------
|
132 |
FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
|
133 |
// —– Initialize simulation run ————————————
|
134 |
run->Init(); |
135 |
Int_t nSteps = -15000;
|
136 |
|
137 |
// —– Runtime database ———————————————
|
138 |
Bool_t kParameterMerged = kTRUE; |
139 |
FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
|
140 |
parOut->open(parFile.Data()); |
141 |
rtdb->setOutput(parOut); |
142 |
rtdb->saveOutput(); |
143 |
rtdb->print(); //
|
144 |
// ----- Run simulation ------------------------------------------------
|
145 |
run->Run(nEvents); |
146 |
|
147 |
// ----- Finish -------------------------------------------------------
|
148 |
timer.Stop(); |
149 |
Double_t rtime = timer.RealTime(); |
150 |
Double_t ctime = timer.CpuTime(); |
151 |
cout << endl << endl; |
152 |
cout << "Macro finished succesfully." << endl;
|
153 |
cout << "Output file is sim.root" << endl;
|
154 |
cout << "Parameter file is par.root" << endl;
|
155 |
cout << "Real time " << rtime << " s, CPU time " << ctime |
156 |
<< "s" << endl << endl;
|
157 |
// cout << "Energy " << momentum << "; mass " << mass << endl;
|
158 |
} |