From 6e950880f3a20e0f21f28e3127ba2506c5ddf2db Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Tue, 2 Dec 2014 12:11:58 +0000
Subject: [PATCH] * Adding Fatima detector from Marc Labiche
---
Inputs/EnergyLoss/Li11[0.0]_CD2.G4table | 574 ++++-----
Inputs/EnergyLoss/Li11[0.0]_Vacuum.G4table | 376 +++---
NPLib/Fatima/Fatima.cxx | 822 +++++++++++++
NPLib/Fatima/Fatima.h | 166 +++
NPLib/Fatima/Makefile | 42 +
NPLib/Fatima/TFatimaData.cxx | 107 ++
NPLib/Fatima/TFatimaData.h | 243 ++++
NPLib/Fatima/TFatimaPhysics.cxx | 356 ++++++
NPLib/Fatima/TFatimaPhysics.h | 91 ++
NPLib/VDetector/DetectorManager.cxx | 24 +-
NPSimulation/Fatima/Fatima.cc | 133 +++
NPSimulation/Fatima/Fatima.hh | 73 ++
NPSimulation/Fatima/FatimaDetector.cc | 1237 ++++++++++++++++++++
NPSimulation/Fatima/FatimaDetector.hh | 191 +++
NPSimulation/Fatima/FatimaModule.cc | 60 +
NPSimulation/Fatima/FatimaModule.hh | 94 ++
NPSimulation/Fatima/FatimaScorers.cc | 514 ++++++++
NPSimulation/Fatima/FatimaScorers.hh | 198 ++++
NPSimulation/src/DetectorConstruction.cc | 29 +-
19 files changed, 4851 insertions(+), 479 deletions(-)
create mode 100755 NPLib/Fatima/Fatima.cxx
create mode 100755 NPLib/Fatima/Fatima.h
create mode 100755 NPLib/Fatima/Makefile
create mode 100755 NPLib/Fatima/TFatimaData.cxx
create mode 100755 NPLib/Fatima/TFatimaData.h
create mode 100755 NPLib/Fatima/TFatimaPhysics.cxx
create mode 100755 NPLib/Fatima/TFatimaPhysics.h
create mode 100755 NPSimulation/Fatima/Fatima.cc
create mode 100755 NPSimulation/Fatima/Fatima.hh
create mode 100755 NPSimulation/Fatima/FatimaDetector.cc
create mode 100755 NPSimulation/Fatima/FatimaDetector.hh
create mode 100755 NPSimulation/Fatima/FatimaModule.cc
create mode 100755 NPSimulation/Fatima/FatimaModule.hh
create mode 100755 NPSimulation/Fatima/FatimaScorers.cc
create mode 100755 NPSimulation/Fatima/FatimaScorers.hh
diff --git a/Inputs/EnergyLoss/Li11[0.0]_CD2.G4table b/Inputs/EnergyLoss/Li11[0.0]_CD2.G4table
index e0692575d..dac0963c5 100644
--- a/Inputs/EnergyLoss/Li11[0.0]_CD2.G4table
+++ b/Inputs/EnergyLoss/Li11[0.0]_CD2.G4table
@@ -1,288 +1,288 @@
Table from Geant4 generate using NPSimulation Particle: Li11[0.0] Material: CD2
-1e-09 0.000597398
-0.001 0.0897023
-0.002 0.0934766
-0.003 0.0959707
-0.004 0.0978906
-0.005 0.0998459
-0.006 0.101762
-0.007 0.103931
-0.008 0.105971
-0.009 0.107957
-0.01 0.11012
-0.011 0.112107
-0.012 0.111561
-0.022 0.110412
-0.032 0.112202
-0.042 0.114378
-0.052 0.117535
-0.062 0.120816
-0.072 0.123961
-0.082 0.126981
-0.092 0.130046
-0.102 0.133266
-0.112 0.136361
-0.122 0.139416
-0.132 0.142438
-0.142 0.145373
-0.152 0.14825
-0.162 0.151072
-0.172 0.153798
-0.182 0.156454
-0.192 0.160143
-0.202 0.164345
-0.212 0.168566
-0.222 0.172794
-0.232 0.177022
-0.242 0.181269
-0.252 0.185542
-0.262 0.189816
-0.272 0.194093
-0.282 0.198388
-0.292 0.202307
-0.302 0.205998
-0.312 0.209665
-0.322 0.213291
-0.332 0.216875
-0.342 0.220418
-0.352 0.223919
-0.362 0.227379
-0.372 0.230796
-0.382 0.234172
-0.392 0.237507
-0.402 0.240799
-0.412 0.244051
-0.422 0.24727
-0.432 0.250446
-0.442 0.253581
-0.452 0.256688
-0.462 0.259761
-0.472 0.262792
-0.482 0.265782
-0.492 0.26873
-0.502 0.271637
-0.512 0.274504
-0.522 0.27733
-0.532 0.280114
-0.542 0.282859
-0.642 0.308222
-0.742 0.32988
-0.842 0.348068
-0.942 0.363145
-1.042 0.375432
-1.142 0.385231
-1.242 0.392868
-1.342 0.398623
-1.442 0.40275
-1.542 0.405491
-2.542 0.391306
-3.542 0.35674
-4.542 0.32543
-5.542 0.299916
-6.542 0.279077
-7.542 0.261662
-8.542 0.246753
-9.542 0.23374
-10.542 0.222201
-11.542 0.211845
-12.542 0.202465
-13.542 0.193908
-14.542 0.18606
-15.542 0.178817
-16.542 0.172114
-17.542 0.16589
-18.542 0.160092
-19.542 0.154679
-20.542 0.149611
-21.542 0.144859
-22.542 0.140602
-23.542 0.136696
-24.542 0.132988
-25.542 0.129463
-26.542 0.12611
-27.542 0.123087
-28.542 0.120035
-29.542 0.117125
-30.542 0.11435
-31.542 0.111701
-32.542 0.109171
-33.542 0.106753
-34.542 0.104434
-35.542 0.102218
-36.542 0.100096
-37.542 0.0980609
-38.542 0.0961093
-39.542 0.094236
-40.542 0.0924369
-41.542 0.0907079
-42.542 0.0890451
-43.542 0.087445
-44.542 0.0859042
-45.542 0.084475
-46.542 0.0830402
-47.542 0.0816562
-48.542 0.0803203
-49.542 0.0790301
-50.542 0.0777835
-51.542 0.0765782
-52.542 0.0754123
-53.542 0.0742839
-54.542 0.0731912
-55.542 0.0721327
-56.542 0.0711068
-57.542 0.0701119
-58.542 0.0691468
-59.542 0.06821
-60.542 0.0673006
-61.542 0.0664173
-62.542 0.0655588
-63.542 0.0647574
-64.542 0.0639441
-65.542 0.0631528
-66.542 0.0623826
-67.542 0.0616325
-68.542 0.060902
-69.542 0.0601902
-70.542 0.0594963
-71.542 0.0588198
-72.542 0.0581599
-73.542 0.057516
-74.542 0.0568876
-75.542 0.0562741
-76.542 0.055675
-77.542 0.0550896
-78.542 0.0545177
-79.542 0.0539587
-80.542 0.0534121
-81.542 0.0528776
-82.542 0.0523548
-92.542 0.0476912
-102.542 0.043858
-112.542 0.0406452
-122.542 0.0379106
-132.542 0.0355525
-142.542 0.0334963
-152.542 0.0316862
-162.542 0.0300796
-172.542 0.0286429
-182.542 0.0273502
-192.542 0.0261802
-202.542 0.0251158
-212.542 0.024143
-222.542 0.0232503
-232.542 0.0224278
-242.542 0.0216675
-252.542 0.0209624
-262.542 0.0203066
-272.542 0.0196949
-282.542 0.019123
-292.542 0.018587
-302.542 0.0180836
-312.542 0.0176098
-322.542 0.0171631
-332.542 0.0167412
-342.542 0.0163419
-352.542 0.0159635
-362.542 0.0156043
-372.542 0.0152629
-382.542 0.014938
-392.542 0.0146284
-402.542 0.014333
-412.542 0.0140508
-422.542 0.0137809
-432.542 0.0135227
-442.542 0.0132752
-452.542 0.0130378
-462.542 0.01281
-472.542 0.0125912
-482.542 0.0123807
-492.542 0.0121781
-502.542 0.0119831
-512.542 0.0117951
-522.542 0.0116138
-532.542 0.0114387
-542.542 0.0112697
-552.542 0.0111064
-562.542 0.0109485
-572.542 0.0107957
-582.542 0.0106478
-592.542 0.0105045
-602.542 0.0103657
-612.542 0.0102311
-622.542 0.0101006
-632.542 0.00997386
-642.542 0.00985083
-652.542 0.00973131
-662.542 0.00961515
-672.542 0.00950222
-682.542 0.00939237
-692.542 0.0092855
-702.542 0.00918146
-712.542 0.00908016
-722.542 0.00898147
-732.542 0.0088853
-742.542 0.00879156
-752.542 0.00870014
-762.542 0.00861097
-772.542 0.00852396
-782.542 0.00843903
-882.542 0.0076894
-982.542 0.00708404
-1082.54 0.00658454
-1182.54 0.00616516
-1282.54 0.00580797
-1382.54 0.00549998
-1482.54 0.00523163
-1582.54 0.00499573
-1682.54 0.00478669
-1782.54 0.00460016
-1882.54 0.00443271
-1982.54 0.00428155
-2082.54 0.00414442
-2182.54 0.00401947
-2282.54 0.00390514
-2382.54 0.00380014
-2482.54 0.00370341
-2582.54 0.003614
-2682.54 0.00353112
-2782.54 0.00345409
-2882.54 0.00338233
-2982.54 0.00331531
-3082.54 0.00325261
-3182.54 0.00319381
-3282.54 0.00313858
-3382.54 0.0030866
-3482.54 0.00303761
-3582.54 0.00299135
-3682.54 0.00294762
-3782.54 0.00290622
-3882.54 0.00286697
-3982.54 0.00282972
-4082.54 0.00279432
-4182.54 0.00276065
-4282.54 0.00272858
-4382.54 0.00269802
-4482.54 0.00266885
-4582.54 0.002641
-4682.54 0.00261438
-4782.54 0.00258891
-5782.54 0.00238473
-6782.54 0.00224434
-7782.54 0.00214332
-8782.54 0.00206825
-9782.54 0.00200823
-10782.5 0.00196053
-11782.5 0.00192252
-12782.5 0.00189187
-13782.5 0.00186692
-14782.5 0.00184643
-15782.5 0.00182953
-25782.5 0.00175971
-35782.5 0.001755
-135783 0.00189001
-235783 0.00196877
-335783 0.00201955
-435783 0.00205667
-535783 0.00208582
-635783 0.00210979
-735783 0.00213012
+1e-09 0.000597413
+0.001 0.0897045
+0.002 0.0934789
+0.003 0.095973
+0.004 0.097893
+0.005 0.0998483
+0.006 0.101764
+0.007 0.103933
+0.008 0.105974
+0.009 0.107959
+0.01 0.110122
+0.011 0.11211
+0.012 0.111563
+0.022 0.110415
+0.032 0.112205
+0.042 0.11438
+0.052 0.117537
+0.062 0.120819
+0.072 0.123964
+0.082 0.126984
+0.092 0.130049
+0.102 0.133269
+0.112 0.136364
+0.122 0.139419
+0.132 0.14244
+0.142 0.145376
+0.152 0.148253
+0.162 0.151075
+0.172 0.153801
+0.182 0.156457
+0.192 0.160146
+0.202 0.164348
+0.212 0.168569
+0.222 0.172798
+0.232 0.177026
+0.242 0.181273
+0.252 0.185545
+0.262 0.189819
+0.272 0.194096
+0.282 0.198391
+0.292 0.202311
+0.302 0.206002
+0.312 0.209669
+0.322 0.213295
+0.332 0.216879
+0.342 0.220422
+0.352 0.223923
+0.362 0.227383
+0.372 0.230801
+0.382 0.234177
+0.392 0.237511
+0.402 0.240804
+0.412 0.244056
+0.422 0.247274
+0.432 0.250451
+0.442 0.253586
+0.452 0.256693
+0.462 0.259766
+0.472 0.262797
+0.482 0.265787
+0.492 0.268735
+0.502 0.271643
+0.512 0.274509
+0.522 0.277335
+0.532 0.28012
+0.542 0.282864
+0.642 0.308228
+0.742 0.329886
+0.842 0.348074
+0.942 0.363152
+1.042 0.375439
+1.142 0.385238
+1.242 0.392876
+1.342 0.39863
+1.442 0.402758
+1.542 0.405499
+2.542 0.391313
+3.542 0.356747
+4.542 0.325436
+5.542 0.299922
+6.542 0.279083
+7.542 0.261667
+8.542 0.246758
+9.542 0.233745
+10.542 0.222206
+11.542 0.211849
+12.542 0.202469
+13.542 0.193911
+14.542 0.186064
+15.542 0.17882
+16.542 0.172118
+17.542 0.165893
+18.542 0.160095
+19.542 0.154682
+20.542 0.149614
+21.542 0.144862
+22.542 0.140605
+23.542 0.136699
+24.542 0.13299
+25.542 0.129465
+26.542 0.126112
+27.542 0.12309
+28.542 0.120037
+29.542 0.117127
+30.542 0.114352
+31.542 0.111704
+32.542 0.109173
+33.542 0.106755
+34.542 0.104436
+35.542 0.10222
+36.542 0.100097
+37.542 0.0980627
+38.542 0.0961111
+39.542 0.0942378
+40.542 0.0924386
+41.542 0.0907096
+42.542 0.0890468
+43.542 0.0874466
+44.542 0.0859058
+45.542 0.0844765
+46.542 0.0830418
+47.542 0.0816577
+48.542 0.0803218
+49.542 0.0790316
+50.542 0.0777849
+51.542 0.0765796
+52.542 0.0754137
+53.542 0.0742853
+54.542 0.0731926
+55.542 0.0721341
+56.542 0.0711081
+57.542 0.0701132
+58.542 0.069148
+59.542 0.0682113
+60.542 0.0673019
+61.542 0.0664185
+62.542 0.06556
+63.542 0.0647587
+64.542 0.0639453
+65.542 0.063154
+66.542 0.0623837
+67.542 0.0616337
+68.542 0.0609031
+69.542 0.0601913
+70.542 0.0594974
+71.542 0.0588209
+72.542 0.058161
+73.542 0.0575171
+74.542 0.0568887
+75.542 0.0562752
+76.542 0.055676
+77.542 0.0550907
+78.542 0.0545187
+79.542 0.0539597
+80.542 0.0534131
+81.542 0.0528786
+82.542 0.0523558
+92.542 0.047692
+102.542 0.0438588
+112.542 0.040646
+122.542 0.0379113
+132.542 0.0355532
+142.542 0.0334969
+152.542 0.0316868
+162.542 0.0300801
+172.542 0.0286435
+182.542 0.0273507
+192.542 0.0261807
+202.542 0.0251163
+212.542 0.0241435
+222.542 0.0232507
+232.542 0.0224282
+242.542 0.0216679
+252.542 0.0209628
+262.542 0.020307
+272.542 0.0196952
+282.542 0.0191233
+292.542 0.0185873
+302.542 0.0180839
+312.542 0.0176102
+322.542 0.0171635
+332.542 0.0167415
+342.542 0.0163422
+352.542 0.0159638
+362.542 0.0156046
+372.542 0.0152632
+382.542 0.0149383
+392.542 0.0146287
+402.542 0.0143332
+412.542 0.014051
+422.542 0.0137812
+432.542 0.0135229
+442.542 0.0132754
+452.542 0.0130381
+462.542 0.0128103
+472.542 0.0125914
+482.542 0.0123809
+492.542 0.0121784
+502.542 0.0119833
+512.542 0.0117953
+522.542 0.011614
+532.542 0.011439
+542.542 0.0112699
+552.542 0.0111066
+562.542 0.0109487
+572.542 0.0107959
+582.542 0.010648
+592.542 0.0105047
+602.542 0.0103659
+612.542 0.0102313
+622.542 0.0101008
+632.542 0.00997405
+642.542 0.00985101
+652.542 0.00973149
+662.542 0.00961533
+672.542 0.00950239
+682.542 0.00939255
+692.542 0.00928567
+702.542 0.00918164
+712.542 0.00908033
+722.542 0.00898164
+732.542 0.00888547
+742.542 0.00879172
+752.542 0.00870031
+762.542 0.00861113
+772.542 0.00852412
+782.542 0.00843919
+882.542 0.00768954
+982.542 0.00708417
+1082.54 0.00658466
+1182.54 0.00616527
+1282.54 0.00580808
+1382.54 0.00550008
+1482.54 0.00523173
+1582.54 0.00499582
+1682.54 0.00478678
+1782.54 0.00460024
+1882.54 0.0044328
+1982.54 0.00428164
+2082.54 0.0041445
+2182.54 0.00401954
+2282.54 0.00390521
+2382.54 0.00380022
+2482.54 0.00370348
+2582.54 0.00361407
+2682.54 0.00353119
+2782.54 0.00345416
+2882.54 0.00338239
+2982.54 0.00331537
+3082.54 0.00325267
+3182.54 0.00319387
+3282.54 0.00313864
+3382.54 0.00308666
+3482.54 0.00303766
+3582.54 0.00299141
+3682.54 0.00294768
+3782.54 0.00290627
+3882.54 0.00286703
+3982.54 0.00282977
+4082.54 0.00279437
+4182.54 0.0027607
+4282.54 0.00272863
+4382.54 0.00269807
+4482.54 0.0026689
+4582.54 0.00264105
+4682.54 0.00261443
+4782.54 0.00258896
+5782.54 0.00238478
+6782.54 0.00224438
+7782.54 0.00214336
+8782.54 0.00206829
+9782.54 0.00200827
+10782.5 0.00196056
+11782.5 0.00192256
+12782.5 0.00189191
+13782.5 0.00186695
+14782.5 0.00184647
+15782.5 0.00182956
+25782.5 0.00175975
+35782.5 0.00175503
+135783 0.00189004
+235783 0.0019688
+335783 0.00201959
+435783 0.00205671
+535783 0.00208586
+635783 0.00210983
+735783 0.00213016
diff --git a/Inputs/EnergyLoss/Li11[0.0]_Vacuum.G4table b/Inputs/EnergyLoss/Li11[0.0]_Vacuum.G4table
index 653d70581..036450a3e 100644
--- a/Inputs/EnergyLoss/Li11[0.0]_Vacuum.G4table
+++ b/Inputs/EnergyLoss/Li11[0.0]_Vacuum.G4table
@@ -1,34 +1,34 @@
Table from Geant4 generate using NPSimulation Particle: Li11[0.0] Material: Vacuum
-1e-09 3.78241e-16
-0.001 5.58651e-14
-0.002 5.8415e-14
-0.003 5.96957e-14
-0.004 6.09061e-14
-0.005 6.18885e-14
-0.006 6.29769e-14
-0.007 6.40418e-14
-0.008 6.52294e-14
-0.009 6.63049e-14
-0.01 6.72932e-14
-0.011 6.84288e-14
-0.012 6.81706e-14
-0.022 6.74134e-14
-0.032 6.83174e-14
-0.042 7.00442e-14
-0.052 7.18097e-14
-0.062 7.41022e-14
-0.072 7.62818e-14
-0.082 7.83151e-14
-0.092 8.02236e-14
-0.102 8.212e-14
-0.112 8.41639e-14
-0.122 8.6198e-14
-0.132 8.81534e-14
-0.142 9.01199e-14
-0.152 9.20427e-14
-0.162 9.39088e-14
-0.172 9.57495e-14
-0.182 9.75463e-14
+1e-09 3.78244e-16
+0.001 5.58655e-14
+0.002 5.84153e-14
+0.003 5.9696e-14
+0.004 6.09065e-14
+0.005 6.18889e-14
+0.006 6.29773e-14
+0.007 6.40422e-14
+0.008 6.52297e-14
+0.009 6.63053e-14
+0.01 6.72936e-14
+0.011 6.84292e-14
+0.012 6.81709e-14
+0.022 6.74138e-14
+0.032 6.83178e-14
+0.042 7.00446e-14
+0.052 7.181e-14
+0.062 7.41025e-14
+0.072 7.62821e-14
+0.082 7.83155e-14
+0.092 8.0224e-14
+0.102 8.21204e-14
+0.112 8.41643e-14
+0.122 8.61984e-14
+0.132 8.81538e-14
+0.142 9.01203e-14
+0.152 9.20431e-14
+0.162 9.39092e-14
+0.172 9.57499e-14
+0.182 9.75467e-14
0.192 1.00014e-13
0.202 1.02779e-13
0.212 1.05567e-13
@@ -36,187 +36,187 @@ Table from Geant4 generate using NPSimulation Particle: Li11[0.0] Material: Vac
0.232 1.11175e-13
0.242 1.13988e-13
0.252 1.16806e-13
-0.262 1.19629e-13
+0.262 1.1963e-13
0.272 1.22484e-13
0.282 1.25349e-13
-0.292 1.28018e-13
-0.302 1.30514e-13
-0.312 1.32987e-13
-0.322 1.35438e-13
-0.332 1.37873e-13
+0.292 1.28019e-13
+0.302 1.30515e-13
+0.312 1.32988e-13
+0.322 1.35439e-13
+0.332 1.37874e-13
0.342 1.40303e-13
-0.352 1.42709e-13
-0.362 1.45093e-13
-0.372 1.47463e-13
+0.352 1.4271e-13
+0.362 1.45094e-13
+0.372 1.47464e-13
0.382 1.49811e-13
-0.392 1.52135e-13
+0.392 1.52136e-13
0.402 1.54437e-13
-0.412 1.56715e-13
-0.422 1.58971e-13
+0.412 1.56716e-13
+0.422 1.58972e-13
0.432 1.61204e-13
-0.442 1.63413e-13
+0.442 1.63414e-13
0.452 1.656e-13
-0.462 1.67763e-13
-0.472 1.69904e-13
-0.482 1.72021e-13
-0.492 1.74116e-13
-0.502 1.76199e-13
+0.462 1.67764e-13
+0.472 1.69905e-13
+0.482 1.72022e-13
+0.492 1.74117e-13
+0.502 1.762e-13
0.512 1.78263e-13
-0.522 1.80303e-13
-0.532 1.8232e-13
-0.542 1.84315e-13
-0.552 1.86293e-13
-0.562 1.88248e-13
-0.572 1.90181e-13
-0.582 1.92091e-13
-0.592 1.93979e-13
-0.692 2.11679e-13
+0.522 1.80304e-13
+0.532 1.82321e-13
+0.542 1.84316e-13
+0.552 1.86294e-13
+0.562 1.88249e-13
+0.572 1.90182e-13
+0.582 1.92092e-13
+0.592 1.9398e-13
+0.692 2.1168e-13
0.792 2.2734e-13
-0.892 2.41034e-13
-0.992 2.52906e-13
-1.092 2.63151e-13
-1.192 2.71897e-13
-1.292 2.79275e-13
-1.392 2.8544e-13
-1.492 2.90516e-13
-1.592 2.94617e-13
-1.692 2.97854e-13
-1.792 3.00328e-13
-2.792 2.9882e-13
-12.792 1.63184e-13
-22.792 1.17472e-13
-32.792 9.15085e-14
-42.792 7.54953e-14
-52.792 6.4383e-14
-62.792 5.62978e-14
-72.792 5.00863e-14
-82.792 4.52351e-14
-92.792 4.13342e-14
-102.792 3.80767e-14
-112.792 3.5341e-14
-122.792 3.30181e-14
-132.792 3.10022e-14
-142.792 2.92417e-14
-152.792 2.76895e-14
-162.792 2.63096e-14
-172.792 2.50746e-14
-182.792 2.39618e-14
-192.792 2.29536e-14
-202.792 2.20354e-14
-212.792 2.11955e-14
-222.792 2.04239e-14
+0.892 2.41035e-13
+0.992 2.52907e-13
+1.092 2.63152e-13
+1.192 2.71898e-13
+1.292 2.79277e-13
+1.392 2.85442e-13
+1.492 2.90517e-13
+1.592 2.94618e-13
+1.692 2.97855e-13
+1.792 3.0033e-13
+2.792 2.98821e-13
+12.792 1.63185e-13
+22.792 1.17473e-13
+32.792 9.15089e-14
+42.792 7.54956e-14
+52.792 6.43832e-14
+62.792 5.6298e-14
+72.792 5.00865e-14
+82.792 4.52353e-14
+92.792 4.13344e-14
+102.792 3.80769e-14
+112.792 3.53412e-14
+122.792 3.30183e-14
+132.792 3.10023e-14
+142.792 2.92419e-14
+152.792 2.76896e-14
+162.792 2.63097e-14
+172.792 2.50747e-14
+182.792 2.39619e-14
+192.792 2.29537e-14
+202.792 2.20355e-14
+212.792 2.11956e-14
+222.792 2.0424e-14
232.792 1.97125e-14
-242.792 1.90542e-14
-252.792 1.84433e-14
-262.792 1.78746e-14
-272.792 1.73439e-14
+242.792 1.90543e-14
+252.792 1.84434e-14
+262.792 1.78747e-14
+272.792 1.7344e-14
282.792 1.68475e-14
292.792 1.63819e-14
-302.792 1.59443e-14
-312.792 1.55323e-14
+302.792 1.59444e-14
+312.792 1.55324e-14
322.792 1.51436e-14
-332.792 1.47762e-14
+332.792 1.47763e-14
342.792 1.44285e-14
352.792 1.40988e-14
-362.792 1.37857e-14
-372.792 1.3488e-14
+362.792 1.37858e-14
+372.792 1.34881e-14
382.792 1.32046e-14
392.792 1.29344e-14
402.792 1.26765e-14
-412.792 1.243e-14
-422.792 1.21943e-14
+412.792 1.24301e-14
+422.792 1.21944e-14
432.792 1.19686e-14
-442.792 1.17522e-14
-452.792 1.15446e-14
-462.792 1.13453e-14
-472.792 1.11537e-14
+442.792 1.17523e-14
+452.792 1.15447e-14
+462.792 1.13454e-14
+472.792 1.11538e-14
482.792 1.09695e-14
492.792 1.07921e-14
502.792 1.06213e-14
512.792 1.04566e-14
522.792 1.02977e-14
532.792 1.01443e-14
-542.792 9.99614e-15
-552.792 9.85294e-15
-562.792 9.71445e-15
-572.792 9.58042e-15
-582.792 9.45065e-15
-592.792 9.32493e-15
-602.792 9.20307e-15
-612.792 9.08489e-15
-622.792 8.97023e-15
-632.792 8.85892e-15
-642.792 8.75082e-15
-652.792 8.64579e-15
-662.792 8.54371e-15
-672.792 8.44446e-15
-682.792 8.3479e-15
-692.792 8.25394e-15
-702.792 8.16246e-15
-712.792 8.07336e-15
-722.792 7.98655e-15
-732.792 7.90196e-15
-742.792 7.81948e-15
-752.792 7.73905e-15
-762.792 7.66057e-15
-772.792 7.58399e-15
-872.792 6.90856e-15
-972.792 6.36367e-15
-1072.79 5.91456e-15
-1172.79 5.53781e-15
-1272.79 5.21708e-15
-1372.79 4.94068e-15
-1472.79 4.69997e-15
-1572.79 4.48843e-15
-1672.79 4.30104e-15
-1772.79 4.13388e-15
-1872.79 3.98385e-15
-1972.79 3.84844e-15
-2072.79 3.72563e-15
-2172.79 3.61374e-15
-2272.79 3.51138e-15
-2372.79 3.4174e-15
-2472.79 3.33081e-15
-2572.79 3.25078e-15
-2672.79 3.1766e-15
-2772.79 3.10768e-15
-2872.79 3.04347e-15
-2972.79 2.98352e-15
-3072.79 2.92742e-15
-3172.79 2.87482e-15
-3272.79 2.82542e-15
-3372.79 2.77893e-15
-3472.79 2.73512e-15
-3572.79 2.69377e-15
-3672.79 2.65468e-15
-3772.79 2.61767e-15
-3872.79 2.58259e-15
-3972.79 2.5493e-15
-4072.79 2.51767e-15
-4172.79 2.48759e-15
-4272.79 2.45894e-15
-4372.79 2.43163e-15
-4472.79 2.40558e-15
-4572.79 2.3807e-15
-4672.79 2.35692e-15
-5672.79 2.16711e-15
-6672.79 2.03747e-15
-7672.79 1.94471e-15
-8672.79 1.8761e-15
-9672.79 1.82411e-15
-10672.8 1.784e-15
-11672.8 1.75265e-15
-12672.8 1.72791e-15
-13672.8 1.70826e-15
-14672.8 1.69261e-15
-24672.8 1.64202e-15
-34672.8 1.65634e-15
+542.792 9.99618e-15
+552.792 9.85298e-15
+562.792 9.71449e-15
+572.792 9.58047e-15
+582.792 9.45069e-15
+592.792 9.32497e-15
+602.792 9.20311e-15
+612.792 9.08493e-15
+622.792 8.97026e-15
+632.792 8.85896e-15
+642.792 8.75086e-15
+652.792 8.64583e-15
+662.792 8.54375e-15
+672.792 8.44449e-15
+682.792 8.34794e-15
+692.792 8.25397e-15
+702.792 8.16249e-15
+712.792 8.07339e-15
+722.792 7.98659e-15
+732.792 7.90199e-15
+742.792 7.81952e-15
+752.792 7.73908e-15
+762.792 7.66061e-15
+772.792 7.58402e-15
+872.792 6.90859e-15
+972.792 6.36369e-15
+1072.79 5.91458e-15
+1172.79 5.53783e-15
+1272.79 5.2171e-15
+1372.79 4.9407e-15
+1472.79 4.69999e-15
+1572.79 4.48845e-15
+1672.79 4.30105e-15
+1772.79 4.1339e-15
+1872.79 3.98387e-15
+1972.79 3.84846e-15
+2072.79 3.72565e-15
+2172.79 3.61376e-15
+2272.79 3.5114e-15
+2372.79 3.41741e-15
+2472.79 3.33082e-15
+2572.79 3.25079e-15
+2672.79 3.17662e-15
+2772.79 3.10769e-15
+2872.79 3.04348e-15
+2972.79 2.98353e-15
+3072.79 2.92743e-15
+3172.79 2.87483e-15
+3272.79 2.82543e-15
+3372.79 2.77895e-15
+3472.79 2.73514e-15
+3572.79 2.69378e-15
+3672.79 2.65469e-15
+3772.79 2.61768e-15
+3872.79 2.5826e-15
+3972.79 2.54932e-15
+4072.79 2.51769e-15
+4172.79 2.4876e-15
+4272.79 2.45895e-15
+4372.79 2.43164e-15
+4472.79 2.40559e-15
+4572.79 2.38071e-15
+4672.79 2.35693e-15
+5672.79 2.16712e-15
+6672.79 2.03748e-15
+7672.79 1.94472e-15
+8672.79 1.87611e-15
+9672.79 1.82412e-15
+10672.8 1.78401e-15
+11672.8 1.75266e-15
+12672.8 1.72792e-15
+13672.8 1.70827e-15
+14672.8 1.69262e-15
+24672.8 1.64203e-15
+34672.8 1.65635e-15
134673 1.90204e-15
-234673 2.04027e-15
-334673 2.13244e-15
-434673 2.20135e-15
-534673 2.25631e-15
-634673 2.30197e-15
-734673 2.34099e-15
-834673 2.37504e-15
-934673 2.40524e-15
+234673 2.04028e-15
+334673 2.13245e-15
+434673 2.20136e-15
+534673 2.25632e-15
+634673 2.30198e-15
+734673 2.341e-15
+834673 2.37505e-15
+934673 2.40525e-15
diff --git a/NPLib/Fatima/Fatima.cxx b/NPLib/Fatima/Fatima.cxx
new file mode 100755
index 000000000..9ca03bc56
--- /dev/null
+++ b/NPLib/Fatima/Fatima.cxx
@@ -0,0 +1,822 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/2014 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class is mainly an interface to the *
+ * TFatimaPhysics class *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "Fatima.h"
+
+// C++ headers
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <cmath>
+#include <stdlib.h>
+
+// NPL headers
+#include "RootInput.h"
+#include "RootOutput.h"
+
+// ROOT headers
+#include "TChain.h"
+
+using namespace std ;
+
+// Default Constructor
+
+Fatima::Fatima()
+{
+
+ m_NumberOfDetector = 0; //
+ m_EventData = new TFatimaData();
+ m_EventPhysics = new TFatimaPhysics();
+}
+
+
+
+Fatima::~Fatima()
+{
+
+ m_NumberOfDetector = 0;
+ delete m_EventData;
+ delete m_EventPhysics;
+}
+
+
+
+// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+void Fatima::ReadConfiguration(string Path)
+{
+ ifstream ConfigFile ;
+ ConfigFile.open(Path.c_str()) ;
+ string LineBuffer ;
+ string DataBuffer ;
+
+ // A:X1_Y1 --> X:1 Y:1
+ // B:X128_Y1 --> X:128 Y:1
+ // C:X1_Y128 --> X:1 Y:128
+ // D:X128_Y128 --> X:128 Y:128
+
+ double Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Az, Bz, Cz, Dz;
+ TVector3 A, B, C, D;
+ double Theta = 0, Phi = 0, R = 0, beta_u = 0 , beta_v = 0 , beta_w = 0;
+
+ bool check_A = false;
+ bool check_C = false;
+ bool check_B = false;
+ bool check_D = false;
+
+ bool check_Theta = false;
+ bool check_Phi = false;
+ bool check_R = false;
+ bool check_beta = false;
+
+ bool ReadingStatus = false;
+
+ bool isDetector = false;
+
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+
+ // If line is a Fatima bloc, reading toggle to true
+ // and toggle to true flags indicating which shape is treated.
+ if (LineBuffer.compare(0, 14, "FatimaDetector") == 0 ) {
+ cout << "///////////////////////" << endl;
+ cout << "Module found:" << endl;
+
+ if (LineBuffer.compare(0, 14, "FatimaDetector") == 0) isDetector = true;
+ ReadingStatus = true;
+ }
+ // Else don't toggle to Reading Block Status
+ else ReadingStatus = false;
+
+ // Reading Block
+ while (ReadingStatus) {
+ if (isDetector) { // FatimaDetector shape
+ ConfigFile >> DataBuffer ;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {
+ ConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n' );
+ }
+ // Finding another telescope (safety), toggle out
+ else if (DataBuffer.compare(0, 14, "FatimaDetector") == 0) {
+ cout << "WARNING: Another Module is find before standard sequence of Token, Error may occured in Telecope definition" << endl;
+ ReadingStatus = false;
+ }
+
+ // Position method
+ else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az ;
+
+ A = TVector3(Ax, Ay, Az);
+ cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
+ }
+ else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz ;
+
+ B = TVector3(Bx, By, Bz);
+ cout << "X128 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
+ }
+ else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz ;
+
+ C = TVector3(Cx, Cy, Cz);
+ cout << "X1 Y128 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
+ }
+ else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz ;
+
+ D = TVector3(Dx, Dy, Dz);
+ cout << "X128 Y128 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
+ } // End Position Method
+
+ // Angle method
+ else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ Theta = Theta ;
+ cout << "Theta: " << Theta << endl;
+ }
+ else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ Phi = Phi ;
+ cout << "Phi: " << Phi << endl;
+ }
+ else if (DataBuffer.compare(0, 2, "R=") == 0) {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ R = R ;
+ cout << "R: " << R << endl;
+ }
+ else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
+ check_beta = true;
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ beta_u = beta_u ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ beta_v = beta_v ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ beta_w = beta_w ;
+ cout << "Beta: " << beta_u << " " << beta_v << " " << beta_w << endl ;
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) ) {
+ ReadingStatus = false;
+
+ // Add The previously define telescope
+ // With position method
+ if ( check_A && check_B && check_C && check_D ) {
+ AddModuleDummyShape(A ,
+ B ,
+ C ,
+ D ) ;
+ }
+
+ // with angle method
+ else if ( check_Theta && check_Phi && check_R && check_beta ) {
+ AddModuleDummyShape(Theta,
+ Phi,
+ R,
+ beta_u,
+ beta_v,
+ beta_w);
+ }
+
+ // reset boolean flag for point positioning
+ check_A = false;
+ check_B = false;
+ check_C = false;
+ check_D = false;
+
+ // reset boolean flag for angle positioning
+ check_Theta = false;
+ check_Phi = false;
+ check_R = false;
+ check_beta = false;
+
+ // reset boolean flag for shape determination
+ isDetector = false;
+
+ } // end test for adding a module
+ } // end test for FatimaDetector shape
+
+
+ } // end while for reading block
+ } // end while for reading file
+
+ cout << endl << "/////////////////////////////" << endl<<endl;
+}
+
+
+// Read stream at Path and pick-up calibration parameter using Token
+// If argument is "Simulation" no change calibration is loaded
+void Fatima::ReadCalibrationFile(string Path)
+{
+ // Order of Polynom function used for calibration
+ int Calibration_Si_E_Order;
+ int Calibration_Si_T_Order;
+ int Calibration_SiLi_E_Order;
+ int Calibration_CsI_E_Order;
+
+ // Calibration_Si_X_E[DetectorNumber][StripNumber][Order of Coeff]
+ vector< vector< vector< double > > > Calibration_Si_X_E ;
+ vector< vector< vector< double > > > Calibration_Si_X_T ;
+ vector< vector< vector< double > > > Calibration_Si_Y_E ;
+ vector< vector< vector< double > > > Calibration_Si_Y_T ;
+
+ // Calibration_SiLi_E[DetectorNumber][PadNumber][Order of Coeff]
+ vector< vector< vector< double > > > Calibration_SiLi_E ;
+
+ // Calibration_SiLi_E[DetectorNumber][CrystalNumber][Order of Coeff]
+ vector< vector< vector< double > > > Calibration_CsI_E ;
+
+ if (Path == "Simulation") { // Simulation case: data already calibrated
+ Calibration_Si_E_Order = 1;
+ Calibration_Si_T_Order = 1;
+ Calibration_SiLi_E_Order = 1;
+ Calibration_CsI_E_Order = 1;
+
+ vector<double> Coef;
+ // Order 0 Order 1
+ Coef.push_back(0) ; Coef.push_back(1) ;
+
+ vector< vector<double> > StripLine ;
+ StripLine.resize( 128 , Coef) ;
+
+ Calibration_Si_X_E.resize( m_NumberOfDetector , StripLine) ;
+ Calibration_Si_X_T.resize( m_NumberOfDetector , StripLine) ;
+ Calibration_Si_Y_E.resize( m_NumberOfDetector , StripLine) ;
+ Calibration_Si_Y_T.resize( m_NumberOfDetector , StripLine) ;
+
+ Calibration_SiLi_E.resize( m_NumberOfDetector , StripLine) ;
+ Calibration_CsI_E .resize( m_NumberOfDetector , StripLine) ;
+ }
+ else {
+ }
+}
+
+
+
+// Activated associated Branches and link it to the private member DetectorData address
+// In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
+void Fatima::InitializeRootInput()
+{
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("FATIMA", true);
+ inputChain->SetBranchStatus("fFATIMA*", true);
+ inputChain->SetBranchAddress("FATIMA", &m_EventData);
+}
+
+
+
+// Create associated branches and associated private member DetectorPhysics address
+void Fatima::InitializeRootOutput()
+{
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("FATIMA", "TFatimaPhysics", &m_EventPhysics);
+}
+
+
+
+// This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
+void Fatima::BuildPhysicalEvent()
+{
+ m_EventPhysics -> BuildPhysicalEvent(m_EventData);
+}
+
+
+
+// Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
+// This method aimed to be used for analysis performed during experiment, when speed is requiered.
+// NB: This method can eventually be the same as BuildPhysicalEvent.
+void Fatima::BuildSimplePhysicalEvent()
+{
+ m_EventPhysics -> BuildSimplePhysicalEvent(m_EventData);
+}
+
+
+//
+// Used for 3x3 clusters (FatimaCluster):
+//
+void Fatima::AddModuleSquare(TVector3 C_X1_Y1,
+ TVector3 C_X128_Y1,
+ TVector3 C_X1_Y128,
+ TVector3 C_X128_Y128)
+{
+ m_NumberOfDetector++;
+
+ cout << " m_NumberOfDetector " << m_NumberOfDetector << endl;
+
+ // remove warning using C_X128_Y128
+ C_X128_Y128.Unit();
+
+ // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U = C_X128_Y1 - C_X1_Y1;
+ U = U.Unit();
+
+ // Vector V on Module Face (parallele to X Strip)
+ TVector3 V = C_X1_Y128 - C_X1_Y1;
+ V = V.Unit();
+
+ // Position Vector of Strip Center
+ TVector3 StripCenter = TVector3(0,0,0);
+ // Position Vector of X=1 Y=1 Strip
+ TVector3 Strip_1_1;
+
+
+ // !!!!!! for cluster 3x3 only !!!!!
+ //
+ // Geometry Parameter in case of 3x3 cluster
+ double Face = 169; // mm cf: FatimaCluster.hh
+ double NumberOfStrip = 3; // number of crystals per raw or per column
+ double StripPitch = Face/NumberOfStrip; // mm
+ //
+ //double Face = 98; // mm
+ //double NumberOfStrip = 128;
+ //double StripPitch = Face/NumberOfStrip; // mm
+ //
+ //
+
+ // Buffer object to fill Position Array
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModuleStripPositionX;
+ vector< vector< double > > OneModuleStripPositionY;
+ vector< vector< double > > OneModuleStripPositionZ;
+
+ // Moving StripCenter to 1.1 corner:
+ Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2.);
+
+ for (int i = 0; i < NumberOfStrip; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStrip; j++) {
+ StripCenter = Strip_1_1 + StripPitch*( i*U + j*V );
+// StripCenter += -TargetPosition;
+
+ lineX.push_back( StripCenter.X() );
+ lineY.push_back( StripCenter.Y() );
+ lineZ.push_back( StripCenter.Z() );
+ }
+
+ OneModuleStripPositionX.push_back(lineX);
+ OneModuleStripPositionY.push_back(lineY);
+ OneModuleStripPositionZ.push_back(lineZ);
+ }
+
+ m_DetPositionX.push_back( OneModuleStripPositionX );
+ m_DetPositionY.push_back( OneModuleStripPositionY );
+ m_DetPositionZ.push_back( OneModuleStripPositionZ );
+}
+
+
+
+void Fatima::AddModuleSquare(double theta,
+ double phi,
+ double distance,
+ double beta_u,
+ double beta_v,
+ double beta_w)
+{
+ m_NumberOfDetector++;
+
+
+ // convert from degree to radian:
+ double Pi = 3.141592654;
+ theta = theta * Pi/180. ;
+ phi = phi * Pi/180. ;
+
+ // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U ;
+ // Vector V on Module Face (parallele to X Strip)
+ TVector3 V ;
+ // Vector W normal to Module Face (pointing CsI)
+ TVector3 W ;
+ // Vector position of Module Face center
+ TVector3 C ;
+
+ C = TVector3(distance * sin(theta) * cos(phi),
+ distance * sin(theta) * sin(phi),
+ distance * cos(theta));
+
+ TVector3 YperpC = TVector3( cos(theta) * cos(phi),
+ cos(theta) * sin(phi),
+ -sin(theta));
+
+ W = C.Unit();
+ U = W.Cross(YperpC);
+ V = W.Cross(U);
+
+ U = U.Unit();
+ V = V.Unit();
+
+ U.Rotate( beta_u * Pi/180. , U ) ;
+ V.Rotate( beta_u * Pi/180. , U ) ;
+
+ U.Rotate( beta_v * Pi/180. , V ) ;
+ V.Rotate( beta_v * Pi/180. , V ) ;
+
+ U.Rotate( beta_w * Pi/180. , W ) ;
+ V.Rotate( beta_w * Pi/180. , W ) ;
+
+
+ //
+ // !!!!!! for cluster 3x3 only !!!!!
+ //
+ // Geometry Parameter in case of 3x3 cluster
+ double Face = 169; // mm cf: FatimaCluster.hh
+ double NumberOfStrip = 3; // number of crystals per raw or per column
+ double StripPitch = Face/NumberOfStrip; // mm
+ //
+ //double Face = 98; // mm
+ //double NumberOfStrip = 128;
+ //double StripPitch = Face/NumberOfStrip; // mm
+ //
+ //
+
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModuleStripPositionX;
+ vector< vector< double > > OneModuleStripPositionY;
+ vector< vector< double > > OneModuleStripPositionZ;
+
+ double X, Y, Z;
+
+ // Moving C to the 1.1 corner:
+ C.SetX( C.X() - ( Face/2 - StripPitch/2 ) * ( V.X() + U.X() ) ) ;
+ C.SetY( C.Y() - ( Face/2 - StripPitch/2 ) * ( V.Y() + U.Y() ) ) ;
+ C.SetZ( C.Z() - ( Face/2 - StripPitch/2 ) * ( V.Z() + U.Z() ) ) ;
+
+ for (int i = 0; i < NumberOfStrip; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStrip; j++) {
+ X = C.X() + StripPitch * ( U.X()*i + V.X()*j );
+ Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j );
+ Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j );
+
+ lineX.push_back(X);
+ lineY.push_back(Y);
+ lineZ.push_back(Z);
+ }
+
+ OneModuleStripPositionX.push_back(lineX);
+ OneModuleStripPositionY.push_back(lineY);
+ OneModuleStripPositionZ.push_back(lineZ);
+ }
+
+ m_DetPositionX.push_back( OneModuleStripPositionX );
+ m_DetPositionY.push_back( OneModuleStripPositionY );
+ m_DetPositionZ.push_back( OneModuleStripPositionZ );
+
+
+}
+
+
+//
+// Used for single detector (FatimaDetector):
+//
+void Fatima::AddModuleDummyShape(TVector3 C_X1_Y1,
+ TVector3 C_X128_Y1,
+ TVector3 C_X1_Y128,
+ TVector3 C_X128_Y128)
+{
+ m_NumberOfDetector++;
+
+ cout << " m_NumberOfDetector " << m_NumberOfDetector << endl;
+
+ // remove warning using C_X128_Y128
+ C_X128_Y128.Unit();
+
+ // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U = C_X128_Y1 - C_X1_Y1;
+ U = U.Unit();
+
+ // Vector V on Module Face (parallele to X Strip)
+ TVector3 V = C_X1_Y128 - C_X1_Y1;
+ V = V.Unit();
+
+ // Position Vector of Strip Center
+ TVector3 StripCenter = TVector3(0,0,0);
+ // Position Vector of X=1 Y=1 Strip
+ TVector3 Strip_1_1;
+
+ //
+ // !!!!!! for for single phoswich !!!!
+ // Geometry Parameter in case of 3x3 cluster
+ double Face = 57; // mm cf: FatimaCluster.hh
+ double NumberOfStrip = 1; // number of crystals per raw or per column
+ double StripPitch = Face/NumberOfStrip; // mm
+ //
+ //
+
+ // Buffer object to fill Position Array
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModuleStripPositionX;
+ vector< vector< double > > OneModuleStripPositionY;
+ vector< vector< double > > OneModuleStripPositionZ;
+
+ // Moving StripCenter to 1.1 corner:
+ Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2.);
+
+ for (int i = 0; i < NumberOfStrip; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStrip; j++) {
+ StripCenter = Strip_1_1 + StripPitch*( i*U + j*V );
+// StripCenter += -TargetPosition;
+
+ lineX.push_back( StripCenter.X() );
+ lineY.push_back( StripCenter.Y() );
+ lineZ.push_back( StripCenter.Z() );
+ }
+
+ OneModuleStripPositionX.push_back(lineX);
+ OneModuleStripPositionY.push_back(lineY);
+ OneModuleStripPositionZ.push_back(lineZ);
+ }
+
+ m_DetPositionX.push_back( OneModuleStripPositionX );
+ m_DetPositionY.push_back( OneModuleStripPositionY );
+ m_DetPositionZ.push_back( OneModuleStripPositionZ );
+}
+
+
+
+
+void Fatima::AddModuleDummyShape(double theta,
+ double phi,
+ double distance,
+ double beta_u,
+ double beta_v,
+ double beta_w)
+{
+ m_NumberOfDetector++;
+
+ cout << " m_NumberOfDetector " << m_NumberOfDetector << endl;
+
+ // convert from degree to radian:
+ double Pi = 3.141592654;
+ theta = theta * Pi/180. ;
+ phi = phi * Pi/180. ;
+
+ // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U ;
+ // Vector V on Module Face (parallele to X Strip)
+ TVector3 V ;
+ // Vector W normal to Module Face (pointing CsI)
+ TVector3 W ;
+ // Vector position of Module Face center
+ TVector3 C ;
+
+ C = TVector3(distance * sin(theta) * cos(phi),
+ distance * sin(theta) * sin(phi),
+ distance * cos(theta));
+
+ TVector3 YperpW = TVector3( cos(theta) * cos(phi),
+ cos(theta) * sin(phi),
+ -sin(theta));
+
+ W = C.Unit();
+ U = W.Cross(YperpW);
+ V = W.Cross(U);
+
+ U = U.Unit();
+ V = V.Unit();
+
+ U.Rotate( beta_u * Pi/180. , U ) ;
+ V.Rotate( beta_u * Pi/180. , U ) ;
+
+ U.Rotate( beta_v * Pi/180. , V ) ;
+ V.Rotate( beta_v * Pi/180. , V ) ;
+
+ U.Rotate( beta_w * Pi/180. , W ) ;
+ V.Rotate( beta_w * Pi/180. , W ) ;
+
+ double Face = 50; // mm
+ double NumberOfStrip = 25;
+ double StripPitch = Face/NumberOfStrip; // mm
+
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModuleStripPositionX;
+ vector< vector< double > > OneModuleStripPositionY;
+ vector< vector< double > > OneModuleStripPositionZ;
+
+ double X, Y, Z;
+
+ // Moving C to the 1.1 corner:
+ C.SetX( C.X() - ( Face/2 - StripPitch/2 ) * ( V.X() + U.X() ) ) ;
+ C.SetY( C.Y() - ( Face/2 - StripPitch/2 ) * ( V.Y() + U.Y() ) ) ;
+ C.SetZ( C.Z() - ( Face/2 - StripPitch/2 ) * ( V.Z() + U.Z() ) ) ;
+
+ for (int i = 0; i < NumberOfStrip; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStrip; j++) {
+ X = C.X() + StripPitch * ( U.X()*i + V.X()*j );
+ Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j );
+ Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j );
+
+ lineX.push_back(X);
+ lineY.push_back(Y);
+ lineZ.push_back(Z);
+ }
+
+ OneModuleStripPositionX.push_back(lineX);
+ OneModuleStripPositionY.push_back(lineY);
+ OneModuleStripPositionZ.push_back(lineZ);
+ }
+
+ m_DetPositionX.push_back( OneModuleStripPositionX );
+ m_DetPositionY.push_back( OneModuleStripPositionY );
+ m_DetPositionZ.push_back( OneModuleStripPositionZ );
+}
+
+
+
+double Fatima::GetEnergyDeposit()
+{
+ if (m_EventPhysics->FatimaTotalEnergy.size() > 0)
+ return m_EventPhysics->FatimaTotalEnergy[0];
+ else
+ return -1000;
+}
+
+double Fatima::GetEnergyInDeposit() // inner Layer
+{
+ if (m_EventPhysics->FatimaInTotalEnergy.size() > 0)
+ return m_EventPhysics->FatimaInTotalEnergy[0];
+ else
+ return -1000;
+}
+
+double Fatima::GetEnergyOutDeposit() // Outer Layer
+{
+ if (m_EventPhysics->FatimaOutTotalEnergy.size() > 0)
+ return m_EventPhysics->FatimaOutTotalEnergy[0];
+ else
+ return -1000;
+}
+
+
+TVector3 Fatima::GetPositionOfInteraction()
+{
+ TVector3 Position = TVector3(-1000,-1000,-1000);
+
+
+ if(m_EventPhysics->DetectorNumber[0]>100 ){ // ie: cluster
+
+ cout << "Cluster#: " << m_EventPhysics->DetectorNumber[0] << endl;
+ cout << "Crystal column (x) #: " << m_EventPhysics->FirstStage_X[0] << endl;
+ cout << "Crystal raw (y) #: " << m_EventPhysics->FirstStage_Y[0] << endl;
+
+ if(m_EventPhysics->FirstStage_X[0]>-1 && m_EventPhysics->FirstStage_Y[0]>-1){ // ie: LaBr hit first
+ Position = TVector3(GetDetPositionX(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
+ GetDetPositionY(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
+ GetDetPositionZ(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]));
+ }/*
+ else
+ {
+ Position = TVector3(GetCrystPositionX(m_EventPhysics->ModuleNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
+ GetCrystPositionY(m_EventPhysics->ModuleNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
+ GetCrystPositionZ(m_EventPhysics->ModuleNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]));
+
+ }
+ */
+ }
+
+ if(m_EventPhysics->DetectorNumber[0]>0 ){ // ie: phoswich
+
+ cout << "Detector #: " << m_EventPhysics->DetectorNumber[0] << endl;
+ cout << "Crystal column (x) #: " << m_EventPhysics->FirstStage_X[0] << endl;
+ cout << "Crystal raw (y) #: " << m_EventPhysics->FirstStage_Y[0] << endl;
+
+
+ if(m_EventPhysics->FirstStage_X[0]>-1 && m_EventPhysics->FirstStage_Y[0]>-1){ // ie: LaBr hit first
+ Position = TVector3(GetDetPositionX(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
+ GetDetPositionY(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]),
+ GetDetPositionZ(m_EventPhysics->DetectorNumber[0], m_EventPhysics->FirstStage_X[0], m_EventPhysics->FirstStage_Y[0]));
+ }/*else
+ {
+ Position = TVector3(GetDetPositionX(m_EventPhysics->DetectorNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
+ GetDetPositionY(m_EventPhysics->DetectorNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]),
+ GetDetPositionZ(m_EventPhysics->DetectorNumber[0], m_EventPhysics->SecondStage_X[0], m_EventPhysics->SecondStage_Y[0]));
+
+ }*/
+
+ }
+
+
+ return(Position);
+}
+
+
+
+void Fatima::Print()
+{
+ cout << "Number of Detectors: " << m_NumberOfDetector << endl;
+
+ for (int f = 0; f < m_NumberOfDetector; f++) {
+ cout << "Detector " << f+1 << endl;
+
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ cout << i+1 << " "<< j+1 << " "
+ << m_DetPositionX[f][i][j] << " "
+ << m_DetPositionY[f][i][j] << " "
+ << m_DetPositionZ[f][i][j] << " "
+ << endl ;
+ }
+ }
+
+ /*
+ for (int i = 0; i < 128; i++) {
+ for (int j = 0; j < 128; j++) {
+ cout << i+1 << " "<< j+1 << " "
+ << m_StripPositionX[f][i][j] << " "
+ << m_StripPositionY[f][i][j] << " "
+ << m_StripPositionZ[f][i][j] << " "
+ << endl ;
+ }
+ }
+ */
+
+ }
+}
diff --git a/NPLib/Fatima/Fatima.h b/NPLib/Fatima/Fatima.h
new file mode 100755
index 000000000..903740901
--- /dev/null
+++ b/NPLib/Fatima/Fatima.h
@@ -0,0 +1,166 @@
+
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/2014 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class is mainly an interface to the *
+ * TFatimaPhysics class *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#ifndef Fatima_H
+#define Fatima_H
+
+// NPL
+#include "../include/VDetector.h"
+#include "TFatimaData.h"
+#include "TFatimaPhysics.h"
+
+// Root
+#include "TVector3.h"
+
+class Fatima : public NPA::VDetector
+{
+public:
+ Fatima();
+ virtual ~Fatima();
+
+public:
+ /////////////////////////////////////
+ // Innherited from VDetector Class //
+ /////////////////////////////////////
+ // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+ void ReadConfiguration(string);
+
+ // Read stream at CalibFile and pick-up calibration parameter using Token
+ // If argument is "Simulation" no change calibration is loaded
+ void ReadCalibrationFile(string);
+
+ // Activated associated Branches and link it to the private member DetectorData address
+ // In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
+ void InitializeRootInput();
+
+ // Create associated branches and associated private member DetectorPhysics address
+ void InitializeRootOutput();
+
+ // This method is called at each event read from the Input Tree.
+ // The aim is to build treat Raw dat in order to extract physical parameter.
+ void BuildPhysicalEvent();
+
+ // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
+ // This method aimed to be used for analysis performed during experiment, when speed is requiered.
+ // NB: This method can eventually be the same as BuildPhysicalEvent.
+ void BuildSimplePhysicalEvent();
+
+ // Those two method all to clear the Event Physics or Data
+ void ClearEventPhysics() {m_EventPhysics->Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
+
+public:
+ ////////////////////////////////
+ // Specific to Fatima //
+ ////////////////////////////////
+ // Case of a Square module
+ // Add a Module using Corner Coordinate information
+ void AddModuleSquare(TVector3 C_X1_Y1,
+ TVector3 C_X128_Y1,
+ TVector3 C_X1_Y128,
+ TVector3 C_X128_Y128);
+
+ // Add a Module using R Theta Phi of Si center information
+ void AddModuleSquare(double theta,
+ double phi,
+ double distance,
+ double beta_u,
+ double beta_v,
+ double beta_w);
+
+ // Case of a DummyShape module
+ // Add a Module using Corner Coordinate information
+ void AddModuleDummyShape(TVector3 C_X1_Y1,
+ TVector3 C_X128_Y1,
+ TVector3 C_X1_Y128,
+ TVector3 C_X128_Y128);
+
+ // Add a Module using R Theta Phi of Si center information
+ void AddModuleDummyShape(double theta,
+ double phi,
+ double distance,
+ double beta_u,
+ double beta_v,
+ double beta_w);
+
+ // Getters to retrieve the (X,Y,Z) coordinates of a pixel (crystal) defined by (X,Y) for Detector
+
+ double GetDetPositionX(int N ,int X ,int Y) { return m_DetPositionX[N-1][X][Y]; } // remember index=0 for Fatima Detector (cf FatimaModule.cxx)
+ double GetDetPositionY(int N ,int X ,int Y) { return m_DetPositionY[N-1][X][Y]; }
+ double GetDetPositionZ(int N ,int X ,int Y) { return m_DetPositionZ[N-1][X][Y]; }
+
+ double GetNumberOfDetector() { return m_NumberOfDetector; }
+
+ /*
+ // Getters to retrieve the (X,Y,Z) coordinates of a pixel (crystal) defined by (X,Y) for clusters
+
+ double GetCrystPositionX(int N ,int X ,int Y) { return m_CrystPositionX[N-101][X][Y]; } // remember index=100 for Fatima Cluster (cf FatimaModule.cxx)
+ double GetCrystPositionY(int N ,int X ,int Y) { return m_CrystPositionY[N-101][X][Y]; }
+ double GetCrystPositionZ(int N ,int X ,int Y) { return m_CrystPositionZ[N-101][X][Y]; }
+ */
+
+
+ //double GetStripPositionX(int N ,int X ,int Y) { return m_StripPositionX[N-1][X-1][Y-1]; }
+ //double GetStripPositionY(int N ,int X ,int Y) { return m_StripPositionY[N-1][X-1][Y-1]; }
+ //double GetStripPositionZ(int N ,int X ,int Y) { return m_StripPositionZ[N-1][X-1][Y-1]; }
+
+ //double GetNumberOfModule() { return m_NumberOfModule; }
+
+
+ // Get Root input and output objects
+ TFatimaData* GetEventData() {return m_EventData;}
+ TFatimaPhysics* GetEventPhysics() {return m_EventPhysics;}
+
+ // To be called after a build Physical Event
+ double GetEnergyDeposit();
+ double GetEnergyInDeposit();
+ double GetEnergyOutDeposit();
+ TVector3 GetPositionOfInteraction();
+
+ void Print();
+
+
+private:
+ ////////////////////////////////////////
+ // Root Input and Output tree classes //
+ ////////////////////////////////////////
+ TFatimaData* m_EventData;
+ TFatimaPhysics* m_EventPhysics;
+
+
+private:
+ // Spatial Position of Crystal Calculated on basis of detector(cluster) position
+ /*
+ int m_NumberOfModule; // cluster
+ vector< vector < vector < double > > > m_CrystPositionX;
+ vector< vector < vector < double > > > m_CrystPositionY;
+ vector< vector < vector < double > > > m_CrystPositionZ;
+ */
+ int m_NumberOfDetector; // single detector
+ vector< vector < vector < double > > > m_DetPositionX;
+ vector< vector < vector < double > > > m_DetPositionY;
+ vector< vector < vector < double > > > m_DetPositionZ;
+
+};
+
+#endif
diff --git a/NPLib/Fatima/Makefile b/NPLib/Fatima/Makefile
new file mode 100755
index 000000000..e25e9da81
--- /dev/null
+++ b/NPLib/Fatima/Makefile
@@ -0,0 +1,42 @@
+include ../Makefile.arch
+
+#------------------------------------------------------------------------------
+SHARELIB = libFatima.so
+all: $(SHARELIB)
+#------------------------------------------------------------------------------
+############### Detector ##############
+
+## Fatima ##
+libFatima.so: TFatimaData.o TFatimaDataDict.o Fatima.o TFatimaPhysics.o TFatimaPhysicsDict.o
+ $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
+
+TFatimaDataDict.cxx: TFatimaData.h
+ rootcint -f $@ -c $^
+
+TFatimaPhysicsDict.cxx: TFatimaPhysics.h
+ rootcint -f $@ -c $^
+
+# dependances
+Fatima.o: Fatima.cxx Fatima.h
+TFatimaData.o: TFatimaData.cxx TFatimaData.h
+TFatimaDataDict.o: TFatimaData.cxx TFatimaData.h
+TFatimaPhysics.o: TFatimaPhysics.cxx TFatimaPhysics.h
+TFatimaPhysicsDict.o: TFatimaPhysics.cxx TFatimaPhysics.h
+
+#######################################
+
+############# Clean and More ##########
+clean:
+ @rm -f core *~ *.o *Dict*
+
+distclean:
+ make clean; rm -f *.so
+
+.SUFFIXES: .$(SrcSuf)
+
+###
+
+.$(SrcSuf).$(ObjSuf):
+ $(CXX) $(CXXFLAGS) $(INCLUDE) -c $<
+
+
diff --git a/NPLib/Fatima/TFatimaData.cxx b/NPLib/Fatima/TFatimaData.cxx
new file mode 100755
index 000000000..cd68697d9
--- /dev/null
+++ b/NPLib/Fatima/TFatimaData.cxx
@@ -0,0 +1,107 @@
+#include <iostream>
+#include "TFatimaData.h"
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/2013 *
+ * Last update : 02/07/2014 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class described the raw data of the Fatima detector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+ClassImp(TFatimaData)
+
+TFatimaData::TFatimaData()
+{
+ // Default constructor
+ Clear();
+}
+
+
+
+TFatimaData::~TFatimaData()
+{
+}
+
+
+
+void TFatimaData::Clear()
+{
+
+ /*
+ fFatima_Energy.clear();
+ fFatima_Number.clear();
+ fFatima_Time.clear();
+ */
+
+
+ fFATIMA_LaBr3Stage_E_DetectorNbr.clear();
+ fFATIMA_LaBr3Stage_E_CrystalNbr.clear();
+ fFATIMA_LaBr3Stage_E_Energy.clear();
+ // fFATIMA_LaBr3Stage_Eff_phpeak.clear();
+ // Time
+ fFATIMA_LaBr3Stage_T_DetectorNbr.clear();
+ fFATIMA_LaBr3Stage_T_CrystalNbr.clear();
+ fFATIMA_LaBr3Stage_T_Time.clear();
+
+
+ /*
+ // Second Stage CsI
+ // CsI
+ // Energy
+ fFATIMA_CsIStage_E_DetectorNbr.clear();
+ fFATIMA_CsIStage_E_CrystalNbr.clear();
+ fFATIMA_CsIStage_E_Energy.clear();
+ // fFATIMA_CsIStage_Eff_phpeak.clear();
+ // Time
+ fFATIMA_CsIStage_T_DetectorNbr.clear();
+ fFATIMA_CsIStage_T_CrystalNbr.clear();
+ fFATIMA_CsIStage_T_Time.clear();
+ */
+}
+
+
+
+void TFatimaData::Dump() const
+{
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
+ /*
+ for (unsigned short i = 0; i<fFatima_Energy.size(); i ++) {
+ cout << "Fatima Number " << fFatima_Number[i] << " Energy: " << fFatima_Energy[i] << " Time: "<< fFatima_Time[i] << endl;
+ }
+ */
+
+ for (UShort_t i = 0; i < fFATIMA_LaBr3Stage_E_DetectorNbr.size(); i++)
+ cout << "DetNbr: " << fFATIMA_LaBr3Stage_E_DetectorNbr[i] << " Crystal: " << fFATIMA_LaBr3Stage_E_CrystalNbr[i] << " Energy: " << fFATIMA_LaBr3Stage_E_Energy[i] << endl;
+ // (X,T)
+ cout << "FATIMA_LaBr3Stage_T_Mult = " << fFATIMA_LaBr3Stage_T_DetectorNbr.size() << endl;
+ for (UShort_t i = 0; i < fFATIMA_LaBr3Stage_T_DetectorNbr.size(); i++)
+ cout << "DetNbr: " << fFATIMA_LaBr3Stage_T_DetectorNbr[i] << " Crystal: " << fFATIMA_LaBr3Stage_T_CrystalNbr[i] << " Time: " << fFATIMA_LaBr3Stage_T_Time[i] << endl;
+
+ /*
+ // Second Stage
+ // Energy
+ cout << "FATIMA_CsIStage_E_Mult = " << fFATIMA_CsIStage_E_DetectorNbr.size() << endl;
+ for (UShort_t i = 0; i < fFATIMA_CsIStage_E_DetectorNbr.size(); i++)
+ cout << "Det: " << fFATIMA_CsIStage_E_DetectorNbr[i] << " Crystal: " << fFATIMA_CsIStage_E_CrystalNbr[i] << " Energy: " << fFATIMA_CsIStage_E_Energy[i] << endl;
+ // Time
+ cout << "FATIMA_CsIStage_T_Mult = " << fFATIMA_CsIStage_T_DetectorNbr.size() << endl;
+ for (UShort_t i = 0; i < fFATIMA_CsIStage_T_DetectorNbr.size(); i++)
+ cout << "Det: " << fFATIMA_CsIStage_T_DetectorNbr[i] << " Crystal: " << fFATIMA_CsIStage_T_CrystalNbr[i] << " Time: " << fFATIMA_CsIStage_T_Time[i] << endl;
+
+ */
+
+}
diff --git a/NPLib/Fatima/TFatimaData.h b/NPLib/Fatima/TFatimaData.h
new file mode 100755
index 000000000..1019cc0f2
--- /dev/null
+++ b/NPLib/Fatima/TFatimaData.h
@@ -0,0 +1,243 @@
+#ifndef __FatimaDATA__
+#define __FatimaDATA__
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/2013 *
+ * Last update : 02/07/2014 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class described the raw data of the Fatima detector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include <vector>
+
+#include "TObject.h"
+using namespace std ;
+
+
+class TFatimaData : public TObject {
+
+ protected:
+ // First Stage LaBr
+ // LaBr3
+ // Energy
+ vector<UShort_t> fFATIMA_LaBr3Stage_E_DetectorNbr;
+ vector<UShort_t> fFATIMA_LaBr3Stage_E_CrystalNbr;
+ vector<Double_t> fFATIMA_LaBr3Stage_E_Energy;
+ // vector<Double_t> fFATIMA_LaBr3Stage_Eff_phpeak;
+ // Time
+ vector<UShort_t> fFATIMA_LaBr3Stage_T_DetectorNbr;
+ vector<UShort_t> fFATIMA_LaBr3Stage_T_CrystalNbr;
+ vector<Double_t> fFATIMA_LaBr3Stage_T_Time;
+
+
+ /*
+ // Second Stage CsI
+ // CsI
+ // Energy
+ vector<UShort_t> fFATIMA_CsIStage_E_DetectorNbr;
+ vector<UShort_t> fFATIMA_CsIStage_E_CrystalNbr;
+ vector<Double_t> fFATIMA_CsIStage_E_Energy;
+ // vector<Double_t> fFATIMA_CsIStage_Eff_phpeak;
+ // Time
+ vector<UShort_t> fFATIMA_CsIStage_T_DetectorNbr;
+ vector<UShort_t> fFATIMA_CsIStage_T_CrystalNbr;
+ vector<Double_t> fFATIMA_CsIStage_T_Time;
+ */
+
+ /*
+ private:
+ vector<double> fFatima_Energy;
+ vector<double> fFatima_Time;
+ vector<short> fFatima_Number;
+ */
+
+ public:
+ TFatimaData();
+ virtual ~TFatimaData();
+
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+ ///////////////////// GETTERS ////////////////////////
+
+ //
+ // First stage
+ //
+ // (E)
+ UShort_t GetFATIMALaBr3StageEMult() {
+ return fFATIMA_LaBr3Stage_E_DetectorNbr.size(); // TODO: Maybe change to CrystalNbr for multiplicity
+ }
+ UShort_t GetFATIMALaBr3StageEDetectorNbr(Int_t i) {
+ return fFATIMA_LaBr3Stage_E_DetectorNbr.at(i);
+ }
+ UShort_t GetFATIMALaBr3StageECrystalNbr(Int_t i) {
+ return fFATIMA_LaBr3Stage_E_CrystalNbr.at(i);
+ }
+ Double_t GetFATIMALaBr3StageEEnergy(Int_t i) {
+ return fFATIMA_LaBr3Stage_E_Energy.at(i);
+ }
+ /*
+ Double_t GetFATIMALaBr3StageEffphpeak(Int_t i) {
+ return fFATIMA_LaBr3Stage_Eff_phpeak.at(i);
+ }
+ */
+
+ // (T)
+ UShort_t GetFATIMALaBr3StageTMult() {
+ return fFATIMA_LaBr3Stage_E_DetectorNbr.size();
+ }
+ UShort_t GetFATIMALaBr3StageTDetectorNbr(Int_t i) {
+ return fFATIMA_LaBr3Stage_T_DetectorNbr.at(i);
+ }
+ UShort_t GetFATIMALaBr3StageTCrystalNbr(Int_t i) {
+ return fFATIMA_LaBr3Stage_T_CrystalNbr.at(i);
+ }
+ Double_t GetFATIMALaBr3StageTTime(Int_t i) {
+ return fFATIMA_LaBr3Stage_T_Time.at(i);
+ }
+
+ /*
+ //
+ // Second stage (CsI
+ //
+ // (E)
+ UShort_t GetFATIMACsIStageEMult() {
+ return fFATIMA_CsIStage_E_DetectorNbr.size();
+ }
+ UShort_t GetFATIMACsIStageEDetectorNbr(Int_t i) {
+ return fFATIMA_CsIStage_E_DetectorNbr.at(i);
+ }
+ UShort_t GetFATIMACsIStageECrystalNbr(Int_t i) {
+ return fFATIMA_CsIStage_E_CrystalNbr.at(i);
+ }
+ Double_t GetFATIMACsIStageEEnergy(Int_t i) {
+ return fFATIMA_CsIStage_E_Energy.at(i);
+ }
+ */
+ /*
+ Double_t GetFATIMACsIStageEffphpeak(Int_t i) {
+ return fFATIMA_CsIStage_Eff_phpeak.at(i);
+ }
+ */
+ /*
+ // (T)
+ UShort_t GetFATIMACsIStageTMult() {
+ return fFATIMA_CsIStage_E_DetectorNbr.size();
+ }
+ UShort_t GetFATIMACsIStageTDetectorNbr(Int_t i) {
+ return fFATIMA_CsIStage_T_DetectorNbr.at(i);
+ }
+ UShort_t GetFATIMACsIStageTCrystalNbr(Int_t i) {
+ return fFATIMA_CsIStage_T_CrystalNbr.at(i);
+ }
+ Double_t GetFATIMACsIStageTTime(Int_t i) {
+ return fFATIMA_CsIStage_T_Time.at(i);
+ }
+ */
+
+ /*
+ // (E)
+ //double GetEnergy(int i) {return fFatima_Energy[i];}
+ // (T)
+ //double GetTime(int i) {return fFatima_Time[i];}
+ // (N)
+ int GetFatimaNumber(int i) {return fFatima_Number[i];}
+ double GetEnergySize() {return fFatima_Energy.size();}
+ // (T)
+ double GetTimeSize() {return fFatima_Time.size();}
+ // (N)
+ int GetFatimaNumberSize() {return fFatima_Number.size();}
+ */
+
+ ///////////////////// SETTERS ////////////////////////
+
+ //
+ // First stage
+ //
+ // (E)
+
+ void SetFATIMALaBr3StageEDetectorNbr(UShort_t DetNbr) {
+ fFATIMA_LaBr3Stage_E_DetectorNbr.push_back(DetNbr);
+ }
+ void SetFATIMALaBr3StageECrystalNbr(UShort_t CrystalNbr) {
+ fFATIMA_LaBr3Stage_E_CrystalNbr.push_back(CrystalNbr);
+ }
+ void SetFATIMALaBr3StageEEnergy(Double_t Energy) {
+ fFATIMA_LaBr3Stage_E_Energy.push_back(Energy);
+ }
+ /*
+ void SetFATIMALaBr3StageEffphpeak(Double_t Energy) {
+ fFATIMA_LaBr3Stage_Eff_phpeak.push_back(Energy);
+ }
+ */
+
+ // (T)
+ void SetFATIMALaBr3StageTDetectorNbr(UShort_t DetNbr) {
+ fFATIMA_LaBr3Stage_T_DetectorNbr.push_back(DetNbr);
+ }
+ void SetFATIMALaBr3StageTCrystalNbr(UShort_t CrystalNbr) {
+ fFATIMA_LaBr3Stage_T_CrystalNbr.push_back(CrystalNbr);
+ }
+ void SetFATIMALaBr3StageTTime(Double_t Time) {
+ fFATIMA_LaBr3Stage_T_Time.push_back(Time);
+ }
+
+ /*
+ //
+ // Second stage (CsI
+ //
+ // (E)
+ void SetFATIMACsIStageEDetectorNbr(UShort_t DetNbr) {
+ fFATIMA_CsIStage_E_DetectorNbr.push_back(DetNbr);
+ }
+ void SetFATIMACsIStageECrystalNbr(UShort_t CrystalNbr) {
+ fFATIMA_CsIStage_E_CrystalNbr.push_back(CrystalNbr);
+ }
+ void SetFATIMACsIStageEEnergy(Double_t Energy) {
+ fFATIMA_CsIStage_E_Energy.push_back(Energy);
+ }
+ */
+
+ /*
+ void SetFATIMACsIStageEffphpeak(Double_t Energy) {
+ fFATIMA_CsIStage_Eff_phpeak.push_back(Energy);
+ }
+ */
+
+ /*
+ // (T)
+ void SetFATIMACsIStageTDetectorNbr(UShort_t DetNbr) {
+ fFATIMA_CsIStage_T_DetectorNbr.push_back(DetNbr);
+ }
+ void SetFATIMACsIStageTCrystalNbr(UShort_t CrystalNbr) {
+ fFATIMA_CsIStage_T_CrystalNbr.push_back(CrystalNbr);
+ }
+ void SetFATIMACsIStageTTime(Double_t Time) {
+ fFATIMA_CsIStage_T_Time.push_back(Time);
+ }
+ */
+
+ /*
+ // (E)
+ void SetEnergy(double E) {fFatima_Energy.push_back(E);}
+ void SetTime(double T) {fFatima_Time.push_back(T);}
+ void SetFatimaNumber(int N) {fFatima_Number.push_back(N);}
+ */
+ ClassDef(TFatimaData,1) // FatimaData structure
+};
+
+#endif
diff --git a/NPLib/Fatima/TFatimaPhysics.cxx b/NPLib/Fatima/TFatimaPhysics.cxx
new file mode 100755
index 000000000..199d02201
--- /dev/null
+++ b/NPLib/Fatima/TFatimaPhysics.cxx
@@ -0,0 +1,356 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/2014 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class stores the physical results after NPAnalysis is run*
+ * for the FATIMA detector. *
+ * This class derives from TObject (ROOT) and its aim is to be *
+ * stored in the output TTree of NPAnalysis. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "TFatimaPhysics.h"
+#include <vector>
+#include <iostream>
+#include <cstdlib>
+
+
+ClassImp(TFatimaPhysics)
+
+
+TFatimaPhysics::TFatimaPhysics()
+{
+ //EventMultiplicity = 0;
+}
+
+
+
+TFatimaPhysics::~TFatimaPhysics()
+{
+ Clear();
+}
+
+
+
+void TFatimaPhysics::BuildSimplePhysicalEvent(TFatimaData* Data)
+{
+ BuildPhysicalEvent(Data);
+}
+
+
+
+void TFatimaPhysics::BuildPhysicalEvent(TFatimaData* Data)
+{
+
+ int multLaBrE = Data->GetFATIMALaBr3StageEMult();
+ // int multCsIE = Data->GetFATIMACsIStageEMult();
+
+ cout << "%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+ cout << "multLaBr= " << multLaBrE << endl;
+ //cout << "multCsI= " << multCsIE << endl;
+
+ FatimaEventMult=multLaBrE;
+ //FatimaEventMult=multLaBrE+multCsIE;
+
+ // Get energy from strips and store it
+
+ if(FatimaEventMult>=1)
+ {
+
+ double EnergyTot=0.;
+
+ double HighestEnergy;
+ int DetID_HighestE;
+ double CrystID_HighestE;
+ // Initialisation
+ HighestEnergy=0;
+ DetID_HighestE=-1;
+ CrystID_HighestE=-1;
+
+ /*
+ double HighestEnergy, Outer_HighestEnergy;
+ int DetID_HighestE, OuterDetID_HighestE;
+ double CrystID_HighestE, OuterCrystID_HighestE;
+
+ // Initialisation
+ HighestEnergy=Outer_HighestEnergy=0;
+ DetID_HighestE=OuterDetID_HighestE=-1;
+ CrystID_HighestE=OuterCrystID_HighestE=-1;
+ */
+
+ if(multLaBrE>=1){
+ //cout << "cava1b" <<endl;
+ //cout << Data->GetFATIMALaBr3StageEEnergy(0) <<endl;
+ //cout << "cava1b" <<endl;
+
+ double EnergyTotFirst;
+ double EnergyStripFront;
+ double EnergyStrip;
+ //double HighestEnergyStrip;
+
+ for(int j=0;j<multLaBrE;j++)
+ {
+ EnergyStripFront= Data->GetFATIMALaBr3StageEEnergy(j);
+
+ EnergyStrip = EnergyStripFront;
+ cout << "Energy LaBr=" << EnergyStrip << endl;
+ FatimaLaBr_E.push_back(EnergyStrip);
+
+ EnergyTotFirst += EnergyStrip;
+ EnergyTot += EnergyStrip;
+
+ //Let's find the first detector hit:
+ if(HighestEnergy<EnergyStripFront)
+ {
+ HighestEnergy=EnergyStripFront;
+ DetID_HighestE=Data->GetFATIMALaBr3StageEDetectorNbr(j);
+ CrystID_HighestE=Data->GetFATIMALaBr3StageECrystalNbr(j);
+ }
+
+
+ }
+
+ // Fill total energy in inner shell and the first detector/crystal hit
+
+ FatimaInTotalEnergy.push_back(EnergyTotFirst);
+
+ cout << " Id of Fatima detector hit= " << DetID_HighestE << endl;
+ cout << " EnergyTotFirst= " << EnergyTotFirst << endl;
+ cout << " EnergyTot= " << EnergyTot << endl;
+
+ if(DetID_HighestE>100)// Cluster case
+ {
+ //ModuleNumber.push_back(DetID_HighestE);
+ DetectorNumber.push_back(-1);
+ }else // Detector case
+ {
+ //ModuleNumber.push_back(-1);
+ DetectorNumber.push_back(DetID_HighestE);
+ }
+
+ if(DetID_HighestE<100) // Detector case
+ {
+ FirstStage_X.push_back(0);
+ FirstStage_Y.push_back(0);
+ }
+
+
+ /*
+ if(DetID_HighestE>100) // Cluster case
+ {
+
+ if(CrystID_HighestE==0)
+ {
+ FirstStage_X.push_back(0);
+ FirstStage_Y.push_back(0);
+ }
+ if(CrystID_HighestE==1)
+ {
+ FirstStage_X.push_back(1);
+ FirstStage_Y.push_back(0);
+ }
+ if(CrystID_HighestE==2)
+ {
+ FirstStage_X.push_back(2);
+ FirstStage_Y.push_back(0);
+ }
+ if(CrystID_HighestE==3)
+ {
+ FirstStage_X.push_back(0);
+ FirstStage_Y.push_back(1);
+ }
+ if(CrystID_HighestE==4)
+ {
+ FirstStage_X.push_back(1);
+ FirstStage_Y.push_back(1);
+ }
+ if(CrystID_HighestE==5)
+ {
+ FirstStage_X.push_back(2);
+ FirstStage_Y.push_back(1);
+ }
+ if(CrystID_HighestE==6)
+ {
+ FirstStage_X.push_back(0);
+ FirstStage_Y.push_back(2);
+ }
+ if(CrystID_HighestE==7)
+ {
+ FirstStage_X.push_back(1);
+ FirstStage_Y.push_back(2);
+ }
+ if(CrystID_HighestE==8)
+ {
+ FirstStage_X.push_back(2);
+ FirstStage_Y.push_back(2);
+ }
+ } // end of cluster case
+ */
+
+ }
+
+ /*
+ if(multCsIE>=1){
+ double EnergySecond;
+ double EnergyTotSecond;
+ for(int j=0;j<multCsIE;j++)
+ {
+ EnergySecond = Data->GetFATIMACsIStageEEnergy(j);
+ FatimaCsI_E.push_back(EnergySecond);
+ EnergyTotSecond +=EnergySecond;
+
+ EnergyTot += EnergySecond;
+ cout << "Energy CsI=" << EnergySecond << endl;
+ //cout << "Energytot CsI=" << EnergyTot << endl;
+
+ //Let's find the first outer detector hit:
+ if(Outer_HighestEnergy<EnergySecond)
+ {
+ Outer_HighestEnergy=EnergySecond;
+ OuterDetID_HighestE=Data->GetFATIMACsIStageEDetectorNbr(j);
+ OuterCrystID_HighestE=Data->GetFATIMACsIStageECrystalNbr(j);
+ }
+ }
+
+ // Fill total energy in outer shell
+ FatimaOutTotalEnergy.push_back(EnergyTotSecond);
+
+ cout << " Id of Fatima outer detector hit= " << OuterDetID_HighestE << endl;
+ cout << " EnergyTotSecond= " << EnergyTotSecond << endl;
+ cout << " EnergyTot= " << EnergyTot << endl;
+
+ if(multLaBrE==0 || (EnergyTotSecond>EnergyTot) ) { // = only if Outer layer is hit first
+
+ FirstStage_X.push_back(-1);
+ FirstStage_Y.push_back(-1);
+
+ cout << " Id of Fatima (CsI) detector hit= " << OuterDetID_HighestE << endl;
+
+ if(OuterDetID_HighestE>100)
+ {
+ ModuleNumber.push_back(OuterDetID_HighestE);
+ DetectorNumber.push_back(-1);
+ }else
+ {
+ ModuleNumber.push_back(-1);
+ DetectorNumber.push_back(OuterDetID_HighestE);
+ }
+
+
+ if(OuterDetID_HighestE<100) // Detector case
+ {
+ SecondStage_X.push_back(0);
+ SecondStage_Y.push_back(0);
+ }
+
+
+ if(OuterDetID_HighestE>100) // Cluster case
+ {
+ if(OuterCrystID_HighestE==0)
+ {
+ SecondStage_X.push_back(0);
+ SecondStage_Y.push_back(0);
+ }
+ if(OuterCrystID_HighestE==1)
+ {
+ SecondStage_X.push_back(1);
+ SecondStage_Y.push_back(0);
+ }
+ if(OuterCrystID_HighestE==2)
+ {
+ SecondStage_X.push_back(2);
+ SecondStage_Y.push_back(0);
+ }
+ if(OuterCrystID_HighestE==3)
+ {
+ SecondStage_X.push_back(0);
+ SecondStage_Y.push_back(1);
+ }
+ if(OuterCrystID_HighestE==4)
+ {
+ SecondStage_X.push_back(1);
+ SecondStage_Y.push_back(1);
+ }
+ if(OuterCrystID_HighestE==5)
+ {
+ SecondStage_X.push_back(2);
+ SecondStage_Y.push_back(1);
+ }
+ if(OuterCrystID_HighestE==6)
+ {
+ SecondStage_X.push_back(0);
+ SecondStage_Y.push_back(2);
+ }
+ if(OuterCrystID_HighestE==7)
+ {
+ SecondStage_X.push_back(1);
+ SecondStage_Y.push_back(2);
+ }
+ if(OuterCrystID_HighestE==8)
+ {
+ SecondStage_X.push_back(2);
+ SecondStage_Y.push_back(2);
+ }
+
+ }
+ }
+
+ }
+ */
+
+ // Fill total energy
+ FatimaTotalEnergy.push_back(EnergyTot);
+ }
+}
+
+
+
+void TFatimaPhysics::Clear()
+{
+ //EventMultiplicity= 0;
+ FatimaEventMult= 0;
+ DetectorNumber.clear(); // phoswich #
+
+ //EventType.clear();
+ FatimaInTotalEnergy.clear(); // inner shell
+ //FatimaOutTotalEnergy.clear(); // outter shell
+ FatimaTotalEnergy.clear();
+
+
+ // LaBr
+ FatimaLaBr_E.clear();
+ //FatimaIn_1stDetId.clear();
+ //FatimaIn_1stCrystId.clear();
+ //FirstStage_T.clear();
+ FirstStage_X.clear();
+ FirstStage_Y.clear();
+
+ /* // CsI
+ FatimaCsI_E.clear();
+ //FatimaOut_1stDetId.clear();
+ //FatimaOut_1stCrystId.clear();
+ //SecondStage_T.clear();
+ //SecondStage_N.clear();
+ SecondStage_X.clear();
+ SecondStage_Y.clear();
+ */
+
+ /*
+ // CsI
+ ThirdStage_E.clear();
+ ThirdStage_T.clear();
+ ThirdStage_N.clear();
+ */
+}
diff --git a/NPLib/Fatima/TFatimaPhysics.h b/NPLib/Fatima/TFatimaPhysics.h
new file mode 100755
index 000000000..4a4592765
--- /dev/null
+++ b/NPLib/Fatima/TFatimaPhysics.h
@@ -0,0 +1,91 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/2014 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class stores the physical results after NPAnalysis is run*
+ * for the FATIMA detector. *
+ * This class derives from TObject (ROOT) and its aim is to be *
+ * stored in the output TTree of NPAnalysis. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#ifndef TFATIMAPHYSICS_H
+#define TFATIMAPHYSICS_H
+
+#include <vector>
+#include "TObject.h"
+#include "TFatimaData.h"
+#include <cstdlib>
+
+using namespace std ;
+
+class TFatimaPhysics : public TObject
+{
+public:
+ TFatimaPhysics();
+ ~TFatimaPhysics();
+
+public:
+ void Clear();
+ void Clear(const Option_t*) {};
+ void BuildPhysicalEvent(TFatimaData* Data);
+ void BuildSimplePhysicalEvent(TFatimaData* Data);
+
+public:
+ // Provide Physical Multiplicity
+ Int_t FatimaEventMult;
+
+ // Provide a Classification of Event
+ //vector<int> EventType;
+
+
+ // Telescope
+ //vector<int> ModuleNumber; // Id of cluster
+
+ vector<int> DetectorNumber; // Id of Detector
+
+ // FirstStage
+ vector<double> FatimaLaBr_E;
+ //vector<double> FirstStage_T;
+ vector<int> FirstStage_X; // column # of LaBr crystal in 3x3 cluster
+ vector<int> FirstStage_Y; // raw # of LaBr crystal in 3x3 cluster
+
+ /*
+ // SecondStage
+ vector<double> FatimaCsI_E;
+ //vector<double> SecondStage_T;
+ //vector<int> SecondStage_N;
+ vector<int> SecondStage_X; // column # of CsI crystal in 3x3 cluster
+ vector<int> SecondStage_Y; // raw # of CsI crystal in 3x3 cluster
+ */
+
+
+ /*
+ // ThirdStage
+ vector<double> ThirdStage_E;
+ vector<double> ThirdStage_T;
+ vector<int> ThirdStage_N;
+ */
+
+ // Physical Value
+ vector<double> FatimaTotalEnergy;
+ vector<double> FatimaInTotalEnergy;
+ vector<double> FatimaOutTotalEnergy;
+
+
+ ClassDef(TFatimaPhysics,1) // GaspardTrackerPHysics structure
+};
+
+#endif
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index 030b87339..f6c99e03e 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -1,3 +1,4 @@
+
/*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
@@ -32,6 +33,7 @@
#include "Hyde2Tracker.h"
#include "Paris.h"
#include "Shield.h"
+#include "Fatima.h"
#include "TAnnularS1Physics.h"
#include "TCATSPhysics.h"
#include "TCharissaPhysics.h"
@@ -93,6 +95,7 @@ void DetectorManager::ReadConfigurationFile(string Path) {
bool IonisationChamber = false;
bool LaBr3 = false;
bool MUST2 = false;
+ bool FatimaDet = false;
bool ParisDet = false;
bool S2 = false;
bool SPEG = false;
@@ -252,8 +255,25 @@ void DetectorManager::ReadConfigurationFile(string Path) {
AddDetector("EXOGAM", myDetector);
#endif
}
-
-
+ ////////////////////////////////////////////
+ ///////////// Search for Fatima ////////////
+ ////////////////////////////////////////////
+ else if (LineBuffer.compare(0, 6, "Fatima") == 0 && FatimaDet == false) {
+#ifdef INC_FATIMA
+ FatimaDet = true;
+ cout << "//////// FATIMA ////////" << endl << endl;
+
+ // Instantiate the new array as a VDetector Object
+ VDetector* myDetector = new Fatima();
+ // Read Position of Telescope
+ ConfigFile.close();
+ myDetector->ReadConfiguration(Path);
+ ConfigFile.open(Path.c_str());
+
+ // Add array to the VDetector Vector
+ AddDetector("FATIMA", myDetector);
+#endif
+ }
////////////////////////////////////////////
//////////// Search for Gaspard ////////////
////////////////////////////////////////////
diff --git a/NPSimulation/Fatima/Fatima.cc b/NPSimulation/Fatima/Fatima.cc
new file mode 100755
index 000000000..d7951824e
--- /dev/null
+++ b/NPSimulation/Fatima/Fatima.cc
@@ -0,0 +1,133 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: This class manages different shapes of module for the Fatima *
+ * detector. It allows to have Fatima geometries with an *
+ * heterogeneous set of modules *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <cmath>
+
+// NPTool headers
+#include "Fatima.hh"
+#include "FatimaDetector.hh"
+
+using namespace std;
+
+
+Fatima::Fatima()
+{
+}
+
+
+
+Fatima::~Fatima()
+{
+}
+
+
+
+// Read stream at Configfile to pick-up parameters of detector (Position,...)
+// Called in DetecorConstruction::ReadDetextorConfiguration Method
+void Fatima::ReadConfiguration(string Path)
+{
+ // open configuration file
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+
+ bool bFatimaDetector = false;
+
+ string LineBuffer;
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+ if (LineBuffer.compare(0, 14, "FatimaDetector") == 0 && bFatimaDetector == false) {
+ bFatimaDetector = true;
+
+ // instantiate a new "detector" corresponding to the Square elements
+ FatimaModule* myDetector = new FatimaDetector();
+
+ // read part of the configuration file corresponding to square elements
+ ConfigFile.close();
+ myDetector->ReadConfiguration(Path);
+ ConfigFile.open(Path.c_str());
+
+ // ms_InterCoord comes from VDetector
+ myDetector->SetInterCoordPointer(ms_InterCoord);
+
+ // store FatimaDetector "detector"
+ m_Modules.push_back(myDetector);
+ }
+
+ }
+}
+
+
+
+// Construct detector and initialize sensitive part.
+// Called After DetectorConstruction::AddDetector Method
+void Fatima::ConstructDetector(G4LogicalVolume* world)
+{
+ // loop on sub-detectors belonging to Fatima
+ int nbDetectors = m_Modules.size();
+ for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ConstructDetector(world);
+}
+
+
+
+// Connect the FatimaData class to the output TTree
+// of the simulation
+void Fatima::InitializeRootOutput()
+{
+ // loop on sub-detectors belonging to Fatima
+ int nbDetectors = m_Modules.size();
+ for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeRootOutput();
+}
+
+
+
+// Initialize all scorers necessary for each detector
+void Fatima::InitializeScorers()
+{
+ // loop on sub-detectors belonging to Fatima
+ int nbDetectors = m_Modules.size();
+ for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeScorers();
+}
+
+
+
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAction
+void Fatima::ReadSensitive(const G4Event* event)
+{
+ // Before looping on each sub-detector, clear the static variable
+ // ms_InterCoord
+ // This is done on the first element of the m_Modules vector.
+ // This should be done here since this variable (of type TIneractionCoordinates)
+ // deals with multiplicity of events > 1.
+ m_Modules[0]->GetInterCoordPointer()->Clear();
+
+ // We do the same for the static variable ms_Event
+ m_Modules[0]->GetEventPointer()->Clear();
+
+ // loop on sub-detectors belonging to Fatima
+ int nbDetectors = m_Modules.size();
+ for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ReadSensitive(event);
+}
diff --git a/NPSimulation/Fatima/Fatima.hh b/NPSimulation/Fatima/Fatima.hh
new file mode 100755
index 000000000..545369180
--- /dev/null
+++ b/NPSimulation/Fatima/Fatima.hh
@@ -0,0 +1,73 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: This class manages different shapes of module for the Fatima *
+ * detector. It allows to have Fatima geometries with an *
+ * heterogeneous set of modules *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#ifndef Fatima_h
+#define Fatima_h 1
+
+// C++ headers
+#include <vector>
+
+// NPTool header
+#include "VDetector.hh"
+#include "FatimaModule.hh"
+
+using namespace std;
+
+
+
+class Fatima : public VDetector
+{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+public:
+ Fatima();
+ virtual ~Fatima();
+
+ ////////////////////////////////////////////////////
+ ///////// Inherite from VDetector class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(string Path);
+
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world);
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput();
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event);
+
+public:
+ // Initialize all scorers necessary for each detector
+ void InitializeScorers();
+
+private:
+ vector<FatimaModule*> m_Modules;
+};
+#endif
diff --git a/NPSimulation/Fatima/FatimaDetector.cc b/NPSimulation/Fatima/FatimaDetector.cc
new file mode 100755
index 000000000..98bc7d6a8
--- /dev/null
+++ b/NPSimulation/Fatima/FatimaDetector.cc
@@ -0,0 +1,1237 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/14 *
+ *---------------------------------------------------------------------------*
+ * Decription: Define a detector module for the Fatima detector. *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+
+// C++ headers
+#include <sstream>
+#include <string>
+#include <cmath>
+
+// G4 Geometry headers
+#include "G4Trd.hh"
+#include "G4Box.hh"
+#include "G4Tubs.hh"
+#include "G4Trap.hh"
+
+#include "G4SubtractionSolid.hh"
+
+// G4 various headers
+#include "G4MaterialTable.hh"
+#include "G4Element.hh"
+#include "G4ElementTable.hh"
+#include "G4VisAttributes.hh"
+#include "G4Colour.hh"
+#include "G4RotationMatrix.hh"
+#include "G4Transform3D.hh"
+#include "G4PVPlacement.hh"
+#include "G4PVDivision.hh"
+
+// G4 sensitive
+#include "G4SDManager.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+// NPTool headers
+#include "FatimaDetector.hh"
+#include "ObsoleteGeneralScorers.hh"
+#include "FatimaScorers.hh"
+#include "RootOutput.h"
+
+// CLHEP
+#include "CLHEP/Random/RandGauss.h"
+
+using namespace std;
+using namespace CLHEP;
+using namespace FATIMADETECTOR;
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+FatimaDetector::FatimaDetector()
+{
+ ms_InterCoord = 0;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+FatimaDetector::~FatimaDetector()
+{
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void FatimaDetector::AddModule(G4ThreeVector X1_Y1,
+ G4ThreeVector X128_Y1,
+ G4ThreeVector X1_Y128,
+ G4ThreeVector X128_Y128)
+{
+ m_DefinitionType.push_back(true);
+
+ m_X1_Y1.push_back(X1_Y1);
+ m_X128_Y1.push_back(X128_Y1);
+ m_X1_Y128.push_back(X1_Y128);
+ m_X128_Y128.push_back(X128_Y128);
+
+ m_R.push_back(0);
+ m_Theta.push_back(0);
+ m_Phi.push_back(0);
+ m_beta_u.push_back(0);
+ m_beta_v.push_back(0);
+ m_beta_w.push_back(0);
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void FatimaDetector::AddModule(G4double R,
+ G4double Theta,
+ G4double Phi,
+ G4double beta_u,
+ G4double beta_v,
+ G4double beta_w)
+{
+ G4ThreeVector empty = G4ThreeVector(0, 0, 0);
+
+ m_DefinitionType.push_back(false);
+
+ m_X1_Y1.push_back(empty);
+ m_X128_Y1.push_back(empty);
+ m_X1_Y128.push_back(empty);
+ m_X128_Y128.push_back(empty);
+
+ m_R.push_back(R);
+ m_Theta.push_back(Theta);
+ m_Phi.push_back(Phi);
+ m_beta_u.push_back(beta_u);
+ m_beta_v.push_back(beta_v);
+ m_beta_w.push_back(beta_w);
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void FatimaDetector::VolumeMaker(G4int DetecNumber,
+ G4ThreeVector MMpos,
+ G4RotationMatrix* MMrot,
+ G4LogicalVolume* world)
+{
+ G4double NbrTelescopes = DetecNumber;
+ G4String DetectorNumber;
+ ostringstream Number;
+ Number << NbrTelescopes;
+ DetectorNumber = Number.str();
+
+ /////////////////////////////////////////////////////////////////
+ /////////////////Element Definition ///////////////////////////
+ ////////////////////////////////////////////////////////////////
+ G4String symbol, name;
+ G4double density = 0. , a = 0, z = 0;
+ G4int ncomponents = 0;
+ G4int nel = 0, natoms = 0;
+
+ ////////////////////////////////////////////////////////////////
+ /////////////////Material Definition ///////////////////////////
+ ////////////////////////////////////////////////////////////////
+ a=137.327*g/mole;
+ G4Element* Ba = new G4Element(name="Barium",symbol="Ba",z=56.,a);
+ a=18.9984032*g/mole;
+ G4Element* F = new G4Element(name="Fluorine",symbol="F",z=9.,a);
+ a=22.99*g/mole;
+ G4Element* Na = new G4Element(name="Sodium",symbol="Na",z=11.,a);
+ a=79.904*g/mole;
+ G4Element* Br = new G4Element(name="Bromine",symbol="Br",z=35.,a);
+ a=126.90477*g/mole;
+ G4Element* I = new G4Element(name="Iodine",symbol="I",z=53.,a);
+ a=132.90545*g/mole;
+ G4Element* Cs = new G4Element(name="Cesium",symbol="Cs",z=55.,a);
+ a=138.9055*g/mole;
+ G4Element* La = new G4Element(name="Lanthanum",symbol="La",z=57.,a);
+
+ // Vacuum
+ G4Element* N = new G4Element("Nitrogen" , symbol = "N" , z = 7 , a = 14.01 * g / mole);
+ G4Element* O = new G4Element("Oxigen" , symbol = "O" , z = 8 , a = 16.00 * g / mole);
+
+ density = 0.000000001 * mg / cm3;
+ G4Material* Vacuum = new G4Material("Vacuum", density, ncomponents = 2);
+ Vacuum->AddElement(N, .7);
+ Vacuum->AddElement(O, .3);
+
+ // Germanium
+ density = 5.32*g/cm3;
+ G4Material* Ge;
+ Ge = new G4Material(name="Germanium",32,72.59*g/mole, density);
+
+ // Aluminum
+ density = 2.7*g/cm3;
+ G4Material* Alu;
+ Alu = new G4Material(name="Aluminum",13,26.98154*g/mole, density);
+
+ // lead
+ density =11.35*g/cm3;
+ G4Material* Lead;
+ Lead = new G4Material(name="Lead",82,207.2*g/mole, density);
+
+ // NaI
+ density = 3.67*g/cm3, nel = 2;
+ G4Material* NaI = new G4Material(name="NaI",density,nel);
+ NaI->AddElement(Na, natoms = 1);
+ NaI->AddElement(I, natoms = 1);
+
+ // CsI
+ density = 4.51*g/cm3, nel = 2;
+ G4Material* CsI = new G4Material(name="CsI", density, nel);
+ CsI->AddElement(Cs, natoms = 1);
+ CsI->AddElement(I, natoms = 1);
+
+ // LaBr3
+ density = 5.29*g/cm3, nel = 2;
+ G4Material* LaBr3 = new G4Material(name="LaBr3",density,nel);
+ LaBr3->AddElement(La, natoms = 1);
+ LaBr3->AddElement(Br, natoms = 3);
+
+ // BaF2
+ density = 4.89*g/cm3, nel = 2;
+ G4Material* BaF2 = new G4Material(name="BaF2", density, nel);
+ BaF2->AddElement(Ba, natoms = 1);
+ BaF2->AddElement(F, natoms = 2);
+
+ ////////////////////////////////////////////////////////////////
+ ////////////// Starting Volume Definition //////////////////////
+ ////////////////////////////////////////////////////////////////
+ // Little trick to avoid warning in compilation: Use a PVPlacement "buffer".
+ // If don't you will have a Warning unused variable 'myPVP'
+ G4PVPlacement* PVPBuffer;
+ G4String Name = "FatimaDetector" + DetectorNumber ;
+
+ // Mother Volume
+ //G4Box* solidFatimaDetector = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length);
+ G4Tubs* solidFatimaDetector = new G4Tubs(Name,0, 0.5*FaceFront, 0.5*Length, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicFatimaDetector = new G4LogicalVolume(solidFatimaDetector, Vacuum, Name, 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot, MMpos) ,
+ logicFatimaDetector ,
+ Name ,
+ world ,
+ false ,
+ 0);
+
+ logicFatimaDetector->SetVisAttributes(G4VisAttributes::Invisible);
+ //logicFatimaDetector->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
+ //if (m_non_sensitive_part_visiualisation) logicFatimaDetector->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
+
+ // Detector construction
+ // LaBr3
+ G4ThreeVector positionLaBr3Stage = G4ThreeVector(0, 0, LaBr3Stage_PosZ);
+
+ //G4Box* solidLaBr3Stage = new G4Box("solidLaBr3Stage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*LaBr3Thickness);
+ G4Tubs* solidLaBr3Stage = new G4Tubs("solidLaBr3Stage", 0., 0.5*LaBr3Face, 0.5*LaBr3Thickness, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, LaBr3, "logicLaBr3Stage", 0, 0, 0);
+ //G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, Ge, "logicLaBr3Stage", 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(0,
+ positionLaBr3Stage,
+ logicLaBr3Stage,
+ Name + "_LaBr3Stage",
+ logicFatimaDetector,
+ false,
+ 0);
+
+ // Set LaBr3 sensible
+ logicLaBr3Stage->SetSensitiveDetector(m_LaBr3StageScorer);
+
+ // Visualisation of LaBr3Stage Strip
+ G4VisAttributes* LaBr3VisAtt = new G4VisAttributes(G4Colour(1., 0., 0.));
+ logicLaBr3Stage->SetVisAttributes(LaBr3VisAtt);
+
+
+ // Aluminium can around LaBr3
+ // LaBr3 Can
+ G4ThreeVector positionLaBr3Can = G4ThreeVector(0, 0, LaBr3Can_PosZ);
+
+ G4Tubs* solidLaBr3Can = new G4Tubs("solidLaBr3Can", 0.5*CanInnerDiameter, 0.5*CanOuterDiameter, 0.5*CanLength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLaBr3Can = new G4LogicalVolume(solidLaBr3Can, Alu, "logicLaBr3Can", 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(0,
+ positionLaBr3Can,
+ logicLaBr3Can,
+ Name + "_LaBr3Can",
+ logicFatimaDetector,
+ false,
+ 0);
+
+
+ // Visualisation of LaBr3Can
+ G4VisAttributes* CanVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
+ logicLaBr3Can->SetVisAttributes(CanVisAtt);
+
+ // Aluminium window in front of LaBr3
+ // LaBr3 Window
+ G4ThreeVector positionLaBr3Win = G4ThreeVector(0, 0, LaBr3Win_PosZ);
+
+ G4Tubs* solidLaBr3Win = new G4Tubs("solidLaBr3Win", 0.5*WinInnerDiameter, 0.5*WinOuterDiameter, 0.5*WinLength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLaBr3Win = new G4LogicalVolume(solidLaBr3Win, Alu, "logicLaBr3Win", 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(0,
+ positionLaBr3Win,
+ logicLaBr3Win,
+ Name + "_LaBr3Win",
+ logicFatimaDetector,
+ false,
+ 0);
+ // Visualisation of LaBr3Win
+ logicLaBr3Win->SetVisAttributes(CanVisAtt);
+
+
+ // CsI stage replaced by Aluminum housing (TODO: replace CsI by Alu in variable names to avoid confusion)
+ G4ThreeVector positionCsIStage = G4ThreeVector(0, 0, CsIStage_PosZ);
+
+ //G4Box* solidCsIStage = new G4Box("solidCsIStage", 0.5*CsIFace, 0.5*CsIFace, 0.5*CsIThickness);
+
+ G4Tubs* solidPMout = new G4Tubs("solidPMOut", 0.5*LaBr3Face, 0.5*CsIFace, 0.5*CsIThickness, 0.*deg, 360.*deg);
+ G4Tubs* solidPMin = new G4Tubs("solidPMIn", 0.5*LaBr3Face-0.1*cm, 0.5*CsIFace-0.5*cm, 0.5*(CsIThickness-2.*cm)-0.1*cm, 0.*deg, 360.*deg);
+ G4RotationMatrix* RotMat=NULL;
+ const G4ThreeVector &Trans= G4ThreeVector(0.,0.,1.*cm);
+ G4SubtractionSolid* solidCsIStage = new G4SubtractionSolid("solidCsIStage", solidPMout,solidPMin, RotMat, Trans);
+
+ //G4Tubs* solidCsIStage = new G4Tubs("solidCsIStage", 0.5*LaBr3Face, 0.5*CsIFace, 0.5*CsIThickness, 0.*deg, 360.*deg);
+ //G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, CsI, "logicCsIStage", 0, 0, 0);
+ G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, Alu, "logicCsIStage", 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(0,
+ positionCsIStage,
+ logicCsIStage,
+ Name + "_CsIStage",
+ logicFatimaDetector,
+ false,
+ 0);
+
+
+ // Visualisation of CsIStage Strip
+ G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
+ logicCsIStage->SetVisAttributes(CsIVisAtt);
+
+
+
+ // Lead shielding
+ // A
+ G4ThreeVector positionLeadAShield = G4ThreeVector(0, 0, LeadAShield_PosZ);
+ G4Tubs* solidLeadA = new G4Tubs("solidLead", 0.5*LeadAMinR, 0.5*LeadAMaxR, 0.5*LeadALength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLeadAShield = new G4LogicalVolume(solidLeadA, Lead, "logicLeadAShield", 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(0,
+ positionLeadAShield,
+ logicLeadAShield,
+ Name + "_LeadAShield",
+ logicFatimaDetector,
+ false,
+ 0);
+ // B
+ G4ThreeVector positionLeadBShield = G4ThreeVector(0, 0, LeadBShield_PosZ);
+ G4Tubs* solidLeadB = new G4Tubs("solidLead", 0.5*LeadBMinR, 0.5*LeadBMaxR, 0.5*LeadBLength, 0.*deg, 360.*deg);
+ G4LogicalVolume* logicLeadBShield = new G4LogicalVolume(solidLeadB, Lead, "logicLeadBShield", 0, 0, 0);
+
+ PVPBuffer = new G4PVPlacement(0,
+ positionLeadBShield,
+ logicLeadBShield,
+ Name + "_LeadBShield",
+ logicFatimaDetector,
+ false,
+ 0);
+
+
+ // Visualisation of CsIStage Strip
+ G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(1., 1., 0.));
+ logicLeadAShield->SetVisAttributes(LeadVisAtt);
+ logicLeadBShield->SetVisAttributes(LeadVisAtt);
+
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Virtual Method of VDetector class
+
+// Read stream at Configfile to pick-up parameters of detector (Position,...)
+// Called in DetecorConstruction::ReadDetextorConfiguration Method
+void FatimaDetector::ReadConfiguration(string Path)
+{
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+ string LineBuffer, DataBuffer;
+
+ // A:X1_Y1 --> X:1 Y:1
+ // B:X128_Y1 --> X:128 Y:1
+ // C:X1_Y128 --> X:1 Y:128
+ // D:X128_Y128 --> X:128 Y:128
+
+ G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
+ G4ThreeVector A , B , C , D ;
+ G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
+
+ bool ReadingStatus = false;
+
+ bool check_A = false;
+ bool check_C = false;
+ bool check_B = false;
+ bool check_D = false;
+
+ bool check_Theta = false;
+ bool check_Phi = false;
+ bool check_R = false;
+ bool check_beta = false;
+
+ bool checkVis = false;
+
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+ if (LineBuffer.compare(0, 14, "FatimaDetector") == 0) {
+ G4cout << "///" << G4endl ;
+ G4cout << "Detector element found: " << G4endl ;
+ ReadingStatus = true ;
+ }
+
+ while (ReadingStatus) {
+ ConfigFile >> DataBuffer;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;}
+
+ // Position method
+ else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax * mm ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay * mm ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az * mm ;
+
+ A = G4ThreeVector(Ax, Ay, Az);
+ cout << "X1 Y1 corner position : " << A << endl;
+ }
+ else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx * mm ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By * mm ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz * mm ;
+
+ B = G4ThreeVector(Bx, By, Bz);
+ cout << "X128 Y1 corner position : " << B << endl;
+ }
+ else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx * mm ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy * mm ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz * mm ;
+
+ C = G4ThreeVector(Cx, Cy, Cz);
+ cout << "X1 Y128 corner position : " << C << endl;
+ }
+ else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx * mm ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy * mm ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz * mm ;
+
+ D = G4ThreeVector(Dx, Dy, Dz);
+ cout << "X128 Y128 corner position : " << D << endl;
+ }
+
+ // Angle method
+ else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ Theta = Theta * deg;
+ cout << "Theta: " << Theta / deg << endl;
+ }
+ else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ Phi = Phi * deg;
+ cout << "Phi: " << Phi / deg << endl;
+ }
+ else if (DataBuffer.compare(0, 2, "R=") == 0) {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ R = R * mm;
+ cout << "R: " << R / mm << endl;
+ }
+ else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
+ check_beta = true;
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ beta_u = beta_u * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ beta_v = beta_v * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ beta_w = beta_w * deg ;
+ G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
+ }
+
+ else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
+ checkVis = true ;
+ ConfigFile >> DataBuffer;
+ if (DataBuffer.compare(0, 3, "all") == 0) m_non_sensitive_part_visiualisation = true;
+ }
+
+ else G4cout << "WARNING: Wrong Token, FatimaDetector: Detector Element not added" << G4endl;
+
+ // Add The previously define telescope
+ // With position method
+ if ((check_A && check_B && check_C && check_D && checkVis) &&
+ !(check_Theta && check_Phi && check_R)) {
+ ReadingStatus = false;
+ check_A = false;
+ check_C = false;
+ check_B = false;
+ check_D = false;
+ checkVis = false;
+
+ AddModule(A, B, C, D);
+ }
+
+ // With angle method
+ if ((check_Theta && check_Phi && check_R && checkVis) &&
+ !(check_A && check_B && check_C && check_D)) {
+ ReadingStatus = false;
+ check_Theta = false;
+ check_Phi = false;
+ check_R = false;
+ check_beta = false;
+ checkVis = false;
+
+ AddModule(R, Theta, Phi, beta_u, beta_v, beta_w);
+ }
+ }
+ }
+}
+
+
+
+// Construct detector and inialise sensitive part.
+// Called After DetectorConstruction::AddDetector Method
+void FatimaDetector::ConstructDetector(G4LogicalVolume* world)
+{
+ G4RotationMatrix* MMrot = NULL ;
+ G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ;
+
+ G4int NumberOfTelescope = m_DefinitionType.size() ;
+
+ for (G4int i = 0; i < NumberOfTelescope; i++) {
+ // By Point
+ if (m_DefinitionType[i]) {
+ // (u,v,w) unitary vector associated to telescope referencial
+ // (u,v) // to silicon plan
+ // w perpendicular to (u,v) plan and pointing ThirdStage
+ MMu = m_X128_Y1[i] - m_X1_Y1[i];
+ MMu = MMu.unit();
+
+ MMv = m_X1_Y128[i] - m_X1_Y1[i];
+ MMv = MMv.unit();
+
+ // G4double MMscal = MMu.dot(MMv);
+
+ MMw = MMu.cross(MMv);
+// if (MMw.z() > 0) MMw = MMv.cross(MMu) ;
+ MMw = MMw.unit();
+
+ MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4;
+
+ // Passage Matrix from Lab Referential to Telescope Referential
+ // MUST2
+ MMrot = new G4RotationMatrix(MMu, MMv, MMw);
+ // translation to place Telescope
+ MMpos = MMw * Length * 0.5 + MMCenter;
+ }
+
+ // By Angle
+ else {
+ G4double Theta = m_Theta[i] ;
+ G4double Phi = m_Phi[i] ;
+
+ // (u,v,w) unitary vector associated to telescope referencial
+ // (u,v) // to silicon plan
+ // w perpendicular to (u,v) plan and pointing ThirdStage
+ // Phi is angle between X axis and projection in (X,Y) plan
+ // Theta is angle between position vector and z axis
+ G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
+ G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
+ G4double wZ = m_R[i] * cos(Theta / rad);
+ MMw = G4ThreeVector(wX, wY, wZ);
+
+ // vector corresponding to the center of the module
+ G4ThreeVector CT = MMw;
+
+ // vector parallel to one axis of silicon plane
+ G4double ii = cos(Theta / rad) * cos(Phi / rad);
+ G4double jj = cos(Theta / rad) * sin(Phi / rad);
+ G4double kk = -sin(Theta / rad);
+ G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
+
+ MMw = MMw.unit();
+ MMu = MMw.cross(Y);
+ MMv = MMw.cross(MMu);
+ MMv = MMv.unit();
+ MMu = MMu.unit();
+
+ // Passage Matrix from Lab Referential to Telescope Referential
+ // MUST2
+ MMrot = new G4RotationMatrix(MMu, MMv, MMw);
+ // Telescope is rotate of Beta angle around MMv axis.
+ MMrot->rotate(m_beta_u[i], MMu);
+ MMrot->rotate(m_beta_v[i], MMv);
+ MMrot->rotate(m_beta_w[i], MMw);
+ // translation to place Telescope
+ MMpos = MMw * Length * 0.5 + CT ;
+
+ G4cout << "%%%%%%%%%%%%%%" << G4endl;
+ G4cout << MMpos << " " << CT << G4endl;
+ }
+
+ VolumeMaker(i + 1, MMpos, MMrot, world);
+ }
+
+ delete MMrot ;
+}
+
+
+
+// Connect the FatimaData class to the output TTree
+// of the simulation
+void FatimaDetector::InitializeRootOutput()
+{
+}
+
+
+
+// Set the TinteractionCoordinates object from VDetector to the present class
+void FatimaDetector::SetInterCoordPointer(TInteractionCoordinates* interCoord)
+{
+ ms_InterCoord = interCoord;
+}
+
+
+
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAvtion
+void FatimaDetector::ReadSensitive(const G4Event* event)
+{
+ //////////////////////////////////////////////////////////////////////////////////////
+ //////////////////////// Used to Read Event Map of detector //////////////////////////
+ //////////////////////////////////////////////////////////////////////////////////////
+
+ momentum = event->GetPrimaryVertex()->GetPrimary()->GetMomentum();
+ G4double EGamma;
+ EGamma = momentum.getR(); // for photon E=p
+ //G4double EGammaMin = EGamma-4*ResoFirstStage;
+ //G4double EGammaMax = EGamma+4*ResoFirstStage;
+
+
+ // First Stage (LaBr3)
+ std::map<G4int, G4int*>::iterator DetectorNumber_itr;
+ std::map<G4int, G4int*>::iterator CrystalNumber_itr;
+ std::map<G4int, G4double*>::iterator Energy_itr;
+ std::map<G4int, G4double*>::iterator Time_itr;
+ //std::map<G4int, G4double*>::iterator X_itr;
+ //std::map<G4int, G4double*>::iterator Y_itr;
+ //std::map<G4int, G4double*>::iterator Pos_X_itr;
+ //std::map<G4int, G4double*>::iterator Pos_Y_itr;
+ //std::map<G4int, G4double*>::iterator Pos_Z_itr;
+ //std::map<G4int, G4double*>::iterator Ang_Theta_itr;
+ //std::map<G4int, G4double*>::iterator Ang_Phi_itr;
+
+ G4THitsMap<G4int>* DetectorNumberHitMap;
+ G4THitsMap<G4int>* CrystalNumberHitMap;
+ G4THitsMap<G4double>* EnergyHitMap;
+ G4THitsMap<G4double>* TimeHitMap;
+ // G4THitsMap<G4double>* XHitMap;
+ // G4THitsMap<G4double>* YHitMap;
+ //G4THitsMap<G4double>* PosXHitMap;
+ //G4THitsMap<G4double>* PosYHitMap;
+ //G4THitsMap<G4double>* PosZHitMap;
+ //G4THitsMap<G4double>* AngThetaHitMap;
+ //G4THitsMap<G4double>* AngPhiHitMap;
+
+ // NULL pointer are given to avoid warning at compilation
+ // Second Stage (CsI)
+ /*
+ std::map<G4int, G4int*>::iterator CsIDetectorNumber_itr;
+ std::map<G4int, G4int*>::iterator CsICrystalNumber_itr; // added by Marc
+ std::map<G4int, G4double*>::iterator CsIStageEnergy_itr ;
+ G4THitsMap<G4int>* CsIDetectorNumberHitMap = NULL ;
+ G4THitsMap<G4int>* CsICrystalNumberHitMap = NULL ; // added by Marc
+ G4THitsMap<G4double>* CsIStageEnergyHitMap = NULL ;
+
+ */
+
+ // Read the Scorer associate to the LaBr
+ //Detector Number
+ G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/DetectorNumber") ;
+ DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ;
+ DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ;
+
+ //Crystal Number
+ G4int CrystDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/CrystalNumber") ;
+ CrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CrystDetCollectionID)) ;
+ CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin() ;
+
+ //Energy
+ G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripEnergy") ;
+ EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
+ Energy_itr = EnergyHitMap->GetMap()->begin() ;
+
+ //Time of Flight
+ G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripTime") ;
+ TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ;
+ Time_itr = TimeHitMap->GetMap()->begin() ;
+
+ /*
+ //Strip Number X
+ G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripIDFront") ;
+ XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ;
+ X_itr = XHitMap->GetMap()->begin() ;
+
+ //Strip Number Y
+ G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/StripIDBack");
+ YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ;
+ Y_itr = YHitMap->GetMap()->begin() ;
+
+
+ //Interaction Coordinate X
+ G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordX") ;
+ PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ;
+ Pos_X_itr = PosXHitMap->GetMap()->begin() ;
+
+ //Interaction Coordinate Y
+ G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordY") ;
+ PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ;
+ Pos_Y_itr = PosYHitMap->GetMap()->begin() ;
+
+ //Interaction Coordinate Z
+ G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordZ") ;
+ PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ;
+ Pos_Z_itr = PosXHitMap->GetMap()->begin() ;
+
+ //Interaction Coordinate Angle Theta
+ G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordAngTheta") ;
+ AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ;
+ Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ;
+
+ //Interaction Coordinate Angle Phi
+ G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerFatimaDetector/InterCoordAngPhi") ;
+ AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)) ;
+ Ang_Phi_itr = AngPhiHitMap->GetMap()->begin() ;
+ */
+
+ /*
+ // Read the Scorer associate to the SecondStage
+ //CsI Detector Number
+ G4int CsIDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerFatimaDetector/CsIDetectorNumber") ;
+ CsIDetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsIDetCollectionID)) ;
+ CsIDetectorNumber_itr = CsIDetectorNumberHitMap->GetMap()->begin() ;
+
+ //CsI Crystal Number // added by Marc
+ G4int CsICryCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerFatimaDetector/CsICrystalNumber") ; // added by Marc
+ CsICrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICryCollectionID)) ; // added by Anna
+ CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin() ; // added by Anna
+
+ //CsI Energy
+ G4int CsIStageEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerFatimaDetector/CsIStageEnergy") ;
+ CsIStageEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIStageEnergyCollectionID)) ;
+ CsIStageEnergy_itr = CsIStageEnergyHitMap->GetMap()->begin() ;
+ */
+
+ // Check the size of different map
+ G4int sizeN = DetectorNumberHitMap->entries(); // number of objects hit by trackID=1 (can be the same object hit several time)
+ //G4int sizeC = CrystalNumberHitMap->entries();
+ G4int sizeE = EnergyHitMap->entries(); // = number of steps with edep non null
+ G4int sizeT = TimeHitMap->entries();
+ //G4int sizeX = PosXHitMap->entries();
+ //G4int sizeY = PosYHitMap->entries();
+ //G4int sizeX = XHitMap->entries();
+ //G4int sizeY = YHitMap->entries();
+
+ /*
+ G4int sizeNCsI= CsIDetectorNumberHitMap->entries();
+ //G4int sizeCCsI= CsICrystalNumberHitMap->entries();
+ //G4int sizeECsI= CsIStageEnergyHitMap->entries();
+ */
+
+ //sizeC *= 1; // remove warning at compilation added by Marc
+ //sizeECsI *= 1; // remove warning at compilation added by Marc
+ //G4cout <<"SizeN=" << sizeN << endl;
+ //G4cout <<"SizeC=" << sizeC << endl;
+ //G4cout <<"SizeN CsI =" << sizeNCsI << endl;
+ //G4cout <<"SizeE CsI =" << sizeECsI << endl;
+
+ //DetectorNumberHitMap->PrintAllHits();
+
+
+ if (sizeE != sizeT) {
+ G4cout << "No match size FATIMA Event Map: sE:"
+ << sizeE << " sT:" << sizeT << endl ;
+
+ // if (sizeE != sizeX) {
+ //G4cout << "No match size FATIMA Event Map: sE:"
+ //<< sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << endl ;
+ return;
+ }
+
+
+ //G4cout <<"SizeN=" << sizeN << G4endl;
+
+
+ //if at least 1 detector is hit:
+ if(sizeN>0)
+ {
+
+ // Deal with trackID=1:
+ G4int N_first= *(DetectorNumber_itr->second); // ID of first det hit
+ G4int NTrackID = DetectorNumber_itr->first - N_first; // first trackID dealt with (not always =1)
+ G4double E = *(Energy_itr->second);
+ G4double T = *(Time_itr->second);
+ G4int NCryst= *(CrystalNumber_itr->second);
+
+ NCryst *= 1; //remove warning at compilation // added by Marc
+ /*
+ G4cout <<"NTrackID=" << NTrackID << G4endl;
+ G4cout <<"N_first=" << N_first << G4endl;
+ G4cout <<"CrystalNumber_first=" << NCryst << G4endl;
+ G4cout <<"Energy first=" << E << G4endl;
+ G4cout <<"Time first =" << T << G4endl;
+ cout<<"*******"<<endl;
+
+ //added by Nicolas on the model of gaspard scorers 8/7/11
+ CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin();
+
+ for (G4int l = 0 ; l < sizeC ; l++) {
+ G4double NCryst = *(CrystalNumber_itr->second);
+ G4int CTrackID = CrystalNumber_itr->first - N_first - NCryst;
+ G4cout << N_first << " " << NTrackID << " " << CTrackID << " " << NCryst << G4endl;
+ if (CTrackID == NTrackID) {
+ G4cout << "****************** NCryst " << NCryst << G4endl;
+ //ms_Event->SetGPDTrkFirstStageFrontEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ //ms_Event->SetGPDTrkFirstStageBackEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ }
+ CrystalNumber_itr++;
+ }
+ //rewind:
+ for (G4int l = 0 ; l < sizeC ; l++) {
+ CrystalNumber_itr--;
+ }
+ */
+
+ //if there is more than 1 interaction in crystals (1 interaction = 1 hit)
+ if(sizeN>1)
+ {
+
+ //At clusters'level
+ for (G4int l = 1; l < sizeN ; l++) { // loop on all the other tracks
+
+ Energy_itr++;
+ Time_itr++;
+ DetectorNumber_itr++;
+ CrystalNumber_itr++;
+
+ G4int N= *(DetectorNumber_itr->second); // ID of Detector hit
+ G4int Nc= *(CrystalNumber_itr->second); // ID of hit crystal (always =0 for Detector)
+
+ NTrackID = DetectorNumber_itr->first - N; // ID of the track
+
+ //G4cout <<"l=" << l << G4endl;
+ //G4cout <<"N=" << N << G4endl;
+ //G4cout <<"DetectorNumber_itr->first =" << DetectorNumber_itr->first << G4endl;
+ //G4cout <<"NTrackID=" << NTrackID << G4endl;
+
+ if(N==N_first)
+ {
+ // At crystals'level:
+ if(Nc==NCryst) // if we are still in the same crystal than the first hit crystal
+ {
+ E += *(Energy_itr->second);
+
+ }else // we fill the tree for the first detector hit and move to the next detector hit
+ {
+ if(E!=0)
+ {
+ //G4cout <<" -1- filling tree with energy=" << E << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first);
+ ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first);
+ ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst); // added by Marc
+ // Fill Energy
+ // ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ E=RandGauss::shoot(E, ResoFirstStage);
+ ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree
+ //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
+
+ // Fill Time
+ ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
+
+ }
+
+ NCryst=Nc; // continue with the next hit crystal
+ E=*(Energy_itr->second);
+ }
+
+
+ }else // if it's a new detector, one fills the tree with first hit and inspect this new cluster
+ {
+
+ if(E!=0)
+ {
+ //G4cout <<" -2- filling tree with energy=" << E << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first );
+ ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first );
+ ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst); // added by Anna
+ // Fill Energy
+ //ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ E=RandGauss::shoot(E, ResoFirstStage);
+ ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree with energy deposited in first detector
+ //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
+
+ // Fill Time
+ ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
+
+ }
+
+ N_first=N;
+ NCryst=Nc; // and continue with the new hit crystal in this new cluster
+ E=*(Energy_itr->second);
+
+ }
+
+
+ //G4cout <<"Energy=" << E << G4endl;
+ //G4cout <<"Time =" << T << G4endl;
+
+ // Always fill the tree at the end of the loop:
+ if(l==(sizeN-1) && E!=0)
+ {
+ //G4cout <<"-3- filling tree with energy=" << E << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first);
+ ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first);
+ ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst); // added by Anna
+ //cout<<NTrackID<<" filled at the end "<<NCryst<<endl;
+ // Fill Energy
+ // ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ E=RandGauss::shoot(E, ResoFirstStage);
+ ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree
+ //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
+
+ // Fill Time
+ ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
+ }
+
+ }
+
+ }else
+ {
+ // Fill the tree if sizeN=1:
+ if(E!=0)
+ {
+ //G4cout <<" -4- filling tree with energy=" << E << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMALaBr3StageEDetectorNbr(m_index["Detector"] + N_first);
+ ms_Event->SetFATIMALaBr3StageTDetectorNbr(m_index["Detector"] + N_first);
+ ms_Event->SetFATIMALaBr3StageECrystalNbr(NCryst);
+ //cout<<NTrackID<<" filled with sizeN=1 "<<NCryst<<endl;
+ // Fill Energy
+ // ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ E=RandGauss::shoot(E, ResoFirstStage);
+ ms_Event->SetFATIMALaBr3StageEEnergy(E); // Fill the tree
+ //if(E>EGammaMin && E<EGammaMax)ms_Event->SetFATIMALaBr3StageEffphpeak(EGamma);
+
+ // Fill Time
+ ms_Event->SetFATIMALaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
+ }
+ }
+
+
+
+ }
+
+
+ ///////// For CsI stage:
+ //EGammaMin = EGamma-4*ResoSecondStage;
+ //EGammaMax = EGamma+4*ResoSecondStage;
+ /* if(sizeNCsI>0)
+ {
+ // Deal with trackID=1:
+ G4int NCsI_first= *(CsIDetectorNumber_itr->second); // ID of first det hit
+ G4int NCsITrackID = CsIDetectorNumber_itr->first - NCsI_first; // first trackID dealt with (not always =1)
+ G4double E_CsI=*(CsIStageEnergy_itr->second); // Energy second stage
+ G4int NCsICryst= *(CsICrystalNumber_itr->second); // added by Marc
+
+ //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
+ //G4cout <<"NCsI_first=" << NCsI_first << G4endl;
+ //G4cout <<"CsI Energy first=" << E_CsI << G4endl;
+
+ CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin(); // added by Marc
+
+ //for (G4int l = 0 ; l < sizeCCsI ; l++) {
+ // G4double NCsICryst = *(CsICrystalNumber_itr->second);
+ // G4int CCsITrackID = CsICrystalNumber_itr->first - NCsI_first - NCsICryst;
+ // G4cout << NCsI_first << " " << NCsITrackID << " " << CCsITrackID << " " << NCsICryst << G4endl;
+ // if (CCsITrackID == NCsITrackID) {
+ // G4cout << "****************** NCsICryst " << NCsICryst << G4endl;
+ // }
+ // CsICrystalNumber_itr++;
+ // }
+ ////rewind:
+ //for (G4int l = 0 ; l < sizeCCsI ; l++) {
+ // CsICrystalNumber_itr--;
+ //}
+
+
+ if(sizeNCsI>1)
+ {
+ for (G4int l = 1; l < sizeNCsI ; l++) { // loop on all the other tracks
+
+ CsIStageEnergy_itr++;
+ CsIDetectorNumber_itr++;
+ CsICrystalNumber_itr++;
+
+
+ G4int NCsI= *(CsIDetectorNumber_itr->second); // ID of Detector hit
+ G4int NcCsI= *(CsICrystalNumber_itr->second); // ID of crystal hit (always =0 for detector)
+ NCsITrackID = CsIDetectorNumber_itr->first - NCsI; // ID of the track
+
+ //G4cout <<"l=" << l << G4endl;
+ //G4cout <<"NCsI=" << NCsI << G4endl;
+ //G4cout <<"DetectorNumber_itr->first =" << CsIDetectorNumber_itr->first << G4endl;
+ //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
+
+ // At detectors'level:
+ if(NCsI==NCsI_first)
+ {
+ // At crystals'level:
+ if(NcCsI==NCsICryst) // if we are still in the same crystal than the first hit CsI crystal
+ {
+
+ E_CsI += *(CsIStageEnergy_itr->second);
+
+ }else // we fill the tree for the first detector hit and move to the next detector hit
+ {
+ if(E_CsI!=0)
+ {
+ // Fill detector number
+ ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first);
+ ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first);
+ ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst);
+ // Fill Energy
+ //ms_Event->SetFATIMACsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
+ E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
+ ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree
+ //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
+
+ // Fill Time
+ //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
+ }
+
+ NCsICryst=NcCsI;
+ E_CsI=*(CsIStageEnergy_itr->second);
+ }
+
+ }else // if it's a new cluster, one fills the tree with first hit and inspect the new cluster
+ {
+ if(E_CsI!=0)
+ {
+ //G4cout <<" -2- filling tree with energy in CsI =" << E_CsI << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first );
+ ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first );
+ ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst); // added by Marc
+ // Fill Energy
+ //ms_Event->SetFATIMALaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
+ E_CsI=RandGauss::shoot(E_CsI, ResoFirstStage);
+ ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree with energy deposited in first detector
+ //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
+
+ // Fill Time
+ //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
+
+ }
+
+ NCsI_first=NCsI; // continue with the new hit cluster
+ NCsICryst=NcCsI; // and continue with the new hit CsI crystal in this new cluster
+ E_CsI=*(CsIStageEnergy_itr->second);
+
+ }
+
+ // Always fill the tree at the end of the loop:
+
+ if(l==(sizeNCsI-1) && E_CsI!=0)
+ {
+ //G4cout <<" -3- filling tree with energy in CsI =" << E_CsI << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first);
+ ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first);
+ ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst);
+ // Fill Energy
+ // ms_Event->SetFATIMACsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
+ E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
+ ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree
+ //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
+ // Fill Time
+ //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
+
+ }
+
+ }
+
+ }else
+ {
+ // Fill the tree if sizeN=1:
+ if(E_CsI!=0)
+ {
+ //////G4cout <<" -4- filling tree with energy in CsI =" << E_CsI << G4endl;
+ // Fill detector number
+ ms_Event->SetFATIMACsIStageEDetectorNbr(m_index["Detector"] + NCsI_first);
+ ms_Event->SetFATIMACsIStageTDetectorNbr(m_index["Detector"] + NCsI_first);
+ ms_Event->SetFATIMACsIStageECrystalNbr(NCsICryst);
+ // Fill Energy
+ //ms_Event->SetFATIMACsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
+ E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
+ ms_Event->SetFATIMACsIStageEEnergy(E_CsI); // Fill the tree
+ //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetFATIMACsIStageEffphpeak(EGamma);
+ // Fill Time
+ //ms_Event->SetFATIMACsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
+ }
+ }
+
+ }
+
+ */
+
+
+ // clear map for next event
+ DetectorNumberHitMap -> clear();
+ CrystalNumberHitMap -> clear(); // added by Marc
+ EnergyHitMap -> clear();
+ TimeHitMap -> clear();
+ //XHitMap -> clear();
+ //YHitMap -> clear();
+ //PosXHitMap -> clear();
+ //PosYHitMap -> clear();
+ //PosZHitMap -> clear();
+ //AngThetaHitMap -> clear();
+ //AngPhiHitMap -> clear();
+
+ /*
+ CsIDetectorNumberHitMap -> clear();
+ CsICrystalNumberHitMap -> clear(); // added by Marc
+ CsIStageEnergyHitMap -> clear();
+ */
+
+
+}
+
+
+
+void FatimaDetector::InitializeScorers()
+{
+ // LaBr3 Associate Scorer
+ m_LaBr3StageScorer = new G4MultiFunctionalDetector("LaBr3StageScorerFatimaDetector");
+
+ // G4VPrimitiveScorer* DetNbr = new GENERALSCORERS::PSDetectorNumber("DetectorNumber", "FatimaDetector", 0);
+ G4VPrimitiveScorer* DetNbr = new FATIMAScorerLaBr3StageDetectorNumber("DetectorNumber", "FatimaDetector", 0);
+ // G4VPrimitiveScorer* TOF = new GENERALSCORERS::PSTOF("StripTime","FatimaDetector", 0);
+ G4VPrimitiveScorer* TOF = new FATIMAScorerLaBr3StageTOF("StripTime","FatimaDetector", 0);
+ //G4VPrimitiveScorer* InteractionCoordinatesX = new GENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","FatimaDetector", 0);
+ //G4VPrimitiveScorer* InteractionCoordinatesY = new GENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","FatimaDetector", 0);
+ //G4VPrimitiveScorer* InteractionCoordinatesZ = new GENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","FatimaDetector", 0);
+ //G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new GENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","FatimaDetector", 0);
+ //G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new GENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","FatimaDetector", 0);
+
+ G4VPrimitiveScorer* Energy = new FATIMAScorerLaBr3StageEnergy("StripEnergy", "FatimaDetector", 0);
+ G4VPrimitiveScorer* CrystNbr = new FATIMAScorerLaBr3StageCrystal("CrystalNumber", "FatimaDetector", 0);
+
+ // G4VPrimitiveScorer* StripPositionX = new FATIMAcorerLaBr3StageFrontStripDummyShape("StripIDFront", 0, NumberOfStrips);
+ // G4VPrimitiveScorer* StripPositionY = new FATIMAScorerLaBr3StageBackStripDummyShape("StripIDBack", 0, NumberOfStrips);
+
+ //and register it to the multifunctionnal detector
+ m_LaBr3StageScorer->RegisterPrimitive(DetNbr);
+ m_LaBr3StageScorer->RegisterPrimitive(CrystNbr);
+ m_LaBr3StageScorer->RegisterPrimitive(Energy);
+ m_LaBr3StageScorer->RegisterPrimitive(TOF);
+ //m_LaBr3StageScorer->RegisterPrimitive(StripPositionX);
+ //m_LaBr3StageScorer->RegisterPrimitive(StripPositionY);
+ //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesX);
+ //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesY);
+ //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesZ);
+ //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
+ //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);
+
+
+
+ /*
+
+ // Second stage Associate Scorer
+ m_CsIStageScorer = new G4MultiFunctionalDetector("CsIStageScorerFatimaDetector");
+
+ G4VPrimitiveScorer* CsIDetNbr = new FATIMAScorerCsIStageDetectorNumber("CsIDetectorNumber", "FatimaDetector", 0);
+ G4VPrimitiveScorer* CsICryNbr = new FATIMAScorerCsIStageCrystalNumber("CsICrystalNumber", "FatimaDetector", 0); // added by Marc
+ G4VPrimitiveScorer* CsIStageEnergy = new FATIMAScorerCsIStageEnergy("CsIStageEnergy", "FatimaDetector", 0);
+
+ m_CsIStageScorer->RegisterPrimitive(CsIDetNbr);
+ m_CsIStageScorer->RegisterPrimitive(CsICryNbr); // Added by Marc
+ m_CsIStageScorer->RegisterPrimitive(CsIStageEnergy);
+
+ */
+
+ // Add All Scorer to the Global Scorer Manager
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_LaBr3StageScorer);
+ //G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIStageScorer);
+
+}
diff --git a/NPSimulation/Fatima/FatimaDetector.hh b/NPSimulation/Fatima/FatimaDetector.hh
new file mode 100755
index 000000000..d4e8d22a3
--- /dev/null
+++ b/NPSimulation/Fatima/FatimaDetector.hh
@@ -0,0 +1,191 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: Define a phoswich module for the Fatima detector. *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#ifndef FatimaDetector_h
+#define FatimaDetector_h 1
+
+// C++ headers
+#include <vector>
+
+// NPTool header
+#include "FatimaModule.hh"
+#include "TInteractionCoordinates.h"
+
+using namespace CLHEP;
+
+
+
+class FatimaDetector : public FatimaModule
+{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+public:
+ FatimaDetector();
+ virtual ~FatimaDetector();
+
+ ////////////////////////////////////////////////////
+ //////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // By Position Method
+ void AddModule(G4ThreeVector TL,
+ G4ThreeVector BL,
+ G4ThreeVector BR,
+ G4ThreeVector CT);
+
+ // By Angle Method
+ void AddModule(G4double R,
+ G4double Theta,
+ G4double Phi,
+ G4double beta_u,
+ G4double beta_v,
+ G4double beta_w);
+
+ // Effectively construct Volume
+ // Avoid to have two time same code for Angle and Point definition
+ void VolumeMaker(G4int DetecNumber,
+ G4ThreeVector MMpos,
+ G4RotationMatrix* MMrot,
+ G4LogicalVolume* world);
+
+
+ ///////////////////////////////////////////
+ //// Inherite from FatimaModule class /////
+ ///////////////////////////////////////////
+public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(string Path);
+
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world);
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput();
+
+ // Initialize all scorers necessary for the detector
+ void InitializeScorers();
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event);
+
+ // Give the static TInteractionCoordinates from VDetector to the classes
+ // deriving from FatimaModule
+ // This is mandatory since the Fatima*** does not derive from VDetector
+ void SetInterCoordPointer(TInteractionCoordinates* interCoord);
+ TInteractionCoordinates* GetInterCoordPointer() {return ms_InterCoord;};
+
+
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+private:
+ // Interaction Coordinates coming from VDetector through the
+ // SetInteractionCoordinatesPointer method
+ TInteractionCoordinates* ms_InterCoord;
+
+
+
+ // True if Define by Position, False is Define by angle
+ vector<bool> m_DefinitionType ;
+
+ // Used for "By Point Definition"
+ vector<G4ThreeVector> m_X1_Y1 ; // Top Left Corner Position Vector
+ vector<G4ThreeVector> m_X1_Y128 ; // Bottom Left Corner Position Vector
+ vector<G4ThreeVector> m_X128_Y1 ; // Bottom Right Corner Position Vector
+ vector<G4ThreeVector> m_X128_Y128 ; // Center Corner Position Vector
+
+ // Used for "By Angle Definition"
+ vector<G4double> m_R ; // |
+ vector<G4double> m_Theta ; // > Spherical coordinate of Strips Silicium Plate
+ vector<G4double> m_Phi ; // |
+
+ vector<G4double> m_beta_u ; // |
+ vector<G4double> m_beta_v ; // > Tilt angle of the Telescope
+ vector<G4double> m_beta_w ; // |
+
+ G4ThreeVector momentum;
+
+};
+
+
+
+namespace FATIMADETECTOR
+{
+
+
+
+ // Resolution
+// const G4double ResoFirstStage = 0; // = 50 keV of Resolution // Unit is MeV/2.35
+ const G4double ResoFirstStage = 0.0099; // = 3.5% at .662MeV of Resolution // Unit is MeV/2.35
+ // const G4double ResoSecondStage = 0.0; // = 13% at .662 MeV of resolution // Unit is MeV/2.35
+ // const G4double ResoThirdStage = 0.0; // = 50 keV of resolution // Unit is MeV/2.35
+ const G4double ResoTimeGpd = 0.212765957;// = 500ps // Unit is ns/2.35
+
+ // Geometry for the mother volume containing the different layers of your dummy shape module
+ const G4double FaceFront = 7.5*cm; //;
+ const G4double Length = 26.33*cm; //
+
+ // Geometry for the phoswitch
+ // LaBr3
+ const G4double LaBr3Face = 3.81*cm;
+ const G4double LaBr3Thickness = 5.08*cm; // for detector 2x2x2
+
+ // Al Can
+ const G4double CanOuterDiameter = 4.3*cm;
+ const G4double CanInnerDiameter = 4.15*cm;
+ const G4double CanLength = 4.33*cm; // = (5.08-(22-0.85)+0.1 with 0.85=(26.33-0.1-5.08)-22
+ // Al front Window
+ const G4double WinOuterDiameter = 4.15*cm;
+ const G4double WinInnerDiameter = 0*cm;
+ const G4double WinLength = 0.1*cm; //
+
+ // PM (TODO Change CsI to Alu)
+ const G4double CsIFace = FaceFront;
+ const G4double CsIThickness = 22.*cm; // for detector
+ //const G4double CsIThickness = 0.001*cm; // For LaBr long
+
+ // Lead tube (TODO Change CsI to Alu)
+ const G4double LeadAMinR = CanOuterDiameter;
+ const G4double LeadAMaxR = CanOuterDiameter+1*cm; // for detector
+ const G4double LeadALength = 4.33*cm; // For LaBr long
+
+ const G4double LeadBMinR = CanOuterDiameter;
+ const G4double LeadBMaxR = FaceFront; // for detector
+ const G4double LeadBLength = 1.*cm; // For LaBr long
+
+ // Position
+ const G4double LaBr3Stage_PosZ = -Length*0.5 + 0.5*LaBr3Thickness + 0.1*cm;
+ const G4double LaBr3Can_PosZ = -Length*0.5 + 0.5*CanLength;
+ const G4double LaBr3Win_PosZ = -Length*0.5 + 0.5*WinLength;
+
+ const G4double CsIStage_PosZ = -Length*0.5 + (Length-CsIThickness) + 0.5*CsIThickness;
+
+ const G4double LeadAShield_PosZ = -Length*0.5 + 0.5*LeadALength;
+ const G4double LeadBShield_PosZ = -Length*0.5 + LeadALength - 0.5*LeadBLength ;
+
+}
+
+#endif
diff --git a/NPSimulation/Fatima/FatimaModule.cc b/NPSimulation/Fatima/FatimaModule.cc
new file mode 100755
index 000000000..b53bfa426
--- /dev/null
+++ b/NPSimulation/Fatima/FatimaModule.cc
@@ -0,0 +1,60 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/14 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class is an Abstract Base Class (ABC) from which should *
+ * derive all different modules from the Fatima detector. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#include "FatimaModule.hh"
+#include "RootOutput.h"
+
+
+TFatimaData *FatimaModule::ms_Event = 0;
+
+
+
+FatimaModule::FatimaModule()
+{
+ if (ms_Event == 0) ms_Event = new TFatimaData();
+
+ InitializeRootOutput();
+ InitializeIndex();
+}
+
+
+
+FatimaModule::~FatimaModule()
+{
+}
+
+
+
+void FatimaModule::InitializeRootOutput()
+{
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ // if the branch does not exist yet, create it
+ if (!pTree->GetBranch("FATIMA"))
+ pTree->Branch("FATIMA", "TFatimaData", &ms_Event);
+}
+
+
+
+void FatimaModule::InitializeIndex()
+{
+ m_index["LaBr"] = 0; // 24 LaBr max
+}
diff --git a/NPSimulation/Fatima/FatimaModule.hh b/NPSimulation/Fatima/FatimaModule.hh
new file mode 100755
index 000000000..45bd3a5b4
--- /dev/null
+++ b/NPSimulation/Fatima/FatimaModule.hh
@@ -0,0 +1,94 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: This class is an Abstract Base Class (ABC) from which should *
+ * derive all different modules from the Fatima detector. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#ifndef FatimaModule_h
+#define FatimaModule_h 1
+
+// C++ headers
+#include <string>
+#include <vector>
+
+// G4 headers
+#include "G4LogicalVolume.hh"
+#include "G4Event.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+// NPTool - ROOT headers
+#include "TInteractionCoordinates.h"
+#include "TFatimaData.h"
+
+using namespace std;
+
+
+
+class FatimaModule
+{
+public:
+ FatimaModule();
+ virtual ~FatimaModule();
+
+public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ virtual void ReadConfiguration(string Path) = 0;
+
+ // Construct detector and inialise sensitive part.
+ virtual void ConstructDetector(G4LogicalVolume* world) = 0;
+
+ // Read sensitive part of a the FatimaModule detector and fill the Root tree.
+ virtual void ReadSensitive(const G4Event* event) = 0;
+
+ // Add Detector branch to the ROOT output tree
+ virtual void InitializeRootOutput();
+
+ // Initialize all scorers necessary for each detector
+ virtual void InitializeScorers() = 0;
+
+ // Give the static TInteractionCoordinates from VDetector to the classes
+ // deriving from FatimaModule
+ // This is mandatory since the Fatima*** does not derive from VDetector
+ virtual void SetInterCoordPointer(TInteractionCoordinates* interCoord) = 0;
+ virtual TInteractionCoordinates* GetInterCoordPointer() = 0;
+
+ // Initialize the Index map for the different modules of Fatima
+ void InitializeIndex();
+
+public:
+ TFatimaData* GetEventPointer() {return ms_Event;};
+
+protected:
+ // Class to store the result of simulation
+ static TFatimaData* ms_Event;
+
+ // Set to true if you want to see Telescope Frame in your visualisation
+ bool m_non_sensitive_part_visiualisation;
+
+protected:
+ // First stage Associate Scorer
+ G4MultiFunctionalDetector* m_LaBr3StageScorer;
+
+ // Second stage Associate Scorer
+ G4MultiFunctionalDetector* m_CsIStageScorer;
+
+protected:
+ map<string, int> m_index;
+};
+
+#endif
diff --git a/NPSimulation/Fatima/FatimaScorers.cc b/NPSimulation/Fatima/FatimaScorers.cc
new file mode 100755
index 000000000..a869412d6
--- /dev/null
+++ b/NPSimulation/Fatima/FatimaScorers.cc
@@ -0,0 +1,514 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : 02/07/13 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class holds all the scorers needed by the *
+ * Fatima*** objects. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+// G4 headers
+#include "G4UnitsTable.hh"
+
+// NPTool headers
+#include "ObsoleteGeneralScorers.hh"
+#include "FatimaScorers.hh"
+
+#include "FatimaDetector.hh"
+
+using namespace FATIMADETECTOR;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Added by Adrien MATTA:
+// Those Scorer use TrackID as map index. This way ones can rebuild energy deposit,
+// time of flight or position,... particle by particle for each event. Because standard
+// scorer provide by G4 don't work this way but using a global ID for each event you should
+// not use those scorer with some G4 provided ones or being very carefull doing so.
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// FirstStage Detector Number Scorer (deal with multiple particle hit)
+FATIMAScorerLaBr3StageDetectorNumber::FATIMAScorerLaBr3StageDetectorNumber(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerLaBr3StageDetectorNumber::~FATIMAScorerLaBr3StageDetectorNumber()
+{
+}
+
+G4bool FATIMAScorerLaBr3StageDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+ // int DetNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0); // this line does exactly the same than the line above
+ int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+
+ // get energy
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ EvtMap->set(CrystalNbr + DetNbr + index, DetNbr);
+ return TRUE;
+}
+void FATIMAScorerLaBr3StageDetectorNumber::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerLaBr3StageDetectorNumber::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerLaBr3StageDetectorNumber::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerLaBr3StageDetectorNumber::DrawAll()
+{
+}
+
+void FATIMAScorerLaBr3StageDetectorNumber::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+ std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
+ for (; itr != EvtMap->GetMap()->end(); itr++) {
+ G4cout << " copy no.: " << itr->first
+ << " Detector Number : " << G4BestUnit(*(itr->second), "DetectorNumber")
+ << G4endl;
+ }
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// FirstStage Energy Scorer (deal with multiple particle hit)
+FATIMAScorerLaBr3StageEnergy::FATIMAScorerLaBr3StageEnergy(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerLaBr3StageEnergy::~FATIMAScorerLaBr3StageEnergy()
+{
+}
+
+G4bool FATIMAScorerLaBr3StageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+ int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+
+ // get energy
+ G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
+ POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
+
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ // if (edep < 100*keV) return FALSE;
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ EvtMap->add(CrystalNbr + DetNbr + index, edep);
+ return TRUE;
+}
+
+void FATIMAScorerLaBr3StageEnergy::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerLaBr3StageEnergy::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerLaBr3StageEnergy::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerLaBr3StageEnergy::DrawAll()
+{
+}
+
+void FATIMAScorerLaBr3StageEnergy::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+ std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
+ for (; itr != EvtMap->GetMap()->end(); itr++) {
+ G4cout << " copy no.: " << itr->first
+ << " energy deposit: " << G4BestUnit(*(itr->second), "Energy")
+ << G4endl;
+ }
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// FirstStage CrystalNbr Scorer (deal with multiple particle hit)
+FATIMAScorerLaBr3StageCrystal::FATIMAScorerLaBr3StageCrystal(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerLaBr3StageCrystal::~FATIMAScorerLaBr3StageCrystal()
+{
+}
+
+G4bool FATIMAScorerLaBr3StageCrystal::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+
+ //G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1);
+ //Adde by Anna:
+ G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); // daughter crystal volume
+ //depth:0 is the world, 1 is the mother...
+ //AC instead of GetReplicaNumber(1)
+ //cout<<"****replica****"<<endl;
+ //cout<<"copy nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1)<<endl;
+ //cout<<"replica nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1)<<endl;
+ //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2)<<endl;
+ //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0)<<endl;
+
+ // get energy
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ // if (edep < 100*keV) return FALSE;
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ //EvtMap->add(DetNbr + index, CrystalNbr); // Marc
+ //cout<<"index "<<index<<endl;
+ //cout<<"DetNbr + index "<<DetNbr + index<<endl;
+ EvtMap->set(CrystalNbr + DetNbr + index, CrystalNbr);// modified by Anna
+ return TRUE;
+
+}
+
+void FATIMAScorerLaBr3StageCrystal::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerLaBr3StageCrystal::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerLaBr3StageCrystal::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerLaBr3StageCrystal::DrawAll()
+{
+}
+
+void FATIMAScorerLaBr3StageCrystal::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+ std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
+ for (; itr != EvtMap->GetMap()->end(); itr++) {
+ G4cout << " copy no.: " << itr->first
+ << " crystal nbr: " << G4BestUnit(*(itr->second), "Crystal")
+ << G4endl;
+ }
+}
+
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// FirstStage ToF Scorer (deal with multiple particle hit)
+FATIMAScorerLaBr3StageTOF::FATIMAScorerLaBr3StageTOF(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerLaBr3StageTOF::~FATIMAScorerLaBr3StageTOF()
+{
+}
+
+G4bool FATIMAScorerLaBr3StageTOF::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+ int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+
+ // get TOF
+ G4double TOF = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ // if (edep < 100*keV) return FALSE;
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ EvtMap->add(CrystalNbr + DetNbr + index, TOF);
+ return TRUE;
+}
+
+void FATIMAScorerLaBr3StageTOF::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerLaBr3StageTOF::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerLaBr3StageTOF::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerLaBr3StageTOF::DrawAll()
+{
+}
+
+void FATIMAScorerLaBr3StageTOF::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+}
+
+
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// SecondStage (CsI) Detector Number Scorer (deal with multiple particle hit)
+FATIMAScorerCsIStageDetectorNumber::FATIMAScorerCsIStageDetectorNumber(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerCsIStageDetectorNumber::~FATIMAScorerCsIStageDetectorNumber()
+{
+}
+
+G4bool FATIMAScorerCsIStageDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+ // int DetNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0); // this line does exactly the same than the line above
+ int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+
+ // get energy
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ EvtMap->set(CrystalNbr + DetNbr + index, DetNbr);
+ return TRUE;
+}
+void FATIMAScorerCsIStageDetectorNumber::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerCsIStageDetectorNumber::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerCsIStageDetectorNumber::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerCsIStageDetectorNumber::DrawAll()
+{
+}
+
+void FATIMAScorerCsIStageDetectorNumber::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+ std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
+ for (; itr != EvtMap->GetMap()->end(); itr++) {
+ G4cout << " copy no.: " << itr->first
+ << " Detector Number : " << G4BestUnit(*(itr->second), "DetectorNumber")
+ << G4endl;
+ }
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// CsIStage Energy Scorer (deal with multiple particle hit)
+FATIMAScorerCsIStageEnergy::FATIMAScorerCsIStageEnergy(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerCsIStageEnergy::~FATIMAScorerCsIStageEnergy()
+{
+}
+
+G4bool FATIMAScorerCsIStageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+ int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+
+ // get energy
+ G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
+ POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
+
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ //if (edep < 100*keV) return FALSE;
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ EvtMap->add(CrystalNbr + DetNbr + index, edep);
+ return TRUE;
+}
+
+void FATIMAScorerCsIStageEnergy::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerCsIStageEnergy::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerCsIStageEnergy::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerCsIStageEnergy::DrawAll()
+{
+}
+
+void FATIMAScorerCsIStageEnergy::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+ std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
+ for (; itr != EvtMap->GetMap()->end(); itr++) {
+ G4cout << " copy no.: " << itr->first
+ << " energy deposit: " << G4BestUnit(*(itr->second), "Energy")
+ << G4endl;
+ }
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// all added by Anna:
+
+// FirstStage CrystalNbr Scorer (deal with multiple particle hit)
+FATIMAScorerCsIStageCrystalNumber::FATIMAScorerCsIStageCrystalNumber(G4String name, G4String volumeName, G4int depth)
+ : G4VPrimitiveScorer(name, depth), HCID(-1)
+{
+ m_VolumeName = volumeName;
+}
+
+FATIMAScorerCsIStageCrystalNumber::~FATIMAScorerCsIStageCrystalNumber()
+{
+}
+
+G4bool FATIMAScorerCsIStageCrystalNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // get detector number
+ int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
+
+ // get energy
+ G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+ //depth:0 is the world, 1 is the mother...
+ //AC instead of GetReplicaNumber(1)
+ //cout<<"****replica****"<<endl;
+ //cout<<"copy nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1)<<endl;
+ //cout<<"replica nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1)<<endl;
+ //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2)<<endl;
+ //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0)<<endl;
+
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ // if (edep < 100*keV) return FALSE;
+ if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
+ G4int index = aStep->GetTrack()->GetTrackID();
+ //cout<<"index "<<index<<endl;
+ //cout<<"DetNbr + index "<<DetNbr + index<<endl;
+ EvtMap->set(CrystalNbr + DetNbr + index, CrystalNbr);
+ return TRUE;
+
+}
+
+void FATIMAScorerCsIStageCrystalNumber::Initialize(G4HCofThisEvent* HCE)
+{
+ EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+void FATIMAScorerCsIStageCrystalNumber::EndOfEvent(G4HCofThisEvent*)
+{
+}
+
+void FATIMAScorerCsIStageCrystalNumber::Clear()
+{
+ EvtMap->clear();
+}
+
+void FATIMAScorerCsIStageCrystalNumber::DrawAll()
+{
+}
+
+void FATIMAScorerCsIStageCrystalNumber::PrintAll()
+{
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl;
+ std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
+ for (; itr != EvtMap->GetMap()->end(); itr++) {
+ G4cout << " copy no.: " << itr->first
+ << " crystal nbr: " << G4BestUnit(*(itr->second), "Crystal")
+ << G4endl;
+ }
+}
diff --git a/NPSimulation/Fatima/FatimaScorers.hh b/NPSimulation/Fatima/FatimaScorers.hh
new file mode 100755
index 000000000..2db2cb6ca
--- /dev/null
+++ b/NPSimulation/Fatima/FatimaScorers.hh
@@ -0,0 +1,198 @@
+/*****************************************************************************
+ * Copyright (C) 2009 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * *
+ * Creation Date : 04/01/13 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: This class holds all the scorers needed by the *
+ * GaspardTracker*** objects. *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+
+#ifndef FATIMAScorer_h
+#define FATIMAScorer_h 1
+
+#include "G4VPrimitiveScorer.hh"
+#include "G4THitsMap.hh"
+
+namespace FATIMASCORERS
+{
+
+
+class FATIMAScorerLaBr3StageDetectorNumber : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerLaBr3StageDetectorNumber(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerLaBr3StageDetectorNumber();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4int>* EvtMap;
+};
+
+
+
+class FATIMAScorerLaBr3StageEnergy : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerLaBr3StageEnergy(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerLaBr3StageEnergy();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4double>* EvtMap;
+};
+
+
+
+class FATIMAScorerLaBr3StageCrystal : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerLaBr3StageCrystal(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerLaBr3StageCrystal();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4int>* EvtMap;
+};
+
+
+class FATIMAScorerLaBr3StageTOF : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerLaBr3StageTOF(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerLaBr3StageTOF();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4double>* EvtMap;
+};
+
+
+
+class FATIMAScorerCsIStageDetectorNumber : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerCsIStageDetectorNumber(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerCsIStageDetectorNumber();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4int>* EvtMap;
+};
+
+// Added by Anna
+class FATIMAScorerCsIStageCrystalNumber : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerCsIStageCrystalNumber(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerCsIStageCrystalNumber();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4int>* EvtMap;
+};
+
+
+class FATIMAScorerCsIStageEnergy : public G4VPrimitiveScorer
+{
+public: // with description
+ FATIMAScorerCsIStageEnergy(G4String name, G4String volumeName, G4int depth = 0);
+ virtual ~FATIMAScorerCsIStageEnergy();
+
+protected: // with description
+ virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+public:
+ virtual void Initialize(G4HCofThisEvent*);
+ virtual void EndOfEvent(G4HCofThisEvent*);
+ virtual void Clear();
+ virtual void DrawAll();
+ virtual void PrintAll();
+
+private:
+ G4String m_VolumeName;
+ G4int HCID;
+ G4THitsMap<G4double>* EvtMap;
+};
+
+}
+
+using namespace FATIMASCORERS;
+#endif
diff --git a/NPSimulation/src/DetectorConstruction.cc b/NPSimulation/src/DetectorConstruction.cc
index 6c86ace87..784355767 100644
--- a/NPSimulation/src/DetectorConstruction.cc
+++ b/NPSimulation/src/DetectorConstruction.cc
@@ -73,6 +73,10 @@
#include "Eurogam.hh"
#endif
+#ifdef INC_FATIMA
+#include "Fatima.hh"
+#endif
+
#ifdef INC_GASPARD
#include "GaspardTracker.hh"
#endif
@@ -202,10 +206,12 @@ G4VPhysicalVolume* DetectorConstruction::ReadConfigurationFile(){
bool cComptonTelescope = false;
bool cDummy = false;
bool cEurogam = false;
+ bool cFatima = false;
bool cGeneralTarget = false;
bool cGeneralChamber = false;
bool cGPDTracker = false;
- bool cHYD2Tracker = false;
+ bool cHYD2Tracker = false;
+ bool cHelios = false;
bool cMUST2 = false;
bool cPlastic = false;
bool cParis = false;
@@ -215,7 +221,6 @@ G4VPhysicalVolume* DetectorConstruction::ReadConfigurationFile(){
bool cTigress = false;
bool cTiara = false;
bool cW1 = false;
- bool cHelios = false;
int VerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
@@ -301,6 +306,26 @@ G4VPhysicalVolume* DetectorConstruction::ReadConfigurationFile(){
#endif
}
+ ////////////////////////////////////////////
+ //////////// Search for FATIMA ///////////
+ ////////////////////////////////////////////
+ else if (LineBuffer.compare(0, 6, "Fatima") == 0 && cFatima == false) {
+#ifdef INC_FATIMA
+ cFatima = true ;
+ G4cout << "//////// Fatima ////////" << G4endl ;
+
+ // Instantiate the new array as a VDetector Object
+ VDetector* myDetector = new Fatima() ;
+
+ // Read Position of Telescope
+ ConfigFile.close() ;
+ myDetector->ReadConfiguration(Path) ;
+ ConfigFile.open(Path.c_str()) ;
+
+ // Add array to the VDetector Vector
+ AddDetector(myDetector) ;
+#endif
+ }
////////////////////////////////////////////
//////////// Search for Gaspard ////////////
--
GitLab