BmnRoot
Loading...
Searching...
No Matches
BmnSiliconTrackFinder.h
Go to the documentation of this file.
1// Author: Vasilisa Lenivenko (VBLHEP) <vasilisa@jinr.ru> 2019-07-18
2
4// //
5// BmnSiliconTrackFinder //
6// //
7// Implementation of an algorithm developed by //
8// V.Lenivenko and V.Paltchik //
9// to the BmnRoot software //
10// //
11// The algorithm serves for searching for track segments //
12// in the Silicon detectors of the BM@N experiment //
13// //
15
16
17#ifndef BMNSILICONTRACKFINDER_H
18#define BMNSILICONTRACKFINDER_H 1
19
20#include "BmnTrack.h"
21#include "BmnSiliconDigit.h"
22#include "BmnSiliconTrack.h"
23#include "BmnSiliconHit.h"
24
25#include "FairTask.h"
26#include "FairTrackParam.h"
27#include "FairMCPoint.h"
28
29#include <TMath.h>
30#include <TString.h>
31#include <TClonesArray.h>
32#include <TVector2.h>
33#include "TList.h"
34#include "TH1D.h"
35#include "TH2D.h"
36
37#include <fstream>
38#include <vector>
39
40using namespace std;
41using namespace TMath;
42
43struct tracksX {
44 Int_t Nhits = 0;
45 Double_t Chi2 = 0.;
46 Double_t Xv = 999.;
47 Double_t Yv = 999.;
48 Double_t param[4] = {999., 999., 999., 999.};
49 Double_t CoordX[4] = {-999., -999., -999., -999.};
50 Double_t CoordXp[4] = {-999., -999., -999., -999.};
51 Double_t CoordY[4] = {-999., -999., -999., -999.};
52 Double_t CoordZ[4] = {-999., -999., -999., -999.};
53 Double_t SigmaX[4] = {-999., -999., -999., -999.};
54 Double_t SigmaXp[4] = {-999., -999., -999., -999.};
55 Double_t SigmaY[4] = {-999., -999., -999., -999.};
56
57 Int_t ModNum[4] = {-1, -1, -1, -1};
58 Double_t StripNumX[4] = {-999., -999., -999., -999.};
59 Double_t StripNumXp[4]= {-999., -999., -999., -999.};
60 Double_t ClSizeX[4] = {-999., -999., -999., -999.};
61 Double_t ClSizeXp[4] = {-999., -999., -999., -999.};
62 Double_t SumQX[4] = {-999., -999., -999., -999.};
63 Double_t SumQXp[4] = {-999., -999., -999., -999.};
64 Int_t Pdg = -1;
65};
66
67
68class BmnSiliconTrackFinder : public FairTask
69{
70 public:
71 BmnSiliconTrackFinder();
72 BmnSiliconTrackFinder(Bool_t, Int_t, Int_t);
73 virtual ~BmnSiliconTrackFinder();
74
75 virtual InitStatus Init();
76 virtual void Exec(Option_t* opt);
77 virtual void Finish();
78
79 static bool compareSegments(const tracksX &a, const tracksX &b){
80 return a.Chi2 < b.Chi2;
81 }
82
83 struct MC_points{
84 Int_t Id = -1;//mcId
85 Int_t Np = -1;
86 Int_t Pdg = -999;
87 Double_t x[4] = {-999., -999.,-999.,-999.};
88 Double_t y[4] = {-999., -999.,-999.,-999.};
89 Double_t z[4] = {-999., -999.,-999.,-999.};
90 Bool_t wo3st = 0;
91
92 Double_t param[4] = {999., 999., 999., 999.};
93 };
94
95 private:
96 Bool_t fDebug = 0;
97 UInt_t fEventNo = 0; // event counter
98 Int_t fRunNumber;
99 Int_t fRunPeriod;
100 Int_t N;
101 Bool_t expData;
102 TList fList;
103 TClonesArray* fBmnSiDigitsArray;
104 TClonesArray* fSiTracks;
106 TString fOutputHitsBranchName;
107 TClonesArray* fBmnSiliconHitsArray;
108 TClonesArray* fBmnSiliconHitsXArray;
109 TClonesArray* fBmnSiliconHitsXpArray;
110 TClonesArray* fBmnSiliconHitsXYArray;
111 TClonesArray* fBmnHitsArray;
112 TClonesArray* fBmnHitsArray2;
113 TString fOutputFileName;
114 TString fInputBranchName;
115 TString fInputBranchName2;
116
117 //--
118 Int_t fNstations;
119 Int_t fNmodules;
120 Int_t fNmodules1;
121 Int_t fNmodules2;
122 Int_t fNmodules3;
123 Double_t half_target_regionX;
124 Double_t half_target_regionY;
125 Double_t ChiSquare;
126 Double_t Z0_SRC_target;
127 Double_t Zcentr;
128 Double_t kX_target;
129 Double_t kY_target;
130 Int_t *iClustXMod;
131 Int_t *iModX;
132 Int_t *Nsp_st;
133 Double_t *AmatrInverted;
134 Double_t *rhs;
135 Double_t *AmatrArray;
136 Double_t *line;
137 Double_t *shiftStXtoGlob;
138 Double_t *shiftStYtoGlob;
139 Double_t *Zposition;
140 Double_t *SensitiveAreaY;
141 Double_t *Ampl_strX;
142 Double_t *Ampl_strXp;
143 Int_t **NclustX;
144 Int_t **NclustXp;
145 Int_t **NhitsXYmod;
146 Int_t **Nleftoversp;
147 Int_t **NleftoverX;
148 Int_t **NleftoverXp;
149 Double_t **Amatr;
150 Double_t **shiftX;
151 Double_t **shiftY;
152 Double_t **shiftYbelow;
153 Double_t **Zstation;
154 Double_t **Xforglfit;
155 Double_t **Xpforglfit;
156 Double_t **Yforglfit;
157 Double_t **SigmXforglfit;
158 Double_t **SigmXpforglfit;
159 Double_t **SigmYforglfit;
160 Double_t **ZstnCforglfit;
161 Double_t ***DigitsArrayX;
162 Double_t ***DigitsArrayXp;
163 Double_t ***XClust;
164 Double_t ***XpClust;
165 Double_t ***XClust_width;
166 Double_t ***XpClust_width;
167 Double_t ***sumQx;
168 Double_t ***sumQxp;
169 Double_t ***XspClust;
170 Double_t ***XpspClust;
171 Double_t ***XspClust_width;
172 Double_t ***XpspClust_width;
173 Double_t ***sumQxsp;
174 Double_t ***sumQxpsp;
175 Double_t ***XCoord;
176 Double_t ***XpCoord;
177 Double_t ***XspCoord;
178 Double_t ***XpspCoord;
179 Double_t ***YspCoord;
180 Int_t ***Sp_pdg;
181 Double_t ***SigmaX;
182 Double_t ***SigmaXp;
183 Double_t ***SigmspX;
184 Double_t ***SigmspXp;
185 Double_t ***SigmspY;
186 Double_t ***leftoverXsp;
187 Double_t ***leftoverXpsp;
188 Double_t ***leftoverYsp;
189 Double_t ***leftoverXsigsp;
190 Double_t ***leftoverXpsigsp;
191 Double_t ***leftoverYsigsp;
192 Double_t ***leftoverX;
193 Double_t ***leftoverXp;
194 Double_t ***leftoverXsig;
195 Double_t ***leftoverXpsig;
196
197 vector<tracksX> vec_tracks;
198 vector<tracksX> CleanTr;
199 vector<MC_points> vec_points;
200
201 TH1D *hNpoint, *hChiSquareNdf, *hNtracks, *hAxglob, *hBxglob, *hAyglob,*hByglob, *hY_st_1,*hY_st_2,*hY_st_3, *hdY_st_1,*hdY_st_2,*hdY_st_3,
202 *hX_st_1,*hX_st_2,*hX_st_3, *hdX_st_1,*hdX_st_2,*hdX_st_3, *hdXp_st_1,*hdXp_st_2,*hdXp_st_3,
203 *N_eff,* D_eff,* E_eff, *hdX_st1_st2, *hdX_st2_st3, *hdX_st1_st3,*hdY_st1_st2, *hdY_st2_st3, *hdY_st1_st3,*hdXst1_0_st2_0,*hdXst1_0_st2_1,*hdXst1_1_st2_0,*hdXst1_1_st2_1,
204 *hdXst1_2_st2_0,*hdXst1_2_st2_1,*hdXst1_3_st2_0,*hdXst1_3_st2_1,*hdXst2_0_st3_1,*hdXst2_0_st3_2,*hdXst2_1_st3_1,*hdXst2_1_st3_2,*hdYst1_0_st2_0,
205 *hdYst1_0_st2_1,*hdYst1_1_st2_0,*hdYst1_1_st2_1,*hdYst1_2_st2_0,*hdYst1_2_st2_1,*hdYst1_3_st2_0,*hdYst1_3_st2_1,*hdYst2_0_st3_1,*hdYst2_0_st3_2,
206 *hdYst2_1_st3_1,*hdYst2_1_st3_2,*hX13_X2_m0,*hX13_X2_m1, *hXp13_Xp2_m0,*hXp13_Xp2_m1, *hY13_Y2_m0,*hY13_Y2_m1, *hY1m0_Y23, *hY1m1_Y23, *hY1m2_Y23, *hY1m3_Y23,
207 *hAx_first_tr, *hAx_more_first_tr, *hdAx_MC_tr, *hdAy_MC_tr,*hdX_MC_tr, *hdY_MC_tr,
208 *hAxsi_mctrue,*hAysi_mctrue,*hBxsi_mctrue,*hBysi_mctrue,
209 *hdAx_MC_tr_comb, *hdAy_MC_tr_comb, *hdX_MC_tr_comb, *hdY_MC_tr_comb,
210 *hDen_mctrSi, *hNum_mctrSi,*hEff_mctrSi, *hNtrsi_mc, *hNtrsi_reco,
211 *hDen_mcreaction, *hNum_mcreaction,*hEff_mcreaction,
212 *hdXp3_mod1, *hdXp3_mod2, *hdXXp3_mod1, *hdXXp3_mod2, *hXXp12CheckLeftover, *hXXp12CheckLeftover03,
213 *hXSi1_run8, *hXSi2_run8, *hYSi1_run8, *hYSi2_run8;
214 TH2D *hvertexXY, * hvertex_aver_XY,* hprofile_beam_z1,* hprofile_beam_z2,* hprofile_beam_z3, *hdYvsYst_1,*hdYvsYst_2,*hdYvsYst_3,
215 *hdXvsXst_1,*hdXvsXst_2,*hdXvsXst_3,
216 *hdYvsYst1_mod0, *hdYvsYst1_mod1,*hdYvsYst1_mod2,*hdYvsYst1_mod3,*hdYvsYst2_mod0,*hdYvsYst2_mod1,*hdYvsYst3_mod1,*hdYvsYst3_mod2,
217 *hdXvsXst1_0,*hdXvsXst1_1,*hdXvsXst1_2,*hdXvsXst1_3,*hdXvsXst2_0,*hdXvsXst2_1,*hdXvsXst3_1,*hdXvsXst3_2,
218 *hSi_st3mc, *hSi_st2mc, * hSi_st1mc, *hNtrsi_mc_vs_reco, *hY_XSisp1, *hY_XSisp2, *hY_XSisp3, *hY_XSisp1bef, *hY_XSisp2bef,
219 *hYXSi1_run8, *hYXSi2_run8;
220 vector<TH1D*> hoccupancyX1, hoccupancyX2, hoccupancyX3, hoccupancyXp1, hoccupancyXp2, hoccupancyXp3;
221
222 void StripsReading(Double_t ***, Double_t ***);
223 void PrepareArraysToProcessEvent();
224 void Case1(Int_t **, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, vector<tracksX>&,
225 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***);
226 void Case2(Int_t & , Int_t **, Double_t ***, Double_t ***, Double_t ***,Double_t ***, Double_t ***, Double_t ***, Int_t **, Int_t **,Double_t ***, Double_t ***, Double_t ***, Double_t ***, vector<tracksX>&,
227 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***,
228 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***);
229 void Case3(Int_t **, Double_t ***, Double_t ***, Double_t ***,Double_t ***, Double_t ***, Double_t ***, Int_t **, Int_t **,Double_t ***, Double_t ***, Double_t ***, Double_t ***, vector<tracksX>&,
230 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***,
231 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***);
232 void Case4(Int_t & , Int_t **, Double_t ***, Double_t ***, Double_t ***,Double_t ***, Double_t ***, Double_t ***, Int_t **, Int_t **,Double_t ***, Double_t ***, Double_t ***, Double_t ***, vector<tracksX>&,
233 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***,
234 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***);
235
236 void Clustering( Double_t***, Double_t *** , Double_t *** , Double_t *** , Double_t *** , Int_t **, Double_t *** , Double_t *** , Double_t *** , Int_t **);
237 void GlobalFit( Double_t **, Double_t **, Double_t **,Double_t **, Double_t **, Double_t **, Double_t *);
238 bool InvertMatrix(Double_t *, Double_t *);
239 void GetXYspatial(Int_t **NClX, Int_t **NClXp, Double_t ***XCoor, Double_t ***XpCoor, Double_t ***Xsp, Double_t ***Xpsp, Double_t ***Ysp, Int_t **NXYsp, Double_t ***SigmX, Double_t ***SigmXp, Double_t ***SigspX, Double_t ***SigspXp, Double_t ***SigspY,
240 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***,
241 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***);
242 void CoordinateCalculation(Int_t **, Int_t **, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Int_t **,
243 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***,
244 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, vector<MC_points>&, Int_t ***);
245 void calculationOfChi2(Double_t **, Double_t **, Double_t **, Double_t **, Int_t *, Int_t *, Double_t & , Double_t *);
246 void CoordinateAlignment( Int_t **, Double_t ***, Double_t ***, Double_t ***, Int_t **, Int_t **, Double_t ***, Double_t ***);
247 void RecordingTracksAfterSpatialPoints(vector<tracksX> & , vector<tracksX> & );
248 void RecordingTracks(vector<tracksX> & , vector<tracksX> & );
249 void RecordingTracks_case3(vector<tracksX> & , vector<tracksX> & );
250 void CheckPoints(Int_t **, Int_t **, Double_t ***, Double_t ***, Double_t ***, Double_t ***,Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***,
251 Int_t **, vector<tracksX> & ,Int_t **, Int_t **,Int_t **, Double_t ***, Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***);
252 void CheckLeftover(Int_t **, Int_t **,Int_t **, Double_t ***, Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***,Double_t ***, vector<tracksX> &,
253 Int_t **, Double_t ***, Double_t ***, Double_t ***, Double_t ***);
254 void PrintAllTracks(vector<tracksX> & );
255 void CountSpatialPoints(Int_t **, Int_t *);
256 void MCefficiencyCalculation(vector<MC_points>&, vector<tracksX>&);
257 void WriteSiliconHits(Int_t **, Int_t **, Int_t **,
258 Double_t ***, Double_t ***, Double_t ***, Double_t ***,
259 Double_t ***, Double_t ***, Double_t ***, Double_t ***,
260 Double_t ***, Double_t ***, Double_t ***, Double_t ***,
261 Double_t ***, Double_t ***, Double_t ***, Double_t ***,
262 Double_t ***, Double_t ***, Double_t ***, Double_t ***, Double_t ***, Int_t ***);
263 void WriteSiliconTracks(vector<tracksX>&);
264
265 Double_t GetXYspatial_station1_2(Int_t &, Int_t &, Double_t & , Double_t &);
266 Double_t FindClusterCenter(Double_t*, Int_t , Double_t &);
267 Bool_t SkipEvent(Int_t **, Int_t **, Double_t ***, Double_t ***);
268 Int_t Angle(Int_t &, Int_t &);
269
270
271 const float tg2_5 = tan(2.5*TMath::Pi()/180.);
272 const float sin2_5 = sin(2.5*TMath::Pi()/180.);
273 const float cos2_5 = cos(2.5*TMath::Pi()/180.);
274 const Int_t kBig = 300;
275 const Int_t kMaxMC = 100;
276 const Int_t nstripXp = 640;
277 const Int_t nstripXpsm = 614;
278 const Int_t fNstrips = 643;
279 const Double_t sq12 = sqrt(12.);
280
281 const Double_t deltaX = 0.0095; //cm
282 const Double_t deltaXp = 0.01030001;//cm
283 const Double_t half_module = 6.0705;
284 const Double_t shift_after_zigzag = 0.01;//cm
285 const Double_t Xv_av = 0.4;
286 const Double_t Yv_av = 0.4;
287 const Double_t half_roadX1_X2 = 0.11;
288 const Double_t half_roadXp1_Xp2 = 0.11;
289 const Double_t half_roadX1_X3 = 0.7;
290 const Double_t half_roadX1_Xp1 = 0.5;//mc//0.3;
291 const Double_t half_roadX2_Xp2 = 0.5;//mc//0.3;
292 const Double_t half_roadX3_Xp3 = 0.55;
293 const Double_t half_roadY1_Y2 = 0.45;
294 const Double_t half_roadXp1_Xp2_diff_sign = 0.3;
295 const Double_t half_roadX1_Xp2 = .5;
296
297 const Int_t PDG_Li7 = 1000030070;//Li7
298 const Int_t PDG_He4 = 1000020040;//He4
299
300 Double_t Shift_toCenterOfMagnetX ;
301 Double_t Shift_toCenterOfMagnetY ;
302 Double_t Shift_toCenterOfMagnetAX ;
303 Double_t Shift_toCenterOfMagnetAY ;
304
305 //cut for cluster charge
306 Double_t Cut_AmplX;
307 Double_t Cut_AmplXp;
308 Double_t Cut_overflow;
309 //cut for strip charge
310 Double_t Cut_AmplStripX;
311 Double_t Cut_AmplStripXp;
312 Double_t Chi2_Globcut;
313 Int_t N_min_points;
314
315 ClassDef(BmnSiliconTrackFinder, 1)
316};
317
318#endif
sin
friend F32vec4 sin(const F32vec4 &a)
Definition P4_F32vec4.h:124
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
cos
friend F32vec4 cos(const F32vec4 &a)
Definition P4_F32vec4.h:125
BmnSiliconTrackFinder::Exec
virtual void Exec(Option_t *opt)
Definition BmnSiliconTrackFinder.cxx:31
BmnSiliconTrackFinder::compareSegments
static bool compareSegments(const tracksX &a, const tracksX &b)
Definition BmnSiliconTrackFinder.h:79
std
STL namespace.
BmnSiliconTrackFinder::MC_points
Definition BmnSiliconTrackFinder.h:83
BmnSiliconTrackFinder::MC_points::Np
Int_t Np
Definition BmnSiliconTrackFinder.h:85
BmnSiliconTrackFinder::MC_points::z
Double_t z[4]
Definition BmnSiliconTrackFinder.h:89
BmnSiliconTrackFinder::MC_points::wo3st
Bool_t wo3st
Definition BmnSiliconTrackFinder.h:90
BmnSiliconTrackFinder::MC_points::Pdg
Int_t Pdg
Definition BmnSiliconTrackFinder.h:86
BmnSiliconTrackFinder::MC_points::y
Double_t y[4]
Definition BmnSiliconTrackFinder.h:88
BmnSiliconTrackFinder::MC_points::x
Double_t x[4]
Definition BmnSiliconTrackFinder.h:87
BmnSiliconTrackFinder::MC_points::param
Double_t param[4]
Definition BmnSiliconTrackFinder.h:92
BmnSiliconTrackFinder::MC_points::Id
Int_t Id
Definition BmnSiliconTrackFinder.h:84
tracksX::SumQX
Double_t SumQX[4]
Definition BmnSiliconTrackFinder.h:62
tracksX::SigmaX
Double_t SigmaX[4]
Definition BmnSiliconTrackFinder.h:53
tracksX::CoordZ
Double_t CoordZ[4]
Definition BmnSiliconTrackFinder.h:52
tracksX::SigmaXp
Double_t SigmaXp[4]
Definition BmnSiliconTrackFinder.h:54
tracksX::CoordX
Double_t CoordX[4]
Definition BmnSiliconTrackFinder.h:49
tracksX::CoordY
Double_t CoordY[4]
Definition BmnSiliconTrackFinder.h:51
tracksX::SigmaY
Double_t SigmaY[4]
Definition BmnSiliconTrackFinder.h:55
tracksX::StripNumX
Double_t StripNumX[4]
Definition BmnSiliconTrackFinder.h:58
tracksX::param
Double_t param[4]
Definition BmnSiliconTrackFinder.h:48
tracksX::StripNumXp
Double_t StripNumXp[4]
Definition BmnSiliconTrackFinder.h:59
tracksX::ClSizeX
Double_t ClSizeX[4]
Definition BmnSiliconTrackFinder.h:60
tracksX::SumQXp
Double_t SumQXp[4]
Definition BmnSiliconTrackFinder.h:63
tracksX::CoordXp
Double_t CoordXp[4]
Definition BmnSiliconTrackFinder.h:50
tracksX::ClSizeXp
Double_t ClSizeXp[4]
Definition BmnSiliconTrackFinder.h:61