BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
drawTracks.C
Go to the documentation of this file.
1
2#include <FairHit.h>
3#include <Rtypes.h>
4#include <TChain.h>
5#include <TClonesArray.h>
6#include <TFitResultPtr.h>
7#include <TGraph.h>
8#include <TString.h>
9
10R__ADD_INCLUDE_PATH($VMCWORKDIR)
11#include "macro/run/bmnloadlibs.C"
12//
13void drawTracks(TString fileName = "srcdst_3332.root",
14 Bool_t isField = kTRUE,
15 Double_t field = 108.1,
16 Int_t nEvents = 1000)
17{
18 // void drawDchTracks(TString fileName = "$VMCWORKDIR/macro/run/srcdst_2322.root", Bool_t isField = kTRUE,
19 // Double_t field = 55.7, Int_t nEvents = 1000) { void drawDchTracks(TString fileName =
20 // "$VMCWORKDIR/macro/run/srcdst_2325.root", Bool_t isField = kTRUE, Double_t field = 75.1, Int_t nEvents = 1000)
21 // { void drawDchTracks(TString fileName = "$VMCWORKDIR/macro/run/srcdst_2327.root", Bool_t isField = kTRUE,
22 // Double_t field = 93.1, Int_t nEvents = 1000) { void drawDchTracks(TString fileName =
23 // "$VMCWORKDIR/macro/run/srcdst_2332.root", Bool_t isField = kTRUE, Double_t field = 108.1, Int_t nEvents =
24 // 1000) { void drawDchTracks(TString fileName = "$VMCWORKDIR/macro/run/srcdst_2335.root", Bool_t isField =
25 // kTRUE, Double_t field = 120.1, Int_t nEvents = 1) {
26
27 bmnloadlibs(); // load libraries
28
29 TGeoManager::Import("current_geo_file.root");
30 FairRunAna* fRunAna = new FairRunAna();
31 BmnFieldMap* magField = new BmnNewFieldMap("field_sp41v4_ascii_Extrap.root");
32 if (isField)
33 magField->SetScale(field / 55.7);
34 else
35 magField->SetScale(0.0);
36 fRunAna->SetField(magField);
37 magField->Init();
38
39 TChain* out = new TChain("bmndata");
40 out->Add(fileName.Data());
41 cout << "#recorded entries = " << out->GetEntries() << endl;
42
43 TClonesArray* dchTracks = nullptr;
44 out->SetBranchAddress("BmnDchTrack", &dchTracks);
45 TClonesArray* gemTracks = nullptr;
46 out->SetBranchAddress("BmnGlobalTrack", &gemTracks);
47 TClonesArray* silTracks = nullptr;
48 out->SetBranchAddress("BmnSiliconTrack", &silTracks);
49 TClonesArray* mwpcTracks = nullptr;
50 out->SetBranchAddress("BmnMwpcTrack", &mwpcTracks);
51 TClonesArray* mwpcSegments = nullptr;
52 out->SetBranchAddress("BmnMwpcSegment", &mwpcSegments);
53
54 Double_t z_targ = -645.0;
55 Double_t z1_sil = -431.0; //-500.0;
56 Double_t z2_sil = -318.0; //-250.0;
57 Double_t z1_dch1 = +400.0;
58 Double_t z2_dch1 = +600.0;
59 Double_t z1_dch2 = +600.0;
60 Double_t z2_dch2 = +800.0;
61 Double_t z1_gem = -208.0;
62 Double_t z2_gem = +400.0;
63 Double_t z1_mwpc = -850.0;
64 Double_t z2_mwpc = -752.0;
65 Double_t z3_mwpc = -357.0;
66 Double_t z4_mwpc = -208.0;
67
68 TCanvas* cX = new TCanvas("cX", "cX", 1000, 500);
69 cX->Range(-1100., -130., 1100., 110.);
70 TCanvas* cY = new TCanvas("cY", "cY", 1000, 500);
71 cY->Range(-1100., -130., 1100., 110.);
72
73 BmnKalmanFilter* kalman = new BmnKalmanFilter();
74
75 for (Int_t iEv = 0; iEv < nEvents; iEv++) {
76 out->GetEntry(iEv);
77
78 if (iEv % 1000 == 0)
79 printf("event: %d\n", iEv);
80
81 if (silTracks) {
82 for (Int_t iTrack = 0; iTrack < silTracks->GetEntriesFast(); iTrack++) {
83 BmnSiliconTrack* track = (BmnSiliconTrack*)silTracks->At(iTrack);
84
85 FairTrackParam* par1 = track->GetParamFirst();
86 Double_t x01 = par1->GetX();
87 Double_t y01 = par1->GetY();
88 Double_t z01 = par1->GetZ();
89 Double_t tx1 = par1->GetTx();
90 Double_t ty1 = par1->GetTy();
91
92 FairTrackParam* par = track->GetParamLast();
93 Double_t x02 = par->GetX();
94 Double_t y02 = par->GetY();
95 Double_t z02 = par->GetZ();
96 Double_t tx2 = par->GetTx();
97 Double_t ty2 = par->GetTy();
98
99 Double_t x1 = x01 + tx1 * (z_targ - z01);
100 Double_t y1 = y01 + ty1 * (z_targ - z01);
101 Double_t x2 = x02 + tx2 * (z2_sil - z02);
102 Double_t y2 = y02 + ty2 * (z2_sil - z02);
103
104 cX->cd();
105 TLine* lX = new TLine(z_targ, x1, z2_sil, x2);
106 lX->SetLineColor(kGreen + 2);
107 lX->SetLineWidth(1);
108 lX->SetLineColorAlpha(kGreen + 2, 0.05);
109 lX->Draw();
110 cY->cd();
111 TLine* lY = new TLine(z_targ, y1, z2_sil, y2);
112 lY->SetLineColor(kGreen + 2);
113 lY->SetLineWidth(1);
114 lY->SetLineColorAlpha(kGreen + 2, 0.05);
115 lY->Draw();
116 }
117 }
118
119 if (gemTracks) {
120 for (Int_t iTrack = 0; iTrack < gemTracks->GetEntriesFast(); iTrack++) {
121 BmnGlobalTrack* track = (BmnGlobalTrack*)gemTracks->At(iTrack);
122 if (track->GetNHits() < 6)
123 continue;
124
125 Int_t nSteps = (isField) ? 20 : 1;
126 Double_t step = Int_t((z2_gem - z1_gem) / nSteps);
127 FairTrackParam parPrev(*(track->GetParamFirst()));
128 // FairTrackParam* parPrev = track->GetParamLast();
129
130 Double_t xPrev = parPrev.GetX();
131 Double_t yPrev = parPrev.GetY();
132 Double_t zPrev = parPrev.GetZ();
133
134 if (isField) {
135 kalman->TGeoTrackPropagate(&parPrev, z1_gem, 2212, NULL, NULL, kTRUE);
136 xPrev = parPrev.GetX();
137 yPrev = parPrev.GetY();
138 zPrev = parPrev.GetZ();
139 } else {
140 xPrev += parPrev.GetTx() * (zPrev - z1_gem);
141 yPrev += parPrev.GetTy() * (zPrev - z1_gem);
142 zPrev = z1_gem;
143 }
144
145 for (Int_t iStep = 0; iStep < nSteps; iStep++) {
146 Double_t z = zPrev + step;
147 if (isField)
148 kalman->TGeoTrackPropagate(&parPrev, z, 2212, NULL, NULL, kTRUE);
149 Double_t x = xPrev + parPrev.GetTx() * step;
150 Double_t y = yPrev + parPrev.GetTy() * step;
151
152 cX->cd();
153 TLine* lX = new TLine(zPrev, xPrev, z, x);
154 lX->SetLineColor(kBlue);
155 lX->SetLineWidth(1);
156 lX->SetLineColorAlpha(kBlue, 0.05);
157 lX->Draw();
158 cY->cd();
159 TLine* lY = new TLine(zPrev, yPrev, z, y);
160 lY->SetLineColor(kBlue);
161 lY->SetLineWidth(1);
162 lY->SetLineColorAlpha(kBlue, 0.05);
163 lY->Draw();
164
165 xPrev = x;
166 yPrev = y;
167 zPrev = z;
168 }
169 }
170 }
171
172 if (mwpcTracks) {
173 for (Int_t iTr = 0; iTr < mwpcTracks->GetEntriesFast(); ++iTr) {
174 BmnTrack* track = (BmnTrack*)mwpcTracks->At(iTr);
175
176 Int_t ip = -1;
177 Double_t Z = track->GetParamFirst()->GetZ();
178 Int_t itr = track->GetFlag();
179
180 if (Z >= -809. && Z <= -807.)
181 ip = 0;
182 if (Z >= -290. && Z <= -288.)
183 ip = 1;
184
185 if (ip < 0)
186 continue;
187
188 Int_t Nhits = track->GetNHits();
189
190 Double_t x0 = track->GetParamFirst()->GetX();
191 Double_t y0 = track->GetParamFirst()->GetY();
192 Double_t z0 = track->GetParamFirst()->GetZ();
193 Double_t tx = track->GetParamFirst()->GetTx();
194 Double_t ty = track->GetParamFirst()->GetTy();
195
196 if (ip == 0) {
197 Double_t x1 = x0 + tx * (z1_mwpc - z0);
198 Double_t y1 = y0 + ty * (z1_mwpc - z0);
199 Double_t x2 = x0 + tx * (z_targ - z0);
200 Double_t y2 = y0 + ty * (z_targ - z0);
201
202 cX->cd();
203 TLine* lX = new TLine(z1_mwpc, x1, z_targ, x2);
204 lX->SetLineColor(kMagenta + 2);
205 lX->SetLineWidth(1);
206 lX->SetLineColorAlpha(kMagenta + 2, 0.05);
207 lX->Draw();
208 cY->cd();
209 TLine* lY = new TLine(z1_mwpc, y1, z_targ, y2);
210 lY->SetLineColor(kMagenta + 2);
211 lY->SetLineWidth(1);
212 lY->SetLineColorAlpha(kMagenta + 2, 0.05);
213 lY->Draw();
214 }
215 if (ip == 1) {
216 Double_t x1 = x0 + tx * (z3_mwpc - z0);
217 Double_t y1 = y0 + ty * (z3_mwpc - z0);
218 Double_t x2 = x0 + tx * (z4_mwpc - z0);
219 Double_t y2 = y0 + ty * (z4_mwpc - z0);
220
221 cX->cd();
222 TLine* lX = new TLine(z3_mwpc, x1, z4_mwpc, x2);
223 lX->SetLineColor(kMagenta + 2);
224 lX->SetLineWidth(1);
225 lX->SetLineColorAlpha(kMagenta + 2, 0.05);
226 lX->Draw();
227 cY->cd();
228 TLine* lY = new TLine(z3_mwpc, y1, z4_mwpc, y2);
229 lY->SetLineColor(kMagenta + 2);
230 lY->SetLineWidth(1);
231 lY->SetLineColorAlpha(kMagenta + 2, 0.05);
232 lY->Draw();
233 }
234 } // iTr
235 }
236
237 if (dchTracks) {
238 for (Int_t iTrack = 0; iTrack < dchTracks->GetEntriesFast(); iTrack++) {
239 BmnDchTrack* track = (BmnDchTrack*)dchTracks->At(iTrack);
240 if (track->GetNHits() < 8)
241 continue;
242 // if (track->GetFlag() != 7) continue;
243
244 FairTrackParam* par = track->GetParamFirst();
245 Double_t z0 = par->GetZ();
246 if (z0 < 550) {
247 Double_t x0 = par->GetX();
248 Double_t y0 = par->GetY();
249 Double_t tx = par->GetTx();
250 Double_t ty = par->GetTy();
251
252 Double_t x1 = x0 + tx * (z1_dch1 - z0);
253 Double_t y1 = y0 + ty * (z1_dch1 - z0);
254 Double_t x2 = x0 + tx * (z2_dch1 - z0);
255 Double_t y2 = y0 + ty * (z2_dch1 - z0);
256
257 cX->cd();
258 TLine* lX = new TLine(z1_dch1, x1, z2_dch1, x2);
259 lX->SetLineColor(kRed);
260 lX->SetLineWidth(1);
261 lX->SetLineColorAlpha(kRed, 0.05);
262 lX->Draw();
263 cY->cd();
264 TLine* lY = new TLine(z1_dch1, y1, z2_dch1, y2);
265 lY->SetLineColor(kRed);
266 lY->SetLineWidth(1);
267 lY->SetLineColorAlpha(kRed, 0.05);
268 lY->Draw();
269 } else if (z0 > 650) {
270 Double_t x0 = par->GetX();
271 Double_t y0 = par->GetY();
272 Double_t tx = par->GetTx();
273 Double_t ty = par->GetTy();
274
275 Double_t x1 = x0 + tx * (z1_dch2 - z0);
276 Double_t y1 = y0 + ty * (z1_dch2 - z0);
277 Double_t x2 = x0 + tx * (z2_dch2 - z0);
278 Double_t y2 = y0 + ty * (z2_dch2 - z0);
279
280 cX->cd();
281 TLine* lX = new TLine(z1_dch2, x1, z2_dch2, x2);
282 lX->SetLineColor(kRed);
283 lX->SetLineWidth(1);
284 lX->SetLineColorAlpha(kRed, 0.05);
285 lX->Draw();
286 cY->cd();
287 TLine* lY = new TLine(z1_dch2, y1, z2_dch2, y2);
288 lY->SetLineColor(kRed);
289 lY->SetLineWidth(1);
290 lY->SetLineColorAlpha(kRed, 0.05);
291 lY->Draw();
292 }
293 }
294 }
295 }
296
297 delete kalman;
298
299 cX->cd();
300 TGaxis* axisX1 = new TGaxis(-900, -100, -900, 100, -100, 100, 510, "");
301 axisX1->SetTitle("X [cm]");
302 axisX1->SetLineWidth(1);
303 axisX1->SetLineColor(kBlack);
304 axisX1->SetLineColorAlpha(kBlack, 0.5);
305 axisX1->Draw();
306
307 TGaxis* axisY1 = new TGaxis(-900, -100, 1000, -100, -900, 1000, 510, "");
308 axisY1->SetTitle("Z [cm]");
309 axisY1->SetLineWidth(1);
310 axisY1->SetLineColor(kBlack);
311 axisY1->SetLineColorAlpha(kBlack, 0.5);
312 axisY1->Draw();
313
314 const Int_t nPlane = 16; // 11;
315
316 // PC0 PC1 t Si1/2 PC3 Si3 PC4
317 Double_t zPos[nPlane] = {-850, -752, -645, -433, -429, -357, -318, -207.4, 65, 89, 113, 137, 161, 185, 510, 710};
318 Double_t xPos[nPlane] = {14, 14, 3, 6, 6, 14, 12, 14, 80, 80, 80, 80, 80, 80, 90, 90};
319 Double_t yPos[nPlane] = {14, 14, 3, 6, 6, 14, 12, 14, 23, 23, 23, 23, 23, 23, 90, 90};
320 for (Int_t i = 0; i < nPlane; ++i) {
321 TLine* lst = new TLine(zPos[i], -xPos[i], zPos[i], xPos[i]);
322 lst->SetLineColor(kBlack);
323 lst->SetLineColorAlpha(kBlack, 0.5);
324 lst->Draw();
325 }
326
327 cY->cd();
328 TGaxis* axisX2 = new TGaxis(-870, -100, -870, 100, -100, 100, 510, "");
329 axisX2->SetTitle("Y [cm]");
330 axisX2->SetLineColor(kBlack);
331 axisX2->SetLineColorAlpha(kBlack, 0.5);
332 axisX2->Draw();
333
334 TGaxis* axisY2 = new TGaxis(-900, -100, 850, -100, -900, 900, 510, "");
335 axisY2->SetTitle("Z [cm]");
336 axisY2->SetLineWidth(1);
337 axisY2->SetLineColor(kBlack);
338 axisY2->SetLineColorAlpha(kBlack, 0.5);
339 axisY2->Draw();
340
341 for (Int_t i = 0; i < nPlane; ++i) {
342 TLine* lst = new TLine(zPos[i], -yPos[i], zPos[i], yPos[i]);
343 lst->SetLineWidth(1);
344 lst->SetLineColor(kBlack);
345 lst->SetLineColorAlpha(kBlack, 0.5);
346 lst->Draw();
347 }
348
349 cX->SaveAs("tracks_ZX.pdf");
350 cY->SaveAs("tracks_ZY.pdf");
351}
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
R__ADD_INCLUDE_PATH
#define R__ADD_INCLUDE_PATH(PATH)
kRed
kRed
kBlack
kBlack
kGreen
kGreen
kMagenta
kMagenta
kBlue
kBlue
BmnFieldMap::SetScale
void SetScale(Double_t factor)
Set a global field scaling factor.
Definition BmnFieldMap.h:55
BmnFieldMap::Init
virtual void Init()
Definition BmnFieldMap.cxx:120
BmnGlobalTrack
BmnGlobalTrack.h.
Definition BmnGlobalTrack.h:14
BmnKalmanFilter::TGeoTrackPropagate
BmnStatus TGeoTrackPropagate(FairTrackParam *par, Double_t zOut, Int_t pdg, vector< Double_t > *F, Double_t *length, Bool_t isField)
Definition BmnKalmanFilter.cxx:438
BmnNewFieldMap
Field Map does not have any symmetries w.r.t coordinate axises.
Definition BmnNewFieldMap.h:16
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:63
BmnTrack::GetNHits
Int_t GetNHits() const
Accessors.
Definition BmnTrack.h:49
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:65
BmnTrack::GetFlag
Int_t GetFlag() const
Definition BmnTrack.h:53
TAttLine::SetLineWidth
virtual void SetLineWidth(Width_t lwidth)
TAttLine::SetLineColorAlpha
virtual void SetLineColorAlpha(Color_t lcolor, Float_t lalpha)
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
TChain::SetBranchAddress
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
TChain::Add
virtual Int_t Add(const char *name, Long64_t nentries=TTree::kMaxEntries)
TChain::GetEntries
Long64_t GetEntries() const override
TChain::GetEntry
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
TGaxis::SetTitle
virtual void SetTitle(const char *title="")
TGeoManager::Import
static TGeoManager * Import(const char *filename, const char *name="", Option_t *option="")
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::At
TObject * At(Int_t idx) const override
TObject::Draw
virtual void Draw(Option_t *option="")
TPad::Range
void Range(Double_t x1, Double_t y1, Double_t x2, Double_t y2) override
TPad::SaveAs
void SaveAs(const char *filename="", Option_t *option="") const override
drawTracks
void drawTracks(TString fileName="srcdst_3332.root", Bool_t isField=kTRUE, Double_t field=108.1, Int_t nEvents=1000)
Definition drawTracks.C:13