Fork me on GitHub

source: git/external/TrackCovariance/AcceptanceClx.cc@ 0983d48

Last change on this file since 0983d48 was ebf40fd, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Major update to handle highly displaced tracks
  • Property mode set to 100644
File size: 9.8 KB
RevLine 
[ebf40fd]1#include <TVector3.h>
2#include "AcceptanceClx.h"
3//
4// Pt splitting routine
5//
6void AcceptanceClx::VecInsert(Int_t i, Float_t x, TVectorF& Vec)
7{
8 // Insert a new element x in location i+1 of vector Vec
9 //
10 Int_t N = Vec.GetNrows(); // Get vector nitial size
11 Vec.ResizeTo(N + 1); // Increase size
12 for (Int_t j = N - 1; j > i; j--)Vec(j + 1) = Vec(j); // Shift all elements above i
13 Vec(i+1) = x;
14}
15//
16void AcceptanceClx::SplitPt(Int_t i, TVectorF &AccPt)
17{
18 Int_t Nrows = fAcc.GetNrows();
19 Int_t Ncols = fAcc.GetNcols();
20 TMatrixF AccMod(Nrows + 1, Ncols); // Size output matrix
21 for (Int_t ith = 0; ith < Ncols; ith++)
22 {
23 AccMod(i + 1, ith) = AccPt(ith);
24 for (Int_t ipt = 0; ipt <= i; ipt++) AccMod(ipt, ith) = fAcc(ipt, ith);
25 for (Int_t ipt = i + 1; ipt < Nrows; ipt++) AccMod(ipt + 1, ith) = fAcc(ipt, ith);
26 }
27 //
28 fAcc.ResizeTo(Nrows + 1, Ncols);
29 fAcc = AccMod;
30}
31//
32// Theta splitting routine
33void AcceptanceClx::SplitTh(Int_t i, TVectorF &AccTh)
34{
35 Int_t Nrows = fAcc.GetNrows();
36 Int_t Ncols = fAcc.GetNcols();
37 TMatrixF AccMod(Nrows, Ncols + 1); // Size output matrix
38 for (Int_t ipt = 0; ipt < Nrows; ipt++)
39 {
40 AccMod(ipt, i + 1) = AccTh(ipt);
41 for (Int_t ith = 0; ith <= i; ith++) AccMod(ipt, ith) = fAcc(ipt, ith);
42 for (Int_t ith = i + 1; ith < Ncols; ith++) AccMod(ipt, ith + 1) = fAcc(ipt, ith);
43 }
44 //
45 fAcc.ResizeTo(Nrows, Ncols + 1);
46 fAcc = AccMod;
47}
48//
49// Constructors
50//
51AcceptanceClx::AcceptanceClx(TString InFile)
52{
53 ReadAcceptance(InFile);
54}
55//
56AcceptanceClx::AcceptanceClx(SolGeom* InGeo)
57{
58 // Initializations
59 //
60 //cout << "Entered constructor of AccpeptanceClx" << endl;
61 // Setup grid
62 // Start grid parameters
63 //
64 // Pt nodes
65 const Int_t NpPtInp = 5;
66 Float_t PtInit[NpPtInp] = { 0., 1., 10., 100., 250. };
67 TVectorF Pta(NpPtInp, PtInit);
68 Int_t NpPt = Pta.GetNrows(); // Nr. of starting pt points
69 // Theta nodes
70 const Int_t NpThInp = 15;
71 Float_t ThInit[NpThInp] = { 0.,5.,10.,20.,30.,40.,50.,90.,
72 130.,140.,150.,160.,170.,175., 180. };
73 TVectorF Tha(NpThInp, ThInit);
74 Int_t NpTh = Tha.GetNrows(); // Nr. of starting theta points
75 //cout << "AcceptanceClv:: Pta and Tha arrays defined" << endl;
76 //
77 // Grid splitting parameters
78 Float_t dPtMin = 0.2; // Grid Pt resolution (GeV)
79 Float_t dThMin = 2.0; // Grid Theta resolution (degrees)
80 Float_t dAmin = 1.0; // Minimum # hits step
81 //
82 fAcc.ResizeTo(NpPt, NpTh);
83 //
84 //
85 // Event loop: fill matrix starting nodes
86 //
87 TVector3 xv(0., 0., 0.);
88 for (Int_t ipt = 0; ipt < NpPt; ipt++) // Scan pt bins
89 {
90 // Momentum in GeV
91 Double_t pt = Pta(ipt);
92 //
93 for (Int_t ith = 0; ith < NpTh; ith++) // Scan theta bins
94 {
95 //
96 // Theta in from degrees to radians
97 Double_t th = TMath::Pi() * Tha(ith) / 180.;
98 Double_t pz = pt / TMath::Tan(th);
99 TVector3 tp(pt, 0., pz);
100 //
101 // Get number of measurement hits
102 //
103 SolTrack* gTrk = new SolTrack(xv, tp, InGeo); // Generated track
104 Int_t Mhits = gTrk->nmHit(); // Nr. Measurement hits
105 fAcc(ipt, ith) = (Float_t)Mhits;
106 }
107 }
108 //
109 // Scan nodes and split if needed
110 //
111 Int_t Nsplits = 1; // Number of split per iteration
112 Int_t Ncycles = 0; // Number of iterations
113 Int_t MaxSplits = 200; // Maximum number of splits
114 Int_t NsplitCnt = 0;
115 //
116 while (Nsplits > 0)
117 {
118 Nsplits = 0;
119 // Scan nodes
120 for (Int_t ipt = 0; ipt < NpPt - 1; ipt++) // Scan pt bins
121 {
122 for (Int_t ith = 0; ith < NpTh - 1; ith++) // Scan theta bins
123 {
124 Float_t dAp = TMath::Abs(fAcc(ipt + 1, ith) - fAcc(ipt, ith));
125 Float_t dPt = TMath::Abs(Pta(ipt + 1) - Pta(ipt));
126 //
127 // Pt split
128 if (dPt > dPtMin && dAp > dAmin && Nsplits < MaxSplits) {
129 NsplitCnt++; // Total splits counter
130 Nsplits++; // Increase splits/cycle
131 NpPt++; // Increase #pt points
132 Float_t newPt = 0.5 * (Pta(ipt + 1) + Pta(ipt));
133 VecInsert(ipt, newPt, Pta);
134 TVectorF AccPt(NpTh);
135 for (Int_t i = 0; i < NpTh; i++)
136 {
137 Double_t pt = newPt;
138 Double_t th = TMath::Pi() * Tha[i] / 180.;
139 Double_t pz = pt / TMath::Tan(th);
140 TVector3 tp(pt, 0., pz);
141 SolTrack* gTrk = new SolTrack(xv, tp, InGeo); // Generated track
142 Int_t Mhits = gTrk->nmHit(); // Nr. Measurement hits
143 AccPt(i) = (Float_t)Mhits;
144 }
145 SplitPt(ipt, AccPt);
146 // Completed Pt split
147 }
148 //
149 Float_t dAt = TMath::Abs(fAcc(ipt, ith + 1) - fAcc(ipt, ith));
150 Float_t dTh = TMath::Abs(Tha(ith + 1) - Tha(ith));
151 //
152 // Theta split
153 if (dTh > dThMin && dAt > dAmin && Nsplits < MaxSplits) {
154 //cout << "Th(" << ith << ") = " << Tha(ith) << ", dAt = " << dAt << endl;
155 NsplitCnt++; // Total splits counter
156 Nsplits++; // Increase splits
157 NpTh++; // Increase #pt points
158 Float_t newTh = 0.5 * (Tha(ith + 1) + Tha(ith));
159 VecInsert(ith, newTh, Tha);
160 TVectorF AccTh(NpPt);
161 for (Int_t i = 0; i < NpPt; i++)
162 {
163 Double_t pt = Pta(i);
164 Double_t th = TMath::Pi() * newTh / 180.;
165 Double_t pz = pt / TMath::Tan(th);
166 TVector3 tp(pt, 0., pz);
167 SolTrack* gTrk = new SolTrack(xv, tp, InGeo); // Generated track
168 Int_t Mhits = gTrk->nmHit(); // Nr. Measurement hits
169 AccTh(i) = (Float_t)Mhits;
170 }
171 SplitTh(ith, AccTh);
172 // Theta splits completed
173 }
174 } // End loop on theta nodes
175 } // End loop on pt nodes
176 Ncycles++;
177 } // End loop on iteration
178 //
179 std::cout<<"AcceptanceClx:: Acceptance generation completed after "<<Ncycles<<" cycles"
180 <<" and a total number of "<<NsplitCnt<< " splits"<<std::endl;
181 //
182 // Store final variables and parameters
183 //
184 fNPtNodes = NpPt;
185 Int_t PtCk = Pta.GetNrows();
186 if (PtCk != NpPt)
187 std::cout << "AcceptanceClx:: Error in grid generation NpPt=" << NpPt << ", Pta size= " << PtCk << std::endl;
188 fPtArray.ResizeTo(NpPt);
189 fPtArray = Pta; // Array of Pt nodes
190 //
191 fNThNodes = NpTh;
192 Int_t ThCk = Tha.GetNrows();
193 if (ThCk != NpTh)
194 std::cout << "AcceptanceClx:: Error in grid generation NpTh=" << NpTh << ", Tha size= " << ThCk << std::endl;
195 fThArray.ResizeTo(NpTh);
196 fThArray = Tha; // Array of Theta nodes
197 //
198 std::cout << "AcceptanceClx:: Acceptance encoding with " << fNPtNodes
199 <<" pt nodes and "<< fNThNodes <<" theta nodes"<< std::endl;
200 Int_t Nrows = fAcc.GetNrows();
201 Int_t Ncols = fAcc.GetNcols();
202}
203
204// Destructor
205AcceptanceClx::~AcceptanceClx()
206{
207 fNPtNodes = 0;
208 fNThNodes = 0;
209 fAcc.Clear();
210 fPtArray.Clear();
211 fThArray.Clear();
212}
213//
214void AcceptanceClx::WriteAcceptance(TFile *fout)
215{
216 //
217 // Write out data
218 TTree* tree = new TTree("treeAcc", "Acceptance tree");
219 TMatrixF* pntMatF = &fAcc;
220 TVectorF* pntVecP = &fPtArray;
221 TVectorF* pntVecT = &fThArray;
222 tree->Branch("AcceptanceMatrix", "TMatrixF", &pntMatF, 64000, 0);
223 tree->Branch("AcceptancePtVec", "TVectorF", &pntVecP, 64000, 0);
224 tree->Branch("AcceptanceThVec", "TVectorF", &pntVecT, 64000, 0);
225 tree->Fill();
226 fout->Write();
227}
228//
229void AcceptanceClx::WriteAcceptance(TString OutFile)
230{
231 //
232 // Write out data
233 TFile* fout = new TFile(OutFile,"RECREATE");
234 WriteAcceptance(fout);
235 fout->Close();
236 delete fout;
237}
238//
239void AcceptanceClx::ReadAcceptance(TString InFile)
240{
241 //
242 // Read in data
243 TFile* f = new TFile(InFile, "READ");
244 //
245 // Import TTree
246 TTree* T = (TTree*)f->Get("treeAcc");
247 //
248 // Get matrix
249 TMatrixF* pAcc = new TMatrixF();
250 T->SetBranchAddress("AcceptanceMatrix", &pAcc);
251 T->GetEntry(0);
252 //
253 fNPtNodes = pAcc->GetNrows();
254 fNThNodes = pAcc->GetNcols();
255 fAcc.ResizeTo(fNPtNodes, fNThNodes);
256 fAcc = *pAcc;
257 //
258 // Get Pt grid
259 TVectorF* pPtArray = new TVectorF();
260 T->SetBranchAddress("AcceptancePtVec", &pPtArray);
261 T->GetEntry(0);
262 fPtArray.ResizeTo(fNPtNodes);
263 fPtArray = *pPtArray;
264 //
265 // Get Theta array
266 TVectorF* pThArray = new TVectorF();
267 T->SetBranchAddress("AcceptanceThVec", &pThArray);
268 T->GetEntry(0);
269 fThArray.ResizeTo(fNThNodes);
270 fThArray = *pThArray;
271 //
272 std::cout << "AcceptanceClx::Read complete: Npt= " << fNPtNodes << ", Nth= " << fNThNodes << std::endl;
273 //
274 f->Close();
275 delete f;
276}
277//
278Double_t AcceptanceClx::HitNumber(Double_t pt, Double_t theta)
279{
280 //
281 // Protect against values out of range
282 Float_t eps = 1.0e-4;
283 Float_t pt0 = (Float_t)pt;
284 if (pt0 <= fPtArray(0)) pt0 = fPtArray(0) + eps;
285 else if (pt0 >= fPtArray(fNPtNodes-1)) pt0 = fPtArray(fNPtNodes-1) - eps;
286 Float_t th0 = (Float_t)theta;
287 if (th0 <= fThArray(0)) th0 = fThArray(0) + eps;
288 else if (th0 >= fThArray(fNThNodes - 1)) th0 = fThArray(fNThNodes - 1) - eps;
289 //
290 // Find cell
291 Float_t* parray = fPtArray.GetMatrixArray();
292 Int_t ip = TMath::BinarySearch(fNPtNodes, parray, pt0);
293 Float_t* tarray = fThArray.GetMatrixArray();
294 Int_t it = TMath::BinarySearch(fNThNodes, tarray, th0);
295 //
296 if (ip<0 || ip >fNPtNodes - 2)
297 {
298 std::cout << "Search error: (ip, pt) = (" << ip << ", " << pt << "), pt0 = " << pt0 << std::endl;
299 std::cout << "Search error: pt nodes = " << fNPtNodes
300 << " , last value = " << fPtArray(fNPtNodes - 1) << std::endl;
301 }
302 if (it<0 || ip >fNThNodes - 2)
303 {
304 std::cout << "Search error: (it, th) = (" << it << ", " << theta << "), th0 = " << th0 << std::endl;
305 std::cout << "Search error: th nodes = " << fNThNodes
306 << " , last value = " << fThArray(fNThNodes - 1) << std::endl;
307 }
308 //
309 // Bilinear interpolation
310 //
311 Double_t spt = (pt0 - fPtArray(ip)) / (fPtArray(ip + 1) - fPtArray(ip));
312 Double_t sth = (th0 - fThArray(it)) / (fThArray(it + 1) - fThArray(it));
313 Double_t A11 = fAcc(ip, it);
314 Double_t A12 = fAcc(ip, it+1);
315 Double_t A21 = fAcc(ip+1, it);
316 Double_t A22 = fAcc(ip+1, it+1);
317 Double_t A = A11 * (1 - spt) * (1 - sth) + A12 * (1 - spt) * sth +
318 A21 * spt * (1 - sth) + A22 * spt * sth;
319 //
320 return A;
321}
322//
323Double_t AcceptanceClx::HitNum(Double_t *x, Double_t *p) // Theta in degrees
324{
325 Double_t pt = x[0];
326 Double_t th = x[1];
327 //
328 return HitNumber(pt, th);
329}
330//
331
Note: See TracBrowser for help on using the repository browser.