Fork me on GitHub

source: git/external/TrackCovariance/AcceptanceClx.cc@ cd7aa16

Last change on this file since cd7aa16 was f17e10d, checked in by michele <michele.selvaggi@…>, 4 years ago

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