|
1
|
void exp1803_full(Int_t nEvents = 100) {
|
|
2
|
// --------------- Telescope T1 -------------------------------------------
|
|
3
|
Double_t T1Dl = 0.5; // [cm]
|
|
4
|
Double_t T1PosZ = 10.; // [cm]
|
|
5
|
Double_t T1D1Thick = 0.002; // [cm]
|
|
6
|
Double_t T1D2Thick = 0.1; // [cm]
|
|
7
|
Double_t T1Side = 6.2; // [cm] this parameter should coincide with SD1 side length from /db/QTelescope/QTelescopeParts.xml
|
|
8
|
Double_t T1Aperture = 3.1; // [cm]
|
|
9
|
// --------------- Telescope D1 -------------------------------------------
|
|
10
|
Double_t D1PosZ = 20.; // [cm]
|
|
11
|
Double_t D1Thick = 0.03; // [cm]
|
|
12
|
// --------------- BeamDet ------------------------------------------------
|
|
13
|
Double_t BeamDetLToF = 1500.; // [cm]
|
|
14
|
Double_t BeamDetPosZToF = -50; // [cm]
|
|
15
|
Double_t BeamDetLMWPC = 32.; // [cm]
|
|
16
|
Double_t BeamDetPosZMWPC = -8; // [cm]
|
|
17
|
// --------------- Beam start position ------------------------------------
|
|
18
|
Double_t beamStartPosition = -1600; // [cm]
|
|
19
|
// --------------- Target -------------------------------------------------
|
|
20
|
Double_t targetH2Thickness = 0.4; // [cm] this parameter should coincide with target H2 thickness in /macro/geo/create_GadastEXP1803_geo.C
|
|
21
|
//---------------------Files-----------------------------------------------
|
|
22
|
TString outFile= "full.root";
|
|
23
|
TString parFile= "par.root";
|
|
24
|
TString workDirPath = gSystem->Getenv("VMCWORKDIR");
|
|
25
|
TString paramFileQTelescope = workDirPath
|
|
26
|
+ "/db/QTelescope/QTelescopeParts.xml";
|
|
27
|
TString paramFileBeamDet = workDirPath
|
|
28
|
+ "/db/BeamDet/BeamDetParts.xml";
|
|
29
|
TString targetGeoFileName = workDirPath + "/geometry/target.h2.geo.root";
|
|
30
|
TString gadastGeoFileName = workDirPath + "/geometry/partOfGadast.v1.geo.root";
|
|
31
|
TString ndGeoFileName = workDirPath + "/geometry/ND.geo.root";
|
|
32
|
TString magnetGeoFileName = workDirPath + "/geometry/magnet.geo.root";
|
|
33
|
// ------------------------------------------------------------------------
|
|
34
|
|
|
35
|
// ----- Timer --------------------------------------------------------
|
|
36
|
TStopwatch timer;
|
|
37
|
timer.Start();
|
|
38
|
// ------------------------------------------------------------------------
|
|
39
|
|
|
40
|
// ----- Create simulation run ----------------------------------------
|
|
41
|
ERRunSim* run = new ERRunSim();
|
|
42
|
/** Select transport engine
|
|
43
|
* TGeant3
|
|
44
|
* TGeant4
|
|
45
|
**/
|
|
46
|
run->SetName("TGeant4"); // Transport engine
|
|
47
|
run->SetOutputFile(outFile.Data()); // Output file
|
|
48
|
// ------------------------------------------------------------------------
|
|
49
|
// ----- Runtime database ---------------------------------------------
|
|
50
|
FairRuntimeDb* rtdb = run->GetRuntimeDb();
|
|
51
|
//-------- Set MC event header --------------------------------------------
|
|
52
|
ERDecayMCEventHeader* decayMCheader = new ERDecayMCEventHeader();
|
|
53
|
run->SetMCEventHeader(decayMCheader);
|
|
54
|
// ----- Create media -------------------------------------------------
|
|
55
|
run->SetMaterials("media.geo"); // Materials
|
|
56
|
//-------- Set MC event header --------------------------------------------
|
|
57
|
// ERMCEventHeader* header = new ERMCEventHeader();
|
|
58
|
// run->SetMCEventHeader(header);
|
|
59
|
// ----- Create detectors ----------------------------------------------
|
|
60
|
FairModule* cave= new ERCave("CAVE");
|
|
61
|
cave->SetGeometryFileName("cave.geo");
|
|
62
|
run->AddModule(cave);
|
|
63
|
|
|
64
|
Int_t verbose = 1;
|
|
65
|
// ----- BeamDet Setup ---------------------------------------------------
|
|
66
|
ERBeamDetSetup* setupBeamDet = ERBeamDetSetup::Instance();
|
|
67
|
setupBeamDet->SetXmlParametersFile(paramFileBeamDet);
|
|
68
|
|
|
69
|
// ----- BeamDet parameters ----------------------------------------------
|
|
70
|
setupBeamDet->AddToF("ToF1", BeamDetPosZToF - BeamDetLToF); //
|
|
71
|
setupBeamDet->AddToF("ToF1", BeamDetPosZToF); // BeamDet parts should be added in ascending order
|
|
72
|
setupBeamDet->AddMWPC("MWPC1", BeamDetPosZMWPC - BeamDetLMWPC); // of Z-coordinate of part.
|
|
73
|
setupBeamDet->AddMWPC("MWPC1", BeamDetPosZMWPC); //
|
|
74
|
//setupBeamDet->SetSensitiveTarget();
|
|
75
|
|
|
76
|
// ----- Create target -------------------------------------------------
|
|
77
|
FairModule* target = new ERTarget("targetH2", kTRUE, 1);
|
|
78
|
target->SetGeometryFileName(targetGeoFileName);
|
|
79
|
run->AddModule(target);
|
|
80
|
|
|
81
|
// ----- Create Part of Gadast ------------------------------------------
|
|
82
|
ERGadast* gadast = new ERGadast("PartofGadast", kTRUE, 1);
|
|
83
|
gadast->SetGeometryFileName(gadastGeoFileName);
|
|
84
|
run->AddModule(gadast);
|
|
85
|
|
|
86
|
// ----- Create Stilbene wall -------------------------------------------
|
|
87
|
ERND* neutronDetector = new ERND("StilbeneWall", kTRUE, 1);
|
|
88
|
neutronDetector->SetGeometryFileName(ndGeoFileName);
|
|
89
|
run->AddModule(neutronDetector);
|
|
90
|
|
|
91
|
// ----- Create Magnet geometry -----------------------------------------
|
|
92
|
FairModule* magnet = new ERTarget("magnet", 1, kTRUE);
|
|
93
|
magnet->SetGeometryFileName(magnetGeoFileName);
|
|
94
|
run->AddModule(magnet);
|
|
95
|
|
|
96
|
// ----- QTelescope Setup ------------------------------------------------
|
|
97
|
ERQTelescopeSetup* setupQTelescope = ERQTelescopeSetup::Instance();
|
|
98
|
setupQTelescope->SetXmlParametersFile(paramFileQTelescope);
|
|
99
|
|
|
100
|
// ----- T1 parameters ----------------------------------------------------
|
|
101
|
// ----- T1.1--------------------------------------------------------------
|
|
102
|
setupQTelescope->AddSi("DoubleSi_SD1", TVector3( T1Side/2 + T1Aperture/2,
|
|
103
|
T1Side/2 - T1Aperture/2,
|
|
104
|
T1PosZ + T1D1Thick/2), "X");
|
|
105
|
setupQTelescope->AddSi("DoubleSi_SD2", TVector3( T1Side/2 + T1Aperture/2,
|
|
106
|
T1Side/2 - T1Aperture/2,
|
|
107
|
T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
|
|
108
|
// ----- T1.2--------------------------------------------------------------
|
|
109
|
setupQTelescope->AddSi("DoubleSi_SD1", TVector3( T1Side/2 - T1Aperture/2,
|
|
110
|
-T1Side/2 - T1Aperture/2,
|
|
111
|
T1PosZ + T1D1Thick/2), "X");
|
|
112
|
setupQTelescope->AddSi("DoubleSi_SD2", TVector3( T1Side/2 - T1Aperture/2,
|
|
113
|
-T1Side/2 - T1Aperture/2,
|
|
114
|
T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
|
|
115
|
// ----- T1.3 -------------------------------------------------------------
|
|
116
|
setupQTelescope->AddSi("DoubleSi_SD1", TVector3(-T1Side/2 - T1Aperture/2,
|
|
117
|
-T1Side/2 + T1Aperture/2,
|
|
118
|
T1PosZ + T1D1Thick/2), "X");
|
|
119
|
setupQTelescope->AddSi("DoubleSi_SD2", TVector3(-T1Side/2 - T1Aperture/2,
|
|
120
|
-T1Side/2 + T1Aperture/2,
|
|
121
|
T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
|
|
122
|
// ----- T1.4--------------------------------------------------------------
|
|
123
|
setupQTelescope->AddSi("DoubleSi_SD1", TVector3(-T1Side/2 + T1Aperture/2,
|
|
124
|
T1Side/2 + T1Aperture/2,
|
|
125
|
T1PosZ + T1D1Thick/2), "X");
|
|
126
|
setupQTelescope->AddSi("DoubleSi_SD2", TVector3(-T1Side/2 + T1Aperture/2,
|
|
127
|
T1Side/2 + T1Aperture/2,
|
|
128
|
T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
|
|
129
|
|
|
130
|
// ----- D1 parameters ----------------------------------------------------
|
|
131
|
setupQTelescope->AddSi("DoubleSi_D1", TVector3( 0,
|
|
132
|
0,
|
|
133
|
D1PosZ + D1Thick/2), "X");
|
|
134
|
|
|
135
|
// ------QTelescope -------------------------------------------------------
|
|
136
|
ERQTelescope* qtelescope= new ERQTelescope("ERQTelescope", kTRUE,verbose);
|
|
137
|
run->AddModule(qtelescope);
|
|
138
|
|
|
139
|
// ------BeamDet ----------------------------------------------------------
|
|
140
|
ERBeamDet* beamdet= new ERBeamDet("ERBeamDet", kTRUE,verbose);
|
|
141
|
run->AddModule(beamdet);
|
|
142
|
|
|
143
|
FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
|
|
144
|
|
|
145
|
Double_t kinE_MevPerNucleon = 40.;
|
|
146
|
// Int_t Z = 1, A = 3, Q = 1;
|
|
147
|
// TString ionName = "3H";
|
|
148
|
Int_t Z = 2, A = 6, Q = 2;
|
|
149
|
TString ionName = "6He";
|
|
150
|
ERIonMixGenerator* generator = new ERIonMixGenerator(ionName, Z, A, Q, 1);
|
|
151
|
Double32_t kin_energy = kinE_MevPerNucleon * 1e-3 * A; //GeV
|
|
152
|
generator->SetKinE(kin_energy);
|
|
153
|
generator->SetPSigmaOverP(0);
|
|
154
|
Double32_t sigmaTheta = 0.004*TMath::RadToDeg();
|
|
155
|
generator->SetThetaSigma(0, sigmaTheta);
|
|
156
|
// generator->SetThetaRange(0., 5.);
|
|
157
|
generator->SetPhiRange(0, 360);
|
|
158
|
generator->SetBoxXYZ(0, 0, 0, 0, beamStartPosition);
|
|
159
|
|
|
160
|
primGen->AddGenerator(generator);
|
|
161
|
run->SetGenerator(primGen);
|
|
162
|
// ------- Decayer --------------------------------------------------------
|
|
163
|
ERDecayer* decayer = new ERDecayer();
|
|
164
|
ERDecayEXP1803* targetDecay = new ERDecayEXP1803();
|
|
165
|
targetDecay->SetTargetVolumeName("tubeH2");
|
|
166
|
targetDecay->SetTargetThickness(targetH2Thickness);
|
|
167
|
decayer->AddDecay(targetDecay);
|
|
168
|
run->SetDecayer(decayer);
|
|
169
|
|
|
170
|
// ------- Magnetic field -------------------------------------------------
|
|
171
|
ERFieldMap* magField = new ERFieldMap("testField","A");
|
|
172
|
run->SetField(magField,"magnet");
|
|
173
|
// ------- QTelescope Digitizer -------------------------------------------
|
|
174
|
ERQTelescopeDigitizer* qtelescopeDigitizer = new ERQTelescopeDigitizer(verbose);
|
|
175
|
qtelescopeDigitizer->SetSiElossThreshold(0);
|
|
176
|
qtelescopeDigitizer->SetSiElossSigma(0);
|
|
177
|
qtelescopeDigitizer->SetSiTimeSigma(0);
|
|
178
|
|
|
179
|
qtelescopeDigitizer->SetCsIElossThreshold(0);
|
|
180
|
qtelescopeDigitizer->SetCsIElossSigma(0);
|
|
181
|
qtelescopeDigitizer->SetCsITimeSigma(0);
|
|
182
|
run->AddTask(qtelescopeDigitizer);
|
|
183
|
// ------ Gadast Digitizer -----------------------------------------------
|
|
184
|
ERGadastDigitizer* gadastDigitizer = new ERGadastDigitizer(verbose);
|
|
185
|
run->AddTask(gadastDigitizer);
|
|
186
|
// ----- BeamDet Digitizer ----------------------------------------------
|
|
187
|
ERBeamDetDigitizer* beamDetDigitizer = new ERBeamDetDigitizer(verbose);
|
|
188
|
// beamDetDigitizer->SetMWPCElossThreshold(0.006);
|
|
189
|
// beamDetDigitizer->SetToFElossThreshold(0.006);
|
|
190
|
// beamDetDigitizer->SetToFElossSigmaOverEloss(0);
|
|
191
|
// beamDetDigitizer->SetToFTimeSigma(1e-10);
|
|
192
|
run->AddTask(beamDetDigitizer);
|
|
193
|
// -----------------------BeamDetTrackFinder------------------------------
|
|
194
|
ERBeamDetTrackFinder* trackFinder = new ERBeamDetTrackFinder(1);
|
|
195
|
run->AddTask(trackFinder);
|
|
196
|
// // ------------------------------------------------------------------------
|
|
197
|
// -----------------------BeamDetTrackPID-------------------------------
|
|
198
|
ERBeamDetPID* pid = new ERBeamDetPID(1);
|
|
199
|
pid->SetIon(ionName, Z, A, Q);
|
|
200
|
//pid->SetBoxPID(203., 206., 0.005, 0.12);
|
|
201
|
pid->SetBoxPID(0., 1000., 0., 1000.);
|
|
202
|
pid->SetOffsetToF(0.);
|
|
203
|
pid->SetProbabilityThreshold(0);
|
|
204
|
|
|
205
|
run->AddTask(pid);
|
|
206
|
//-------Set visualisation flag to true------------------------------------
|
|
207
|
run->SetStoreTraj(kTRUE);
|
|
208
|
//-------Set LOG verbosity -----------------------------------------------
|
|
209
|
FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
|
|
210
|
// ----- Initialize simulation run ------------------------------------
|
|
211
|
run->Init();
|
|
212
|
Int_t nSteps = -15000;
|
|
213
|
// ----- Runtime database ---------------------------------------------
|
|
214
|
Bool_t kParameterMerged = kTRUE;
|
|
215
|
FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
|
|
216
|
parOut->open(parFile.Data());
|
|
217
|
rtdb->setOutput(parOut);
|
|
218
|
rtdb->saveOutput();
|
|
219
|
rtdb->print();
|
|
220
|
//gMC->SetMaxNStep(nSteps);
|
|
221
|
// ----- Run simulation ------------------------------------------------
|
|
222
|
run->Run(nEvents);
|
|
223
|
|
|
224
|
// ----- Finish -------------------------------------------------------
|
|
225
|
timer.Stop();
|
|
226
|
Double_t rtime = timer.RealTime();
|
|
227
|
Double_t ctime = timer.CpuTime();
|
|
228
|
cout << endl << endl;
|
|
229
|
cout << "Macro finished succesfully." << endl;
|
|
230
|
cout << "Output file is sim.root" << endl;
|
|
231
|
cout << "Parameter file is par.root" << endl;
|
|
232
|
cout << "Real time " << rtime << " s, CPU time " << ctime
|
|
233
|
<< "s" << endl << endl;
|
|
234
|
}
|