SND@LHC Software
Loading...
Searching...
No Matches
Scifi.cxx
Go to the documentation of this file.
1#include "Scifi.h"
2#include "ScifiPoint.h"
3
4#include "TGeoManager.h"
5#include "FairRun.h" // for FairRun
6#include "FairRuntimeDb.h" // for FairRuntimeDb
7#include "TList.h" // for TListIter, TList (ptr only)
8#include "TObjArray.h" // for TObjArray
9#include "TString.h" // for TString
10
11#include "TClonesArray.h"
12#include "TVirtualMC.h"
13
14#include "TGeoBBox.h"
15#include "TGeoTrd1.h"
16#include "TGeoCompositeShape.h"
17#include "TGeoTube.h"
18#include "TGeoMaterial.h"
19#include "TGeoMedium.h"
20#include "TParticle.h"
21
22#include "FairVolume.h"
23#include "FairGeoVolume.h"
24#include "FairGeoNode.h"
25#include "FairRootManager.h"
26#include "FairGeoLoader.h"
27#include "FairGeoInterface.h"
28#include "FairGeoMedia.h"
29#include "FairGeoBuilder.h"
30#include "FairRun.h"
31#include "FairRuntimeDb.h"
32
33#include "ShipDetectorList.h"
34#include "ShipUnit.h"
35#include "ShipStack.h"
36
37#include "TGeoUniformMagField.h"
38#include <stddef.h> // for NULL
39#include <iostream> // for operator<<, basic_ostream, etc
40
41using std::cout;
42using std::endl;
43
44using namespace ShipUnit;
45
47: FairDetector("Scifi", "", kTRUE),
48fTrackID(-1),
49fVolumeID(-1),
50fPos(),
51fMom(),
52fTime(-1.),
53fLength(-1.),
54fELoss(-1),
55eventHeader(0),
56last_run_time(-1),
57last_run_pos(-1),
58last_time_alignment_tag(""),
59last_position_alignment_tag(""),
60alignment_init(false),
61fScifiPointCollection(new TClonesArray("ScifiPoint"))
62{
63}
64
65Scifi::Scifi(const char* name, Bool_t Active, const char* Title)
66: FairDetector(name, true, kLHCScifi),
67fTrackID(-1),
68fVolumeID(-1),
69fPos(),
70fMom(),
71fTime(-1.),
72fLength(-1.),
73fELoss(-1),
74eventHeader(0),
75last_run_time(-1),
76last_run_pos(-1),
77last_time_alignment_tag(""),
78last_position_alignment_tag(""),
79alignment_init(false),
80fScifiPointCollection(new TClonesArray("ScifiPoint"))
81{
82}
83
85{
87 fScifiPointCollection->Delete();
89 }
90}
91
93{
94 FairDetector::Initialize();
95}
96
97// ----- Private method InitMedium
98Int_t Scifi::InitMedium(const char* name)
99{
100 static FairGeoLoader *geoLoad=FairGeoLoader::Instance();
101 static FairGeoInterface *geoFace=geoLoad->getGeoInterface();
102 static FairGeoMedia *media=geoFace->getMedia();
103 static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder();
104
105 FairGeoMedium *ShipMedium=media->getMedium(name);
106
107 if (!ShipMedium)
108 {
109 Fatal("InitMedium","Material %s not defined in media file.", name);
110 return -1111;
111 }
112 TGeoMedium* medium=gGeoManager->GetMedium(name);
113 if (medium!=NULL)
114 return ShipMedium->getMediumIndex();
115 return geoBuild->createMedium(ShipMedium);
116}
117
119{
120 InitMedium("CarbonComposite");
121 TGeoMedium *CarbonComposite = gGeoManager->GetMedium("CarbonComposite");
122
123 InitMedium("rohacell");
124 TGeoMedium *rohacell = gGeoManager->GetMedium("rohacell");
125
126 InitMedium("air");
127 TGeoMedium *air = gGeoManager->GetMedium("air");
128
129 InitMedium("Polycarbonate");
130 TGeoMedium *PlasticBase = gGeoManager->GetMedium("Polycarbonate");
131
132 InitMedium("Polystyrene");
133 TGeoMedium *Polystyrene = gGeoManager->GetMedium("Polystyrene");
134
135 InitMedium("PMMA");
136 TGeoMedium *PMMA = gGeoManager->GetMedium("PMMA");
137
138 InitMedium("PMMA2");
139 TGeoMedium *PMMA2 = gGeoManager->GetMedium("PMMA2");
140
141 InitMedium("Epoxy");
142 TGeoMedium *Epoxy = gGeoManager->GetMedium("Epoxy");
143
144 TGeoVolume *volTarget = gGeoManager->GetVolume("volTarget");
145
146 InitMedium("iron");
147 TGeoMedium *Fe =gGeoManager->GetMedium("iron");
148
149// external parameters
150 Double_t fXDimension = conf_floats["Scifi/xdim"];
151 Double_t fYDimension = conf_floats["Scifi/ydim"];
152 Double_t fZDimension = conf_floats["Scifi/zdim"];
153
154 Double_t fWidthScifiMat = conf_floats["Scifi/scifimat_width"];
155 Double_t fLengthScifiMat = conf_floats["Scifi/scifimat_length"];
156 Double_t fZScifiMat = conf_floats["Scifi/scifimat_z"];
157 Double_t fZEpoxyMat = conf_floats["Scifi/epoxymat_z"];
158 Double_t fGapScifiMat = conf_floats["Scifi/scifimat_gap"]; //dead zone between mats
159
160 Double_t fFiberLength = conf_floats["Scifi/fiber_length"];
161 Double_t fScintCore_rmax = conf_floats["Scifi/scintcore_rmax"];
162 Double_t fClad1_rmin = conf_floats["Scifi/clad1_rmin"];
163 Double_t fClad1_rmax = conf_floats["Scifi/clad1_rmax"];
164 Double_t fClad2_rmin = conf_floats["Scifi/clad2_rmin"];
165 Double_t fClad2_rmax = conf_floats["Scifi/clad2_rmax"];
166
167 Double_t fHorPitch = conf_floats["Scifi/horizontal_pitch"]; //Fiber position params
168 Double_t fVertPitch = conf_floats["Scifi/vertical_pitch"];
169 Double_t fOffsetRowS = conf_floats["Scifi/rowlong_offset"];
170 Double_t fOffsetRowL = conf_floats["Scifi/rowshort_offset"];
171
172 Double_t fZCarbonFiber = conf_floats["Scifi/carbonfiber_z"];
173 Double_t fZHoneycomb = conf_floats["Scifi/honeycomb_z"];
174 Double_t fZGlue = conf_floats["Scifi/glue_z"];
175 Double_t fZAirgap = conf_floats["Scifi/airgap_z"];
176
177 Double_t fXPlastBar = conf_floats["Scifi/plastbar_x"]; //Dimension of plastic bar
178 Double_t fYPlastBar = conf_floats["Scifi/plastbar_y"];
179 Double_t fZPlastBar = conf_floats["Scifi/plastbar_z"];
180
181 // Air gaps in baby SciFi modules
182 Double_t fZBabyPlaneGap = conf_floats["Scifi/plane_gap"];
183 Double_t fZBabyTedlarToPlaneGap = conf_floats["Scifi/tedlar_to_plane"];
184 // Offsets between the upstream tedlar sheet and the upstream passive block
185 // It includes the air gap formed btw the baby module frame and the upstream tedlar sheet.
186 Double_t fStationOffset[4] = {0, conf_floats["Scifi/station_offset1"],
187 conf_floats["Scifi/station_offset2"],
188 conf_floats["Scifi/station_offset3"]};
189
190 Int_t fNFibers_Srow = conf_ints["Scifi/nfibers_shortrow"];
191 Int_t fNFibers_Lrow = conf_ints["Scifi/nfibers_longrow"];
192 Int_t fNFibers_z = conf_ints["Scifi/nfibers_z"];
193
194 Double_t fSeparationBrick = conf_floats["Scifi/scifi_separation"]; //Separation between successive SciFi volumes
195 Double_t fZOffset = conf_floats["Scifi/offset_z"];
196 Int_t fNMats = conf_ints["Scifi/nmats"]; //Number of mats in one SciFi plane
197 Int_t fNScifi = conf_ints["Scifi/nscifi"]; //Number of Scifi walls
198 Int_t fNSiPMs = conf_ints["Scifi/nsipm_mat"]; //Number of SiPMs per SciFi mat
199
200 Double_t fWidthChannel = conf_floats["Scifi/channel_width"]; //One channel width
201 Double_t fHeightChannel = conf_floats["Scifi/channel_height"]; //One channel height
202 Double_t fCharr = conf_floats["Scifi/charr_width"]; //Width of an array of 64 channels without gaps
203 Double_t fEdge = conf_floats["Scifi/sipm_edge"]; //Edge at the left and right sides of the SiPM
204 Double_t fCharrGap = conf_floats["Scifi/charr_gap"]; //Gap between two charr
205 Double_t fBigGap = conf_floats["Scifi/sipm_diegap"]; //Gap between two arrays
206 Int_t fNSiPMChan = conf_ints["Scifi/nsipm_channels"]; //Number of channels in each SiPM
207 Double_t firstChannelX = conf_floats["Scifi/firstChannelX"]; //local X Position of first channel in plane
208
209 Int_t PassiveBlockNotCenterred = conf_ints["Scifi/PassiveBlocknotCenterred"]; // flag showing whether the passive material (FeBlocks) are to be centerred wrt the Scifi
210
211//edge positions in TI18 survey system:
212 std::map<int,TVector3> Vedges;
213 Vedges[0]=TVector3(-conf_floats["Scifi/Xpos0"],conf_floats["Scifi/Zpos0"],conf_floats["Scifi/Ypos0"]);
214 Vedges[1]=TVector3(-conf_floats["Scifi/Xpos1"],conf_floats["Scifi/Zpos1"],conf_floats["Scifi/Ypos1"]);
215 Vedges[2]=TVector3(-conf_floats["Scifi/Xpos2"],conf_floats["Scifi/Zpos2"],conf_floats["Scifi/Ypos2"]);
216 Vedges[3]=TVector3(-conf_floats["Scifi/Xpos3"],conf_floats["Scifi/Zpos3"],conf_floats["Scifi/Ypos3"]);
217 Vedges[4]=TVector3(-conf_floats["Scifi/Xpos4"],conf_floats["Scifi/Zpos4"],conf_floats["Scifi/Ypos4"]);
218
219//edge position in Scifi engineering drawing, y down, z towards IP1, pos X left.
220// need y up, z away from IP1, pos X left: y and z need to change sign.
221
222 TVector3 Sedge = TVector3(conf_floats["Scifi/EdgeAX"],-conf_floats["Scifi/EdgeAY"],-conf_floats["Scifi/EdgeAZ"]);
223//first channel position in Scifi engineering drawing:
224 TVector3 SHfirst = TVector3(conf_floats["Scifi/FirstChannelHX"],-conf_floats["Scifi/FirstChannelHY"],-conf_floats["Scifi/FirstChannelHZ"]);
225 TVector3 SVfirst = TVector3(conf_floats["Scifi/FirstChannelVX"],-conf_floats["Scifi/FirstChannelVY"],-conf_floats["Scifi/FirstChannelVZ"]);
226
227//first channel position in sndsw local plane:
228 TVector3 LHfirst = TVector3(conf_floats["Scifi/LfirstChannelHX"],conf_floats["Scifi/LfirstChannelHY"],conf_floats["Scifi/LfirstChannelHZ"]);
229 TVector3 LVfirst = TVector3(conf_floats["Scifi/LfirstChannelVX"],conf_floats["Scifi/LfirstChannelVY"],conf_floats["Scifi/LfirstChannelVZ"]);
230// moving plane to match edges:
231 std::map<int,TVector3> DeltasH;
232 std::map<int,TVector3> DeltasV;
233 for (int i=0;i<5;i++){
234 DeltasH[i] = Vedges[i] - LHfirst + SHfirst - Sedge;
235 DeltasV[i] = Vedges[i] - LVfirst + SVfirst - Sedge;
236 }
237 Double_t totalThickness{};
238
239 //Carbon Fiber Film
240 TGeoVolume *CarbonFiberVolume = gGeoManager->MakeBox("CarbonFiber", CarbonComposite, fXDimension/2, fYDimension/2, fZCarbonFiber/2);
241 CarbonFiberVolume->SetLineColor(kGray - 2);
242 CarbonFiberVolume->SetVisibility(1);
243
244 //Honeycomb Rohacell
245 TGeoVolume *HoneycombVolume = gGeoManager->MakeBox("Honeycomb", rohacell, fXDimension/2, fYDimension/2, fZHoneycomb/2);
246 HoneycombVolume->SetLineColor(kYellow);
247 HoneycombVolume->SetVisibility(1);
248
249 // Glue between Carbon Fiber and Honeycomb or Glue between Carbon Fiber and Scifi mat
250 TGeoVolume *GlueVolume = gGeoManager->MakeBox("Glue", Epoxy, fXDimension/2, fYDimension/2, fZGlue/2);
251 GlueVolume->SetLineColor(kYellow-1);
252 GlueVolume->SetVisibility(1);
253
254 // Air gap in the middle
255 TGeoVolume *AirGapVolume = gGeoManager->MakeBox("Airgap", air, fXDimension/2, fYDimension/2, fZAirgap/2);
256 AirGapVolume->SetLineColor(kGray-1);
257 AirGapVolume->SetVisibility(1);
258
259 //Plastic/Air
260 //Definition of the box containing polycarbonate pieces and an air gap
261 TGeoVolume *PlasticAirVolume = gGeoManager->MakeBox("PlasticAir", air, fXDimension/2, fYDimension/2, fZPlastBar/2);
262 PlasticAirVolume->SetLineColor(kGray-1);
263 PlasticAirVolume->SetVisibility(1);
264 PlasticAirVolume->SetVisDaughters(1);
265
266 //Plastic bars
267 TGeoVolume *PlasticBarVolume = gGeoManager->MakeBox("PlasticBar", PlasticBase, fXPlastBar/2, fYPlastBar/2, fZPlastBar/2);
268 PlasticBarVolume->SetLineColor(kGray-4);
269 PlasticBarVolume->SetVisibility(1);
270
271 PlasticAirVolume->AddNode(PlasticBarVolume, 0, new TGeoTranslation(- fXDimension/2 + fXPlastBar/2, 0, 0)); //bars are placed || to y
272 PlasticAirVolume->AddNode(PlasticBarVolume, 1, new TGeoTranslation(+ fXDimension/2 - fXPlastBar/2, 0, 0));
273
274 //Plastic Glue Air
275 //Definition of the box containing glue between plastic bar and scifi mat, and an air in the middel
276 TGeoVolume *PlasticGlueAirVolume = gGeoManager->MakeBox("PlasticGlueAir", air, fXDimension/2, fYDimension/2, fZGlue/2);
277 PlasticGlueAirVolume->SetLineColor(kGray-1);
278 PlasticGlueAirVolume->SetVisibility(1);
279 PlasticGlueAirVolume->SetVisDaughters(1);
280
281 //Plastic glue bars
282 TGeoVolume *PlasticGlueBarVolume = gGeoManager->MakeBox("PlasticGlueBar", Epoxy, fXPlastBar/2, fYPlastBar/2, fZGlue/2);
283 PlasticGlueBarVolume->SetLineColor(kYellow-1);
284 PlasticGlueBarVolume->SetVisibility(1);
285
286 PlasticGlueAirVolume->AddNode(PlasticGlueBarVolume, 0, new TGeoTranslation(- fXDimension/2 + fXPlastBar/2, 0, 0)); //bars are placed || to y
287 PlasticGlueAirVolume->AddNode(PlasticGlueBarVolume, 1, new TGeoTranslation(+ fXDimension/2 - fXPlastBar/2, 0, 0));
288
289 //Definition of the two air gaps for baby SciFi
290 TGeoVolume *PlaneGapVolume = gGeoManager->MakeBox("PlaneGap", air, fXDimension/2, fYDimension/2, fZBabyPlaneGap/2);
291 PlaneGapVolume->SetLineColor(kOrange-4);
292 PlaneGapVolume->SetTransparency(50);
293 PlaneGapVolume->SetVisibility(1);
294 TGeoVolume *TedlarToPlaneGapVolume = gGeoManager->MakeBox("TedlarToPlaneGap", air, fXDimension/2, fYDimension/2, fZBabyTedlarToPlaneGap/2);
295 TedlarToPlaneGapVolume->SetLineColor(kOrange-4);
296 TedlarToPlaneGapVolume->SetTransparency(50);
297 TedlarToPlaneGapVolume->SetVisibility(1);
298
299 //Fiber volume that contains the scintillating core and double cladding
300 TGeoVolumeAssembly *FiberVolume = new TGeoVolumeAssembly("FiberVolume");
301
302 TGeoVolume *ScintCoreVol = gGeoManager->MakeTube("ScintCoreVol", Polystyrene, 0, fScintCore_rmax, fFiberLength/2);
303 TGeoVolume *Clad1Vol = gGeoManager->MakeTube("Clad1Vol", PMMA, fClad1_rmin, fClad1_rmax, fFiberLength/2);
304 TGeoVolume *Clad2Vol = gGeoManager->MakeTube("Clad2Vol", PMMA2, fClad2_rmin, fClad2_rmax, fFiberLength/2);
305
306 FiberVolume->AddNode(ScintCoreVol, 0);
307 FiberVolume->AddNode(Clad1Vol, 0);
308 FiberVolume->AddNode(Clad2Vol, 0);
309 FiberVolume->SetVisDaughters(kFALSE);
310
311 //Add SciFi fiber as sensitive unit
312 AddSensitiveVolume(ScintCoreVol);
313
314 //Fiber and plane rotations
315 TGeoRotation *rothorfiber = new TGeoRotation("rothorfiber", 90, 90, 0);
316 TGeoRotation *rotvertfiber = new TGeoRotation("rotvertfiber", 0, 90, 0);
317 TGeoRotation *rot = new TGeoRotation("rot", 90, 180, 0);
318
319 //SciFi mats for X and Y fiber planes
320 Double_t MatThickness{};
321 TGeoMedium *MatMaterial;
322 if (fNScifi==5){ //TI18 SciFi
323 MatMaterial = Epoxy;
324 MatThickness = fZEpoxyMat;
325 }
326 if (fNScifi==4){ // no glue in baby SciFi planes
327 MatMaterial = air;
328 MatThickness = fZScifiMat;
329 }
330 TGeoVolume *HorMatVolume = gGeoManager->MakeBox("HorMatVolume", MatMaterial, fLengthScifiMat/2, fWidthScifiMat/2, MatThickness/2);
331 TGeoVolume *VertMatVolume = gGeoManager->MakeBox("VertMatVolume", MatMaterial, fWidthScifiMat/2, fLengthScifiMat/2, MatThickness/2);
332
333 Double_t zPosM;
334 Double_t offsetS = -fWidthScifiMat/2 + fOffsetRowS;
335 Double_t offsetL = -fWidthScifiMat/2 + fOffsetRowL;
336
337 // All fibres will be assigned station number 1 and mat number 1, to keep compatibility with the STMRFFF format.
338 int dummy_station = 1;
339 int dummy_mat = 1;
340
341 //Adding horizontal fibers
342 for (int irow = 0; irow < fNFibers_z; irow++){
343 zPosM = -fZScifiMat/2 + fClad2_rmax + irow*fVertPitch;
344 if (irow%2 == 0){
345 for (int ifiber = 0; ifiber < fNFibers_Srow; ifiber++){
346 HorMatVolume->AddNode(FiberVolume, 1e6*dummy_station + 1e5*0 + 1e4*dummy_mat + 1e3*(irow + 1) + ifiber + 1, new TGeoCombiTrans("rottranshor0", 0, offsetS + ifiber*fHorPitch, zPosM, rothorfiber));
347 }
348 }
349 else{
350 for (int ifiber = 0; ifiber < fNFibers_Lrow; ifiber++){
351 HorMatVolume->AddNode(FiberVolume, 1e6*dummy_station + 1e5*0 + 1e4*dummy_mat + 1e3*(irow + 1) + ifiber + 1, new TGeoCombiTrans("rottranshor1", 0, offsetL + ifiber*fHorPitch, zPosM, rothorfiber));
352 }
353 }
354 }
355
356 //Adding vertical fibers
357 for (int irow = 0; irow < fNFibers_z; irow++){
358 zPosM = -fZScifiMat/2 + fClad2_rmax + irow*fVertPitch;
359 if (irow%2 == 0){
360 for (int ifiber = 0; ifiber < fNFibers_Srow; ifiber++){
361 VertMatVolume->AddNode(FiberVolume, 1e6*dummy_station + 1e5*1 + 1e4*dummy_mat + 1e3*(irow + 1) + ifiber + 1, new TGeoCombiTrans("rottransvert0", offsetS + ifiber*fHorPitch, 0, zPosM, rotvertfiber));
362 }
363 }
364 else{
365 for (int ifiber = 0; ifiber < fNFibers_Lrow; ifiber++){
366 VertMatVolume->AddNode(FiberVolume, 1e6*dummy_station + 1e5*1 + 1e4*dummy_mat + 1e3*(irow + 1) + ifiber + 1, new TGeoCombiTrans("rottransvert1", offsetL + ifiber*fHorPitch, 0, zPosM, rotvertfiber));
367 }
368 }
369 }
370
371 // for testbeam 2023
372 // for now the distinct feature of the testbeam could be the 4 SciFi planes
373 std::map<int, TGeoVolume*> volFeTarget;
374 std::map<int, float> fFeTargetX;
375 std::map<int, float> fFeTargetY;
376 std::map<int, float> fFeTargetZ;
377 for ( int i = 0; i < fNScifi; i++){
378 std::string station = std::to_string(i+1);
379 fFeTargetX[i] = conf_floats[TString("Scifi/FeTargetX"+station)];
380 fFeTargetY[i] = conf_floats[TString("Scifi/FeTargetZ"+station)];
381 fFeTargetZ[i] = conf_floats[TString("Scifi/FeTargetY"+station)];
382 }
383
384 // DetID is of the form:
385 // first digit - station number
386 // second digit - type of the plane: 0-horizontal fiber, 1-vertical fiber
387 // third digit - mat number
388 // fourth digit - row number (in Z direction)
389 // last three digits - fiber number
390 // e.g. DetID = 1021074 -> station 1, horizontal fiber plane, mat 2, row 1, fiber 74
391 for (int istation = 0; istation < fNScifi; istation++){
392 Int_t node = 1e6*(istation+1);
393 std::string station = std::to_string(istation+1);
394 TGeoVolumeAssembly *ScifiVolume = new TGeoVolumeAssembly( TString("ScifiVolume"+station) );
395 TGeoVolumeAssembly *ScifiHorPlaneVol = new TGeoVolumeAssembly( TString("ScifiHorPlaneVol"+station) );
396 TGeoVolumeAssembly *ScifiVertPlaneVol = new TGeoVolumeAssembly( TString("ScifiVertPlaneVol"+station) );
397
398 //TI18 SciFi
399 if (fNScifi==5){
400 //Adding the first half of the SciFi module that contains horizontal fibres
401 ScifiVolume->AddNode(CarbonFiberVolume, 0, new TGeoTranslation(0, 0, fZCarbonFiber/2));
402 ScifiVolume->AddNode(GlueVolume, 0, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue/2));
403 ScifiVolume->AddNode(HoneycombVolume, 0, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb/2));
404 ScifiVolume->AddNode(GlueVolume, 1, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue/2));
405 ScifiVolume->AddNode(CarbonFiberVolume, 1, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber/2));
406 ScifiVolume->AddNode(GlueVolume, 2, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber + fZGlue/2));
407 ScifiVolume->AddNode(ScifiHorPlaneVol, node, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber + fZGlue + fZEpoxyMat/2));
408 ScifiVolume->AddNode(PlasticGlueAirVolume, 0, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber + fZGlue + fZEpoxyMat + fZGlue/2));
409 ScifiVolume->AddNode(PlasticAirVolume, 0, new TGeoTranslation(0, 0, fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber + fZGlue + fZEpoxyMat + fZGlue + fZPlastBar/2));
410
411 Double_t first_half_z = fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber + fZGlue + fZEpoxyMat + fZGlue + fZPlastBar;
412 ScifiVolume->AddNode(AirGapVolume, 0, new TGeoTranslation(0, 0, first_half_z + fZAirgap/2));
413
414 //Adding the second half of the SciFi module that contains vertical fibres
415 ScifiVolume->AddNode(PlasticAirVolume, 1, new TGeoCombiTrans("rottrans0", 0, 0, first_half_z + fZAirgap + fZPlastBar/2, rot));
416 ScifiVolume->AddNode(PlasticGlueAirVolume, 1, new TGeoCombiTrans("rottrans0",0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue/2, rot));
417 ScifiVolume->AddNode(ScifiVertPlaneVol, node, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat/2));
418 ScifiVolume->AddNode(GlueVolume, 3, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue/2));
419 ScifiVolume->AddNode(CarbonFiberVolume, 2, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber/2));
420 ScifiVolume->AddNode(GlueVolume, 4, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber + fZGlue/2));
421 ScifiVolume->AddNode(HoneycombVolume, 0, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber + fZGlue + fZHoneycomb/2));
422 ScifiVolume->AddNode(GlueVolume, 5, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue/2));
423 ScifiVolume->AddNode(CarbonFiberVolume, 3, new TGeoTranslation(0, 0, first_half_z + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber/2));
424
425 // Double_t totalThickness = fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber + fZGlue + fZEpoxyMat + fZGlue + fZPlastBar + fZAirgap + fZPlastBar + fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber;
426 totalThickness = fZCarbonFiber + fZGlue + fZHoneycomb + fZGlue + fZCarbonFiber +
427 fZGlue + fZEpoxyMat + fZGlue + fZPlastBar + fZAirgap + fZPlastBar +
428 fZGlue + fZEpoxyMat + fZGlue + fZCarbonFiber + fZGlue + fZHoneycomb +
429 fZGlue + fZCarbonFiber;
430 }
431
432 //Baby SciFi
433 // No glue in this case and one uses fZScifiMat and not the fZEpoxyMat for the mat thickness
434 if (fNScifi==4){
435 //Adding the first half of the SciFi module that contains horizontal fibres
436 ScifiVolume->AddNode(TedlarToPlaneGapVolume, 1, new TGeoTranslation(0, 0, fZBabyTedlarToPlaneGap/2));
437 ScifiVolume->AddNode(ScifiHorPlaneVol, node, new TGeoTranslation(0, 0, fZBabyTedlarToPlaneGap + fZScifiMat/2));
438
439 // Add the air gap between the X and Y plane
440 ScifiVolume->AddNode(PlaneGapVolume, 0, new TGeoTranslation(0, 0, fZBabyTedlarToPlaneGap + fZScifiMat + fZBabyPlaneGap/2));
441
442 //Adding the second half of the SciFi module that contains vertical fibres
443 ScifiVolume->AddNode(ScifiVertPlaneVol, node, new TGeoTranslation(0, 0, fZBabyTedlarToPlaneGap + 3*fZScifiMat/2 + fZBabyPlaneGap));
444 ScifiVolume->AddNode(TedlarToPlaneGapVolume, 2, new TGeoTranslation(0, 0, 3*fZBabyTedlarToPlaneGap/2 + 2*fZScifiMat + fZBabyPlaneGap));
445
446 totalThickness = 2*fZBabyTedlarToPlaneGap + 2*fZScifiMat + fZBabyPlaneGap;
447 }
448
449 volTarget->AddNode(ScifiVolume, node,
450 new TGeoTranslation(DeltasV[istation][0], DeltasH[istation][1], DeltasH[istation][2]));
451
452 // std::cout<<"deltas "<<DeltasV[istation][0]<<" "<< DeltasH[istation][1]<<" "<< DeltasH[istation][2]<<" "<<totalThickness<<std::endl;
453 // for 2023 testbeam put iron blocks of variable length in between SciFi planes - the planes are dowstream of the blocks!
454 // for the 2024 testbeam, the iron blocks are aligned to Brick 1 and laid on the concrete (same as SciFi 1)
455 if (fNScifi==4 && istation != 0) {
456 volFeTarget[istation] = gGeoManager->MakeBox(TString("volFeTarget"+station),Fe,fFeTargetX[istation]/2., fFeTargetY[istation]/2., fFeTargetZ[istation]/2.);
457 volFeTarget[istation]->SetLineColor(kGreen-4);
458 volTarget->AddNode(volFeTarget[istation],1,
459 new TGeoTranslation(DeltasV[istation][0] - PassiveBlockNotCenterred*fabs(fXDimension-fFeTargetX[istation])/2.,
460 DeltasH[istation][1] + PassiveBlockNotCenterred*(DeltasH[0][1]-DeltasH[istation][1]
461 +fabs(fYDimension-fFeTargetY[istation])/2.),
462 DeltasH[istation][2] - fStationOffset[istation] - fFeTargetZ[istation]/2.));
463 }
464
465 //Creating Scifi planes by appending fiber mats
466 for (int imat = 0; imat < fNMats; imat++){
467 int N = fNMats==1 ? imat : imat-1;
468
469 //Placing mats along Y
470 ScifiHorPlaneVol->AddNode(HorMatVolume, 1e6*(istation+1) + 1e4*(imat + 1), new TGeoTranslation(0, N*(fWidthScifiMat+fGapScifiMat), 0));
471
472 //Placing mats along X
473 ScifiVertPlaneVol->AddNode(VertMatVolume, 1e6*(istation+1) + 1e5 + 1e4*(imat + 1), new TGeoTranslation(N*(fWidthScifiMat+fGapScifiMat), 0, 0));
474 }
475 }
476
477}
478
480{
481 if (gGeoManager->FindVolumeFast("SiPMmapVol")){return;}
482 Double_t fLengthScifiMat = conf_floats["Scifi/scifimat_length"];
483 Double_t fWidthChannel = conf_floats["Scifi/channel_width"];
484 Double_t fHeightChannel = conf_floats["Scifi/channel_height"];
485 Int_t fNSiPMChan = conf_ints["Scifi/nsipm_channels"];
486 Int_t fNSiPMs = conf_ints["Scifi/nsipm_mat"];
487 Int_t fNMats = conf_ints["Scifi/nmats"];
488 Double_t fEdge = conf_floats["Scifi/sipm_edge"];
489 Double_t fCharr = conf_floats["Scifi/charr_width"];
490 Double_t fCharrGap = conf_floats["Scifi/charr_gap"];
491 Double_t fBigGap = conf_floats["Scifi/sipm_diegap"];
492 Double_t firstChannelX = conf_floats["Scifi/firstChannelX"];
493
494 //Contains all plane SiPMs, defined for horizontal fiber plane
495 //To obtain SiPM map for vertical fiber plane rotate by 90 degrees around Z
496 TGeoVolumeAssembly *SiPMmapVol = new TGeoVolumeAssembly("SiPMmapVol");
497
498 TGeoVolume*ChannelVol = gGeoManager->MakeBox("ChannelVol", 0, fLengthScifiMat/2, fWidthChannel/2, fHeightChannel/2);
499
500 //DetID for each channel:
501 //first digit: mat number (0-2)
502 //second digit: SiPM number (0-3)
503 //last three digits: channel number (0-127)
504
505 Double_t SiPMArray_fullwidth = fEdge+fCharr+fCharrGap+fCharr+fEdge;
506 TGeoVolumeAssembly *SiPMArrayVol;
507 int N = fNMats == 1 ? 1 : 0;
508 Double_t pos = -fEdge+firstChannelX + N*fLengthScifiMat;
509 for (int imat = 0; imat < fNMats; imat++){
510 for (int isipms = 0; isipms < fNSiPMs; isipms++){
511 pos+= fEdge;
512 for (int ichannel = 0; ichannel < fNSiPMChan; ichannel++){
513 SiPMmapVol->AddNode(ChannelVol, imat*10000+isipms *1000 + ichannel, new TGeoTranslation(0, pos, 0));
514 pos += fWidthChannel;
515 if (ichannel==(fNSiPMChan/2-1)){pos += fCharrGap;}
516 }
517 pos+=fEdge+fBigGap;
518 }
519 }
520}
521
523 // get mapping to eventHeader
524 eventHeader = e;
525
526 // Initialize
527 if (not alignment_init) {
528 alignment_init = true;
529 // Get available tags from the geometry file
530 std::string tag_string;
531
532 // Time alignment tags
533 for (auto key : conf_floats){
534 tag_string = key.first.Data();
535 if (tag_string.find("Scifi/station1t_") != std::string::npos){
536 covered_runs_time_alignment.push_back(stoi(tag_string.substr(tag_string.find("t_")+2)));
537 }
538 }
540 // Position alignment tags
541 for (auto key : conf_floats){
542 tag_string = key.first.Data();
543 if (tag_string.find("Scifi/LocM100t_") != std::string::npos){
544 covered_runs_position_alignment.push_back(stoi(tag_string.substr(tag_string.find("t_")+2)));
545 }
546 }
548 }
549};
550
551Bool_t Scifi::ProcessHits(FairVolume* vol)
552{
554 //Set parameters at entrance of volume. Reset ELoss.
555 if ( gMC->IsTrackEntering() )
556 {
557 fELoss = 0.;
558 fTime = gMC->TrackTime() * 1.0e09;
559 fLength = gMC->TrackLength();
560 gMC->TrackPosition(fPos);
561 gMC->TrackMomentum(fMom);
562 TGeoNavigator* nav = gGeoManager->GetCurrentNavigator();
563
564 int fibre_local_id = nav->GetMother()->GetNumber() - 1e6 - 1e4; // Local ID within the mat.
565 int fibre_mat_id = nav->GetMother(2)->GetNumber(); // Get mat ID.
566
567 int fibre_station_number = int( fibre_mat_id / 1e6); // Get the station number from the mat
568
569 int fibre_mat_id_station_removed = fibre_mat_id - fibre_station_number*1e6;
570 int fibre_mat_number = int((fibre_mat_id_station_removed - int(fibre_mat_id_station_removed/1e5)*1e5)/1e4);
571
572 fVolumeID = fibre_local_id + fibre_station_number*1e6 + fibre_mat_number*1e4;
573
574 if (fVolumeID==0){std::cout<<"fiber vol id "<<nav->GetMother()->GetName()<<std::endl;}
575
576 }
577 // Sum energy loss for all steps in the active volume
578 fELoss += gMC->Edep();
579
580 // Create ScifiPoint at exit of active volume
581 if ( gMC->IsTrackExiting() ||
582 gMC->IsTrackStop() ||
583 gMC->IsTrackDisappeared() ) {
584 if (fELoss == 0. ) { return kFALSE; }
585 fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
586/* STMRFFF
587First digit S: station # within the sub-detector
588Second digit T: type of the plane: 0-horizontal fiber plane, 1-vertical fiber plane
589Third digit M: determines the mat number
590Fourth digit R: row number (in Z direction)
591Last three digits F: fiber number
592*/
593 TParticle* p=gMC->GetStack()->GetCurrentTrack();
594 Int_t pdgCode = p->GetPdgCode();
595
596 TLorentzVector Pos;
597 gMC->TrackPosition(Pos);
598 Double_t xmean = (fPos.X()+Pos.X())/2. ;
599 Double_t ymean = (fPos.Y()+Pos.Y())/2. ;
600 Double_t zmean = (fPos.Z()+Pos.Z())/2. ;
601
602 AddHit(fTrackID,fVolumeID, TVector3(xmean, ymean, zmean),
603 TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, fLength,
604 fELoss, pdgCode);
605
606 // Increment number of det points in TParticle
607 ShipStack* stack = (ShipStack*) gMC->GetStack();
608 stack->AddPoint(kLHCScifi);
609 }
610
611 return kTRUE;
612}
613
614Double_t Scifi::GetCorrectedTime(Int_t fDetectorID, Double_t rawTime, Double_t L){
615/* expect time in u.ns and path length to sipm u.cm */
616 TString tag = "t";
617 TString sID;
618
619 if (eventHeader){
620 Int_t fRunNumber = eventHeader->GetRunId();
621 if (fRunNumber != last_run_time){
622 last_run_time = fRunNumber;
623 if (fRunNumber<1) {
624 LOG(ERROR) << "Scifi::GetCorrectedTime: non valid run number "<<fRunNumber;
625 return rawTime;
626 }
627
628 if (covered_runs_time_alignment.size()!=0){
629 tag = "t_"+std::to_string(covered_runs_time_alignment[covered_runs_time_alignment.size()-1]);
630 for (int i=1; i<covered_runs_time_alignment.size(); i++){
631 if (fRunNumber>=covered_runs_time_alignment[i-1] && fRunNumber<covered_runs_time_alignment[i]){
632 tag = "t_"+std::to_string(covered_runs_time_alignment[i-1]);
633 }
634 }
635 // special case
636 if (fRunNumber<5193 && fRunNumber>5174) tag = "t_"+std::to_string(covered_runs_time_alignment[0]);
637 }
638 else{
639 // allow reading older geo files with letter tags i.e. A, B, C
640 tag = "tA";
641 if (fRunNumber>5116 && !(fRunNumber<5193 && fRunNumber>5174) ) {tag = "tB";}
642 }
643 // 2023 and 2024 testbeam data don't have custom tags
644 if (fRunNumber>=1e5) {tag = "t";}
646 }
647 }
648 sID.Form("%i",fDetectorID);
649 Double_t cor = conf_floats["Scifi/station"+TString(sID(0,1))+last_time_alignment_tag];
650 if (sID(1,1)=="0"){
651 // In the teatbeam 2024, SciFi 2H needs internal time corrections per channel
652 if (conf_ints["Scifi/channelTimeAlignment"]==1 && floor(fDetectorID/100000)==20) {
653 cor+=conf_vectors["Scifi/station"+TString(sID(0,4))+"XXX"+last_time_alignment_tag][fDetectorID%1000];
654 }
655 cor+=conf_floats["Scifi/station"+TString(sID(0,1))+"H"+TString(sID(2,1))+last_time_alignment_tag];
656 }
657 else{
658 cor+=conf_floats["Scifi/station"+TString(sID(0,1))+"V"+TString(sID(2,1))+last_time_alignment_tag];
659 }
660 cor += L/conf_floats["Scifi/signalSpeed"];
661 return rawTime-cor;
662}
663
664void Scifi::GetPosition(Int_t fDetectorID, TVector3& A, TVector3& B)
665{
666// TGeoVolumeAssembly *SiPMmapVol = gGeoManager->FindVolumeFast("SiPMmapVol");
667// if(!SiPMmapVol ){SiPMmapVol=SiPMOverlap();}
668
669/* STMRFFF
670 First digit S: station # within the sub-detector
671 Second digit T: type of the plane: 0-horizontal fiber plane, 1-vertical fiber plane
672 Third digit M: determines the mat number
673 Fourth digit R: row number (in Z direction)
674 Last three digits F: fiber number
675*/
676
677 Int_t station_number = int(fDetectorID/1e6);
678 Int_t mat_number = int(fDetectorID/1e4)%int(fDetectorID/1e5);
679
680 Int_t local_fibre_id = fDetectorID - (station_number-1)*1e6 - (mat_number-1)*1e4;
681 TString sLocalID;
682 sLocalID.Form("%i", local_fibre_id);
683
684 TString sID;
685 sID.Form("%i",fDetectorID);
686 TString path = "/cave_1/Detector_0/volTarget_1/ScifiVolume"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000/";
687 if (sID(1,1)=="0"){
688 path+="ScifiHorPlaneVol"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000/";
689 path+="HorMatVolume_"+TString(sID(0,3))+"0000/";
690 }else{
691 path+="ScifiVertPlaneVol"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000/";
692 path+="VertMatVolume_"+TString(sID(0,3))+"0000/";
693 }
694 path+="FiberVolume_"+sLocalID;
695 TGeoNavigator* nav = gGeoManager->GetCurrentNavigator();
696 nav->cd(path);
697 LOG(DEBUG) <<path<<" "<<fDetectorID;
698 TGeoNode* W = nav->GetCurrentNode();
699 TGeoBBox* S = dynamic_cast<TGeoBBox*>(W->GetVolume()->GetShape());
700
701 Double_t top[3] = {0,0,S->GetDZ()};
702 Double_t bot[3] = {0,0,-(S->GetDZ())};
703 Double_t Gtop[3],Gbot[3];
704 nav->LocalToMaster(top, Gtop); nav->LocalToMaster(bot, Gbot);
705 A.SetXYZ(Gtop[0],Gtop[1],Gtop[2]);
706 B.SetXYZ(Gbot[0],Gbot[1],Gbot[2]);
707
708}
709TVector3 Scifi::GetLocalPos(Int_t id, TVector3* glob){
710 TString sID;
711 sID.Form("%i",id);
712 TString path = "/cave_1/Detector_0/volTarget_1/ScifiVolume"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000/";
713 if (sID(1,1)=="0"){
714 path+="ScifiHorPlaneVol"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000";
715 }else{
716 path+="ScifiVertPlaneVol"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000";
717 }
718 TGeoNavigator* nav = gGeoManager->GetCurrentNavigator();
719 nav->cd(path);
720 Double_t aglob[3];
721 Double_t aloc[3];
722 glob->GetXYZ(aglob);
723 nav->MasterToLocal(aglob,aloc);
724 return TVector3(aloc[0],aloc[1],aloc[2]);
725}
726
727void Scifi::GetSiPMPosition(Int_t SiPMChan, TVector3& A, TVector3& B)
728{
729/* STMRFFF
730 First digit S: station # within the sub-detector
731 Second digit T: type of the plane: 0-horizontal fiber plane, 1-vertical fiber plane
732 Third digit M: determines the mat number 0-2
733 Fourth digit S: SiPM number 0-3
734 Last three digits F: local SiPM channel number in one mat 0-127
735*/
736 Int_t locNumber = SiPMChan%100000;
737 Int_t globNumber = int(SiPMChan/100000)*100000;
738 Float_t locPosition = SiPMPos[locNumber]; // local position in plane of reference plane.
739 Double_t fFiberLength = conf_floats["Scifi/fiber_length"];
740 Int_t fNMats = conf_ints["Scifi/nmats"];
741
742 TString tag = "";
743
744 // in case of old data with FairEventHeader, user will be responsible to use the correct geofile.
745 if (eventHeader){
746 Int_t fRunNumber = eventHeader->GetRunId();
747 if (fRunNumber != last_run_pos){
748 last_run_pos = fRunNumber;
749
750 if (fRunNumber<1) {
751 LOG(ERROR) << "Scifi::GetSiPMPosition: non valid run number "<<fRunNumber;
752 return;
753 }
754
755 if (covered_runs_position_alignment.size()!=0){
756 tag = "t_"+std::to_string(covered_runs_position_alignment[covered_runs_position_alignment.size()-1]);
757 for (int i=1; i<covered_runs_position_alignment.size(); i++){
758 if (fRunNumber>=covered_runs_position_alignment[i-1] && fRunNumber<covered_runs_position_alignment[i]){
759 tag = "t_"+std::to_string(covered_runs_position_alignment[i-1]);
760 }
761 }
762 }
763 else{
764 // allow reading older geo files with letter tags i.e. A, B, C
765 tag = "E";
766 if (fRunNumber<4575) {tag = "A";}
767 else if (fRunNumber<4855) {tag = "B";}
768 else if (fRunNumber<5172) {tag = "C";}
769 else if (fRunNumber<5431) {tag = "D";}
770 }
771 // 2023 testbeam data doesn't have a custom tag
772 if (fRunNumber>=1e5) {tag = "";}
774 }
775 }
776 TString sID;
777 sID.Form("%i",SiPMChan);
778 Int_t digits = fNMats==1 ? 2 : 1;
779 locPosition += conf_floats["Scifi/LocM"+TString(sID(0,3))+last_position_alignment_tag];
780 Float_t rotPhi = conf_floats["Scifi/RotPhiS"+TString(sID(0,digits))+last_position_alignment_tag];
781 Float_t rotPsi = conf_floats["Scifi/RotPsiS"+TString(sID(0,digits))+last_position_alignment_tag];
782 Float_t rotTheta = conf_floats["Scifi/RotThetaS"+TString(sID(0,digits))+last_position_alignment_tag];
783
784 Double_t loc[3] = {0,0,0};
785 TString path = "/cave_1/Detector_0/volTarget_1/ScifiVolume"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000/";
786 TGeoNavigator* nav = gGeoManager->GetCurrentNavigator();
787 Double_t glob[3] = {0,0,0};
788
789 if (sID(1,1)=="0"){
790 path+="ScifiHorPlaneVol"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000";
791 loc[0] = -fFiberLength/2 - (rotPhi + rotPsi)*locPosition ;
792 loc[1] = locPosition - fFiberLength/2 * (rotPhi + rotPsi) ;
793 loc[2] = rotTheta*locPosition;
794 nav->cd(path);
795 nav->LocalToMaster(loc, glob);
796 A.SetXYZ( glob[0], glob[1],glob[2] );
797 loc[0] = fFiberLength/2 - (rotPhi + rotPsi)*locPosition ;
798 loc[1] = locPosition + fFiberLength/2 * (rotPhi + rotPsi) ;
799 loc[2] = - rotTheta*locPosition;
800 nav->LocalToMaster(loc, glob);
801 B.SetXYZ( glob[0], glob[1],glob[2] );
802 }else{
803 path+="ScifiVertPlaneVol"+TString(sID(0,1))+"_"+TString(sID(0,1))+"000000";
804 loc[0] = locPosition + fFiberLength/2*(rotPhi + rotPsi);
805 loc[1] = -fFiberLength/2 + locPosition*(rotPhi + rotPsi);
806 loc[2] = -fFiberLength/2*rotTheta;
807 nav->cd(path);
808 nav->LocalToMaster(loc, glob);
809 A.SetXYZ( glob[0], glob[1],glob[2] );
810 loc[0] = locPosition - fFiberLength/2*(rotPhi + rotPsi);
811 loc[1] = fFiberLength/2 + locPosition*(rotPhi + rotPsi);
812 loc[2] = -fFiberLength/2*rotTheta;
813 nav->LocalToMaster(loc, glob);
814 B.SetXYZ( glob[0], glob[1],glob[2] );
815 }
816}
817
818Double_t Scifi::ycross(Double_t a,Double_t R,Double_t x)
819{
820 Double_t y = -1;
821 Double_t A = R*R - (x-a)*(x-a);
822 if ( !(A<0) ){y = TMath::Sqrt(A);}
823 return y;
824}
825Double_t Scifi::integralSqrt(Double_t ynorm)
826{
827 Double_t y = 1./2.*(ynorm*TMath::Sqrt(1-ynorm*ynorm)+TMath::ASin(ynorm));
828 return y;
829}
830Double_t Scifi::fraction(Double_t R,Double_t x,Double_t y)
831{
832 Double_t F= 2*R*R*(integralSqrt(y/R) );
833 F-=(2*x*y);
834 Double_t result = F/(R*R*TMath::Pi());
835 return result;
836}
837Double_t Scifi::area(Double_t a,Double_t R,Double_t xL,Double_t xR)
838{
839 Double_t fracL = -1;
840 Double_t fracR = -1;
841 if (xL<=a-R && xR>=a+R) {return 1;}
842 Double_t leftC = ycross(a,R,xL);
843 Double_t rightC = ycross(a,R,xR);
844 if (leftC<0 && rightC<0) {return -1;}
845 if ( !(rightC<0) ){ fracR = fraction(R,abs(xR-a),rightC);}
846 if ( !(leftC<0) ) { fracL = fraction(R,abs(xL-a),leftC);}
847 Double_t theAnswer = 0;
848 if ( !(leftC<0) ) {
849 if(xL<a){theAnswer += 1-fracL;}
850 else{ theAnswer += fracL;}
851 if ( !(rightC<0) ) {theAnswer -=1;}
852 }
853 if ( !(rightC<0) ){
854 if(xR>a){ theAnswer += 1-fracR;}
855 else{ theAnswer += fracR;}
856 }
857 return theAnswer;
858}
859
861 Float_t fibresRadius = -1;
862 Float_t dSiPM = -1;
863 TGeoNode* vol;
864 TGeoNode* fibre;
865 SiPMOverlap(); // 12 SiPMs per mat, made for horizontal mats, fibres staggered along y-axis.
866 auto sipm = gGeoManager->FindVolumeFast("SiPMmapVol");
867 TObjArray* Nodes = sipm->GetNodes();
868 auto plane = gGeoManager->FindVolumeFast("ScifiHorPlaneVol1");
869 for (int imat = 0; imat < plane->GetNodes()->GetEntriesFast(); imat++){
870 auto mat = static_cast<TGeoNode*>(plane->GetNodes()->At(imat));
871 Float_t t1 = mat->GetMatrix()->GetTranslation()[1];
872 auto vmat = mat->GetVolume();
873 for (int ifibre = 0; ifibre < vmat->GetNodes()->GetEntriesFast(); ifibre++){
874 fibre = static_cast<TGeoNode*>(vmat->GetNodes()->At(ifibre));
875 if (fibresRadius<0){
876 auto tmp = fibre->GetVolume()->GetShape();
877 auto S = dynamic_cast<TGeoBBox*>(tmp);
878 fibresRadius = S->GetDX();
879 }
880 Float_t t2 = fibre->GetMatrix()->GetTranslation()[1];
881 Int_t fID = fibre->GetNumber()%100000 + imat*1e4; // local fibre number, global fibre number = SO+fID
882 Float_t a = t1+t2;
883
884 // check for overlap with any of the SiPM channels in the same mat
885 for(Int_t nChan = 0; nChan< Nodes->GetEntriesFast();nChan++){ // 12 SiPMs total and 4 SiPMs per mat times 128 channels
886 vol = static_cast<TGeoNode*>(Nodes->At(nChan));
887 Int_t N = vol->GetNumber()%100000;
888 if (imat!=int(N/10000)){continue;}
889 Float_t xcentre = vol->GetMatrix()->GetTranslation()[1];
890 if (dSiPM<0){
891 TGeoBBox* B = dynamic_cast<TGeoBBox*>(vol->GetVolume()->GetShape());
892 dSiPM = B->GetDY();
893 }
894 if (TMath::Abs(xcentre-a)>4*fibresRadius){ continue;} // no need to check further
895 Float_t W = area(a,fibresRadius,xcentre-dSiPM,xcentre+dSiPM);
896 if (W<0){ continue;}
897 std::array<float, 2> Wa;
898 Wa[0] = W;
899 Wa[1] = a;
900 fibresSiPM[N][fID] = Wa;
901 }
902 }
903 }
904 // calculate also local SiPM positions based on fibre positions and their fraction
905 // probably an overkill, maximum difference between weighted average and central position < 6 micron.
906 std::map<Int_t,std::map<Int_t,std::array<float, 2>>>::iterator it;
907 std::map<Int_t,std::array<float, 2>>::iterator itx;
908 for (it = fibresSiPM.begin(); it != fibresSiPM.end(); it++)
909 {
910 Int_t N = it->first;
911 Float_t m = 0;
912 Float_t w = 0;
913 for (itx = it->second.begin(); itx != it->second.end(); itx++)
914 {
915 m+=(itx->second)[0]*(itx->second)[1];
916 w+=(itx->second)[0];
917 }
918 SiPMPos[N]=m/w;
919 }
920// make inverse mapping, which fibre is associated to which SiPMs
921 for (it = fibresSiPM.begin(); it != fibresSiPM.end(); it++)
922 {
923 Int_t N = it->first;
924 for (itx = it->second.begin(); itx != it->second.end(); itx++)
925 {
926 Int_t nfibre = itx->first;
927 siPMFibres[nfibre][N]=itx->second;
928 }
929 }
930}
932{
933 fScifiPointCollection->Clear();
934}
935
936
938{
939
946 FairRootManager::Instance()->Register("ScifiPoint", "Scifi",
947 fScifiPointCollection, kTRUE);
948}
949
950TClonesArray* Scifi::GetCollection(Int_t iColl) const
951{
952 if (iColl == 0) { return fScifiPointCollection; }
953 else { return NULL; }
954}
955
957{
958 fScifiPointCollection->Clear();
959}
960
961
962ScifiPoint* Scifi::AddHit(Int_t trackID, Int_t detID,
963 TVector3 pos, TVector3 mom,
964 Double_t time, Double_t length,
965 Double_t eLoss, Int_t pdgCode)
966{
967 TClonesArray& clref = *fScifiPointCollection;
968 Int_t size = clref.GetEntriesFast();
969 return new(clref[size]) ScifiPoint(trackID, detID, pos, mom,
970 time, length, eLoss, pdgCode);
971}
972
Double_t m
@ kLHCScifi
Definition Scifi.h:20
TVector3 GetLocalPos(Int_t id, TVector3 *glob)
Definition Scifi.cxx:709
Double32_t fTime
momentum at entrance
Definition Scifi.h:111
std::map< Int_t, float > SiPMPos
inverse mapping
Definition Scifi.h:116
ScifiPoint * AddHit(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, Double_t time, Double_t length, Double_t eLoss, Int_t pdgCode)
Definition Scifi.cxx:962
void GetSiPMPosition(Int_t SiPMChan, TVector3 &A, TVector3 &B)
Definition Scifi.cxx:727
Double32_t fELoss
length
Definition Scifi.h:113
TString last_time_alignment_tag
Definition Scifi.h:130
int last_run_time
Definition Scifi.h:129
virtual void EndOfEvent()
Definition Scifi.cxx:931
virtual void Register()
Definition Scifi.cxx:937
virtual TClonesArray * GetCollection(Int_t iColl) const
Definition Scifi.cxx:950
std::map< Int_t, std::map< Int_t, std::array< float, 2 > > > siPMFibres
mapping of fibres to SiPM channels
Definition Scifi.h:115
virtual Bool_t ProcessHits(FairVolume *v=0)
Definition Scifi.cxx:551
virtual void Reset()
Definition Scifi.cxx:956
std::vector< int > covered_runs_position_alignment
Definition Scifi.h:128
Scifi()
Definition Scifi.cxx:46
std::map< TString, std::vector< Float_t > > conf_vectors
Definition Scifi.h:123
void GetPosition(Int_t id, TVector3 &vLeft, TVector3 &vRight)
Definition Scifi.cxx:664
bool alignment_init
Definition Scifi.h:132
std::map< TString, Int_t > conf_ints
Definition Scifi.h:121
std::map< TString, Float_t > conf_floats
Definition Scifi.h:120
void InitEvent(SNDLHCEventHeader *e)
Definition Scifi.cxx:522
TClonesArray * fScifiPointCollection
local SiPM channel position
Definition Scifi.h:118
Double_t ycross(Double_t a, Double_t R, Double_t x)
Definition Scifi.cxx:818
Double_t integralSqrt(Double_t ynorm)
Definition Scifi.cxx:825
int last_run_pos
Definition Scifi.h:129
virtual void Initialize()
Definition Scifi.cxx:92
std::map< Int_t, std::map< Int_t, std::array< float, 2 > > > fibresSiPM
energy loss
Definition Scifi.h:114
Double_t fraction(Double_t R, Double_t x, Double_t y)
Definition Scifi.cxx:830
virtual void SiPMOverlap()
Definition Scifi.cxx:479
Int_t InitMedium(const char *name)
Definition Scifi.cxx:98
std::vector< int > covered_runs_time_alignment
Definition Scifi.h:127
TString last_position_alignment_tag
Definition Scifi.h:131
TLorentzVector fPos
volume id
Definition Scifi.h:109
Double32_t fLength
time
Definition Scifi.h:112
void ConstructGeometry()
Definition Scifi.cxx:118
Double_t GetCorrectedTime(Int_t fDetectorID, Double_t rawTime, Double_t L)
Definition Scifi.cxx:614
virtual ~Scifi()
Definition Scifi.cxx:84
void SiPMmapping()
Definition Scifi.cxx:860
ClassDef(Scifi, 3) private Int_t fVolumeID
track index
Definition Scifi.h:100
TLorentzVector fMom
position at entrance
Definition Scifi.h:110
SNDLHCEventHeader * eventHeader
Definition Scifi.h:124
Double_t area(Double_t a, Double_t R, Double_t xL, Double_t xR)
Definition Scifi.cxx:837
ClassImp(ecalContFact) ecalContFact