BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnGemFastDigitize.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- BmnGemFastDigitize source file -----
3// ----- Created 4/07/19 by A. Zinchenko -----
4// -------------------------------------------------------------------------
5
6#include "BmnGemFastDigitize.h"
7
8#include "CbmMCTrack.h"
9#include "CbmStsPoint.h"
10#include "CbmStsSensor.h"
11#include "FairRootManager.h"
12#include "TFile.h"
13
14#include <TCut.h>
15#include <TDatabasePDG.h>
16#include <TEntryList.h>
17#include <TGeoManager.h>
18#include <TH1I.h>
19#include <TH2I.h>
20#include <TParticlePDG.h>
21#include <TROOT.h>
22#include <TRandom.h>
23#include <TSystem.h>
24#include <TTree.h>
25// #include <TH1F.h>
26// #include <TH3F.h>
27
28#include <fstream>
29#include <iostream>
30#include <queue>
31
32using std::cout;
33using std::endl;
34using std::ifstream;
35using std::iterator;
36using std::map;
37using std::multimap;
38using std::pair;
39using std::set;
40using std::vector;
41
42// ----- Default constructor -------------------------------------------
43BmnGemFastDigitize::BmnGemFastDigitize()
44 : FairTask("GemFastDigitize")
45 , fhInd(NULL)
46{}
47// -------------------------------------------------------------------------
48
49// ----- Destructor ----------------------------------------------------
50BmnGemFastDigitize::~BmnGemFastDigitize() {}
51// -------------------------------------------------------------------------
52
53// ----- Public method Init --------------------------------------------
54InitStatus BmnGemFastDigitize::Init()
55{
56
57 // Get RootManager
58 FairRootManager* ioman = FairRootManager::Instance();
59 if (!ioman) {
60 cout << "-E- BmnGemFastDigitize::Init: "
61 << "RootManager not instantiated!" << endl;
62 return kFATAL;
63 }
64
65 // Get input arrays
66 fPoints = (TClonesArray*)ioman->GetObject("StsPoint");
67 if (!fPoints) {
68 // cout << "-W- BmnGetFastDigitize::Init: "<< "No STSPoint array!" << endl;
69 return kERROR;
70 }
71 fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack");
72
73 // Get list of files with digit trees
74 if (gSystem->Exec("du -s clusLib")) {
75 LOG(fatal) << " !!! No clusLib directory found - exit !!!";
76 }
77 // gSystem->Exec("ls -1 clusLib/thx*histo.root > aaa");
78 // gSystem->Exec("ls -1 thxp10thyp00.histo.root > aaa");
79
80 // ifstream fin("aaa");
81 ifstream fin("clusLib/file_list.txt");
82 TString fileName;
83 // const Int_t nbinsx = 15, nbinsy = 5;
84 // Double_t xlim[nbinsx] = {-75, 75}; // theta_x (deg)
85 const Int_t nbinsx = 29, nbinsy = 5;
86 Double_t xlim[nbinsx] = {-72.5, 72.5}; // theta_x (deg)
87 Double_t ylim[nbinsy] = {-25, 25}; // theta_y
88 fhInd = new TH2I("hInd", "", nbinsx, xlim[0], xlim[1], nbinsy, ylim[0], ylim[1]);
89 fhInd->StatOverflows(kTRUE); // correct numbering of underflow / overflow bins
90
91 while (1) {
92 fin >> fileName;
93 if (fin.eof())
94 break;
95
96 // Decode file name to get theta_X and theta_Y
97 // Int_t ithx = TString(fileName(4,2)).Atoi(), ithy = TString(fileName(10,2)).Atoi();
98 // if (fileName(4) == 'm') ithx = -ithx;
99 // if (fileName(10) == 'm') ithy = -ithy;
100 Int_t ithx = TString(fileName(12, 2)).Atoi(), ithy = TString(fileName(18, 2)).Atoi();
101 if (fileName(11) == 'm')
102 ithx = -ithx;
103 if (fileName(17) == 'm')
104 ithy = -ithy;
105 TString hname = "h";
106 hname += ithx;
107 hname += "_";
108 hname += ithy;
109 /*
110 //Double_t xmin[2] = {238.0, 252.5}, xmax[2] = {242.5, 257.0};
111 Double_t xmin[2] = {236.8, 251.4}, xmax[2] = {238.3, 252.8};
112 Int_t nbins = 100;
113 //TH2I *hXY = new TH2I("hXY","Y vs X",nbins,xmin[0],xmax[0],nbins,xmin[1],xmax[1]); // Y vs X (cog[1] vs cog[0])
114 //clus->Project(hXY->GetName(),"cog[1]:cog[0]",cut);
115 */
116 // TFile *file = TFile::Open("checkDigi_st5_new.histo.root");
117 TFile* file = TFile::Open(fileName);
118 TTree* clus = (TTree*)file->Get("clus");
119
120 // Get stored histogram
121 TH2I* hXY = (TH2I*)file->Get("hXY");
122 hXY = (TH2I*)hXY->Clone(hname);
123 Int_t indXY = fhInd->Fill(ithx, ithy);
124 fHistMap[indXY] = hXY;
125 fTreeMap[indXY] = clus;
126
127 clus->SetBranchStatus("*", 0); // disable branches
128 clus->SetBranchStatus("qtot*", 1); // enable branches
129 clus->SetBranchStatus("ibeg*", 1); // enable branches
130 clus->SetBranchStatus("strin*", 1); // enable branches
131 clus->SetBranchStatus("dpitin*", 1); // enable branches
132 clus->SetBranchStatus("strout*", 1); // enable branches
133 clus->SetBranchStatus("cog*", 1); // enable branches
134 clus->SetBranchStatus("digis*", 1); // enable branches
135
136 // TList *list = file.GetList();
137 // Int_t ntot = list->GetEntries();
138
139 // TCut cut0("abs(imax[0]-240)<5&&strout[0]==244"); // side 0
140 // TCut cut1("abs(imax[1]-254)<5&&strout[1]==259&&imax[1]>252"); // side 1
141 TCut cut0("abs(ibeg[0]-237)<5&&strout[0]==244"); // side 0
142 TCut cut1("abs(ibeg[1]-252)<5&&strout[1]==259"); // side 1
143 TCut cut = ""; // cut0 && cut1;
144 if (file->Get("CUT"))
145 cut = *((TCut*)file->Get("CUT"));
146
147 TEntryList* entryList = NULL;
148 clus->Draw(">>entryList", cut, "entrylist", 9000000000);
149 entryList = (TEntryList*)gDirectory->Get("entryList");
150 clus->SetEntryList(entryList);
151 cout << " xxxxxx " << fileName << " " << entryList->GetN() << " " << ithx << " " << ithy << " " << indXY
152 << endl;
153
154 fvDig0 = &fvDigis[0];
155 fvDig1 = &fvDigis[1];
156
157 clus->SetBranchAddress("qtot[2]", &fQtot);
158 clus->SetBranchAddress("ibeg[2]", &fIbeg);
159 clus->SetBranchAddress("strin[2]", &fStrin);
160 clus->SetBranchAddress("dpitin[2]", &fDpitin);
161 clus->SetBranchAddress("strout[2]", &fStrout);
162 clus->SetBranchAddress("cog[2]", &fCog);
163 clus->SetBranchAddress("digis0", &fvDig0);
164 clus->SetBranchAddress("digis1", &fvDig1);
165
166 Int_t entries = entryList->GetN();
167 entries = TMath::Min(10000, entries);
168 map<Int_t, Int_t> strin[2]; // entrance strips
169 std::multimap<Float_t, Int_t> mmm;
170 fXYmap[0][indXY] = mmm;
171 fXYmap[1][indXY] = mmm;
172
173 for (Int_t j = 0; j < entries; ++j) {
174 Int_t treeEn = entryList->GetEntry(j);
175 clus->GetEntry(treeEn);
176 // cout << j << " " << entryList->GetEntry(j) << " " << cog[0] << " " << cog[1] << endl;
177
178 for (Int_t side = 0; side < 2; ++side) {
179 fXYmap[side][indXY].insert(pair<Float_t, Int_t>(fCog[side], treeEn));
180 if (strin[side].find(fStrin[side]) == strin[side].end())
181 strin[side][fStrin[side]] = 0;
182 ++strin[side][fStrin[side]];
183 }
184 }
185
186 // Find reference strip
187 for (Int_t side = 0; side < 2; ++side) {
188 map<Int_t, Int_t>::iterator it = strin[side].begin();
189 fStrRef[side][indXY] = it->first;
190 Int_t maxStrip = it->second;
191
192 for (; it != strin[side].end(); ++it) {
193 if (it->second > maxStrip) {
194 maxStrip = it->second;
195 fStrRef[side][indXY] = it->first;
196 }
197 }
198 }
199
200 // file.Close();
201 }
202
203 cout << "-I- BmnGemFastDigitize: Initialisation successful" << endl;
204 return kSUCCESS;
205}
206// -------------------------------------------------------------------------
207
208// ----- Public method Exec --------------------------------------------
209void BmnGemFastDigitize::Exec(Option_t* opt) {}
210// -------------------------------------------------------------------------
211
212// ----- Public method Finish ------------------------------------------
213void BmnGemFastDigitize::Finish() {}
214// -------------------------------------------------------------------------
215
216// ----- Public method ProduceHitResponse ------------------------------
217void BmnGemFastDigitize::ProduceHitResponseFast(CbmStsSensor* sensor,
218 std::set<Int_t>& pSet,
219 Double_t* stripSignalF,
220 Double_t* stripSignalB,
221 map<Int_t, set<Int_t>>& chanPointMapF,
222 map<Int_t, set<Int_t>>& chanPointMapB)
223{
224 // Produce sensor response
225
226 Int_t iPoint = -1, stationNr = Int_t(sensor->GetStationNr());
227 CbmStsPoint* point = NULL;
228 set<Int_t>::iterator it1;
229 Double_t dPitch = 0;
230
231 for (it1 = pSet.begin(); it1 != pSet.end(); it1++) {
232
233 // AZ Apply overall efficiency
234 // const Double_t eff0 = 0.65;
235 // if (gRandom->Rndm() > eff0) continue;
236
237 iPoint = (*it1);
238 point = (CbmStsPoint*)fPoints->UncheckedAt(iPoint);
239 TVector3 mom;
240 point->Momentum(mom);
241
242 CbmMCTrack* mcTr = (CbmMCTrack*)fMCTracks->UncheckedAt(point->GetTrackID()); // ES
243 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(mcTr->GetPdgCode());
244 Double_t mnoc = 1.0, mass = 3.7283999; // He4
245 if (part)
246 mass = part->Mass();
247 if (mass > 0.0001) {
248 TLorentzVector lorv;
249 lorv.SetVectM(mom, mass);
250 Double_t beta = lorv.Beta();
251 Double_t gamma = lorv.Gamma();
252 Double_t charge = (part) ? part->Charge() / 3 : 1;
253 mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
254 if (mnoc < 0.5)
255 continue; // AZ
256 // cout << " ***clus*** " << stationNr << " " << mass << " " << beta << " " << gamma << " " << mnoc << endl;
257 }
258 Double_t coeff = mnoc / 29.2; // 29.2 - mnoc for test particles (pion with p = 1.6 GeV/c)
259 // AZ 31.07.20 coeff *= 2; // to correct ADC conversion factor
260 coeff *= 1.; // to correct ADC conversion factor
261
262 // Double_t moduleThickness = 0.3 + 0.25 + 0.2 + 0.15 ; // DriftGapThickness + FirstTransferGapThickness +
263 // SecondTransferGapThickness + InductionGapThickness;
264
265 TVector3 posIn, posOut;
266 point->PositionIn(posIn);
267 point->PositionOut(posOut);
268 TVector3 dpos = posOut;
269 dpos -= posIn;
270 dpos *= 0.001;
271 posIn += dpos;
272 posOut -= dpos;
273
274 Double_t xrand, yrand;
275
276 // Get theX, theY data set
277 if (mom.Z() < 0)
278 mom *= -1; // particle going backward
279 Double_t theX = TMath::ATan2(mom.X(), mom.Z()) * TMath::RadToDeg();
280 Double_t theY = TMath::ATan2(mom.Y(), mom.Z()) * TMath::RadToDeg();
281 if (stationNr % 2 == 0)
282 theY = -theY; // 28-09-19 - OK
283 Int_t indXY = FindDataSet(theX, theY);
284 // cout << " ------------------ Data set: " << indXY << " " << posIn.Z() << " " << theX << " " << theY << endl;
285 TH2I* hxy = fHistMap[indXY];
286
287 // TTree *clus = (TTree*) gROOT->FindObjectAny("clus");
288 //((TH2I*)gROOT->FindObjectAny("hXY"))->GetRandom2(xrand,yrand);
289 hxy->GetRandom2(xrand, yrand);
290 Float_t xy[2];
291 xy[0] = xrand;
292 xy[1] = yrand;
293
294 Int_t istrip0[2] = {0};
295 multimap<Float_t, Int_t>::iterator mit, mit1;
296
297 for (Int_t ixy = 0; ixy < 2; ++ixy) {
298 map<Int_t, Float_t> charges[2]; // 2 clusters
299 // for (Int_t iclus = 0; iclus < 2; ++iclus) charges[iclus].clear();
300 mit = fXYmap[ixy][indXY].lower_bound(xy[ixy]);
301 if (mit == fXYmap[ixy][indXY].end())
302 --mit;
303 // First cluster
304 Int_t entry = mit->second;
305 fTreeMap[indXY]->GetEntry(entry);
306 Double_t delt = TMath::Abs(fCog[ixy] - xy[ixy]);
307 Int_t ndig = fvDigis[ixy].size();
308 Double_t qTot0 = fQtot[ixy];
309 Double_t dpitch0 = fDpitin[ixy];
310 for (Int_t j = 0; j < ndig; ++j) {
311 charges[0][j + fIbeg[ixy]] = fvDigis[ixy][j];
312 }
313 // cout << cog[ixy] << " " << xy[ixy] << endl;
314 // Second cluster
315 mit1 = mit;
316 if (mit != fXYmap[ixy][indXY].begin())
317 --mit1;
318 else
319 ++mit1;
320 entry = mit1->second;
321 fTreeMap[indXY]->GetEntry(entry);
322 // cout << cog[ixy] << " " << xy[ixy] << endl;
323 Double_t delt1 = TMath::Abs(fCog[ixy] - xy[ixy]);
324 Double_t dist = delt + delt1;
325 // cout << dist << " " << delt << " " << delt1 << endl;
326 Double_t w = (dist - delt) / dist;
327 Double_t w1 = (dist - delt1) / dist;
328 // w1 = 0; w = 1; // !!!!! test
329 // AZ if (gRandom->Rndm() >= 0.5) { w1 = 0; w = 1; }
330 // AZ else { w1 = 1; w = 0; }
331 // cout << w << " " << w1 << " " << w+w1 << " " << (mit==xyMap[ixy].begin()) << endl;
332 ndig = fvDigis[ixy].size();
333 Double_t qTot1 = fQtot[ixy];
334 Double_t dpitch1 = fDpitin[ixy];
335 for (Int_t j = 0; j < ndig; ++j) {
336 charges[1][j + fIbeg[ixy]] = fvDigis[ixy][j];
337 }
338 Int_t i0 = TMath::Min(charges[0].begin()->first, charges[1].begin()->first);
339 // i0 = TMath::Max (i0,0);
340 Int_t i1 = TMath::Max(charges[0].rbegin()->first, charges[1].rbegin()->first);
341 // if (ixy == 0) i1 = TMath::Min (i1,sensor->GetNChannelsFront());
342 // else i1 = TMath::Min (i1,sensor->GetNChannelsBack());
343 /*
344 Double_t cognew[3] = {0}, qqq[3] = {0};
345 for (Int_t iclus = 0; iclus < 2; ++iclus) {
346 for (map<Int_t,Float_t>::iterator it = charges[iclus].begin(); it != charges[iclus].end(); ++it) {
347 cognew[iclus] += it->first * it->second;
348 qqq[iclus] += it->second;
349 }
350 cognew[iclus] /= qqq[iclus];
351 h1[ixy][iclus+1]->Fill(cognew[iclus]);
352 }
353 */
354 // Int_t istr;
355 if (ixy == 0) {
356 if (stationNr % 2)
357 istrip0[ixy] = sensor->GetFrontChannel(posIn.X(), posIn.Y(), posIn.Z(), dPitch);
358 else
359 istrip0[ixy] = sensor->GetFrontChannel(posOut.X(), posOut.Y(), posOut.Z(), dPitch);
360 } else {
361 if (stationNr % 2)
362 istrip0[ixy] = sensor->GetBackChannel(posIn.X(), posIn.Y(), posIn.Z(), dPitch);
363 else
364 istrip0[ixy] = sensor->GetBackChannel(posOut.X(), posOut.Y(), posOut.Z(), dPitch);
365 }
366
367 Double_t pitch = (ixy == 0) ? sensor->GetDx() : sensor->GetDy();
368 dPitch /= pitch;
369 if (dPitch < 0) {
370 --istrip0[ixy];
371 dPitch += 1.0;
372 } // new
373 // if (stationNr % 2 == 0) { istrip0[ixy]-=1; dPitch = 1 - dPitch; }
374 if (stationNr % 2 == 0) {
375 dPitch = 1 - dPitch;
376 }
377 // AZ if (stationNr % 2 == 0) { dPitch = 1 - TMath::Abs(dPitch); }
378 // AZ if (stationNr % 2 == 0) { dPitch = - dPitch; }
379 if (stationNr % 2 == 0)
380 ++istrip0[ixy]; // AZ - how to explain this???
381
382 // Weighted average of 2 clusters
383 Double_t qtot12 = qTot0 * w + qTot1 * w1;
384 Double_t dpitch12 = dpitch0 * w + dpitch1 * w1;
385 // if (stationNr % 2 == 0) { dpitch12 = 1 - dpitch12; } // AZ
386 Int_t ntot = i1 - i0 + 1;
387 Double_t* q12array = new Double_t[ntot];
388 Double_t q12sum0 = 0, q12sum1 = 0;
389
390 for (Int_t j = i0; j <= i1; ++j) {
391 if (charges[0].find(j) == charges[0].end())
392 charges[0][j] = 0;
393 if (charges[1].find(j) == charges[1].end())
394 charges[1][j] = 0;
395 // Double_t q12 = (charges[0][j] * w + charges[1][j] * w1) / (w + w1) * fQtot[ixy];
396 q12array[j - i0] = (charges[0][j] * w + charges[1][j] * w1) * qtot12; // w + w1 = 1
397 q12sum0 += q12array[j - i0];
398 }
399 Double_t shift = dPitch - dpitch12, q12 = 0;
400 // if (stationNr % 2 == 0) shift = -shift; //AZ
401 map<Int_t, Double_t> buffer;
402
403 for (Int_t j = i0; j <= i1; ++j) {
404 if (stationNr % 2)
405 q12 = Interp(q12array, j - i0, ntot, -shift);
406 // else q12 = Interp (q12array, j-i0, ntot, -shift);
407 else
408 q12 = Interp(q12array, j - i0, ntot, -shift);
409
410 // Int_t dj = (ixy == 0) ? j - 241 : j - 256;
411 // Int_t dj = (ixy == 0) ? j - 240 : j - 254;
412 Int_t dj = j - fStrRef[ixy][indXY];
413 buffer[dj] = q12;
414 q12sum1 += q12;
415 // cognew[2] += q12 * j;
416 // qqq[2] += q12;
417 }
418 Double_t scale = q12sum0 / q12sum1;
419 scale *= coeff;
420
421 for (map<Int_t, Double_t>::iterator it = buffer.begin(); it != buffer.end(); ++it) {
422 Int_t jj = it->first;
423 if (stationNr % 2 > 0)
424 jj += istrip0[ixy];
425 else
426 jj = istrip0[ixy] - it->first + 0; // even stations (rotated)
427
428 if (ixy == 0) {
429 if (jj < 0 || jj >= sensor->GetNChannelsFront())
430 continue;
431 stripSignalF[jj] += it->second * scale;
432 chanPointMapF[jj].insert(iPoint);
433 } else {
434 if (jj < 0 || jj >= sensor->GetNChannelsBack())
435 continue;
436 stripSignalB[jj] += it->second * scale;
437 chanPointMapB[jj].insert(iPoint);
438 }
439 }
440 delete[] q12array;
441 // cognew[2] /= qqq[2];
442 // h1[ixy][3]->Fill(cognew[2]);
443 // cogGen[ixy] = cognew[2];
444 } // for (Int_t ixy = 0;
445 } // for (it1 = pSet.begin();
446}
447
448// -------------------------------------------------------------------------
449
450Double_t BmnGemFastDigitize::Interp(Double_t* yy, Int_t indx, Int_t ntot, Double_t dx)
451{
452 static Double_t y0 = 0, y1 = 0;
453 static Int_t indx0 = 0;
454
455 if (indx == 0) {
456 if (dx >= 0) {
457 indx0 = 1;
458 y0 = yy[0];
459 } else {
460 indx0 = 0;
461 y0 = 0.0;
462 }
463 y1 = yy[indx0];
464 } else {
465 y0 = y1;
466 if (indx0 >= ntot)
467 y1 = 0.0;
468 else
469 y1 = yy[indx0];
470 }
471
472 ++indx0;
473 if (dx < 0)
474 dx += 1;
475 // cout << indx << " aaa " << y0 + (y1 - y0) / (x1 - x0) * (x - x0) << " " << yy[indx] << " " << indx0 << endl;
476 // cout << " Interp: " << indx0 << " " << dx << " " << y0 << " " << y0 + (y1 - y0) * dx << endl;
477 // AZ return y0 + (y1 - y0) * dx;
478 return TMath::Max(y0 + (y1 - y0) * dx, 0.0);
479}
480
481// -------------------------------------------------------------------------
482
483Int_t BmnGemFastDigitize::FindDataSet(Double_t theX, Double_t theY)
484{
485 // Find appropriate data set according to theX, theY
486 // (contribution from R.Zinchenko)
487
488 static Int_t first = 1;
489
490 if (first) {
491 first = 0;
492 // Find nearest non-zero bin
493 Int_t nx = fhInd->GetNbinsX(), ny = fhInd->GetNbinsY();
494 vector<Int_t> ia(ny);
495 vector<vector<Int_t>> matri(nx, ia);
496 vector<Float_t> fa(ny);
497 vector<vector<Float_t>> matrf(nx, fa);
498 Int_t ix, iy, /*iz,*/ ix0, iy0, iz0, ixOK, iyOK;
499 std::queue<Int_t> que;
500
501 for (Int_t i = 1; i <= nx; ++i) {
502 for (Int_t j = 1; j <= ny; ++j) {
503 if (fhInd->GetBinContent(i, j)) {
504 matri[i - 1][j - 1] = fhInd->GetBin(i, j);
505 que.push(fhInd->GetBin(i, j));
506 }
507 }
508 }
509
510 while (que.size()) {
511 Int_t bin0 = que.front(); // current cell
512 fhInd->GetBinXYZ(bin0, ix0, iy0, iz0);
513 --ix0;
514 --iy0;
515 Int_t binOK = matri[ix0][iy0]; // active cell
516 fhInd->GetBinXYZ(binOK, ixOK, iyOK, iz0);
517 --ixOK;
518 --iyOK;
519 que.pop();
520
521 // Loop over cell neighbours
522 for (Int_t i = -1; i < 2; ++i) {
523 ix = ix0 + i;
524 if (ix < 0 || ix > nx - 1)
525 continue;
526
527 for (Int_t j = -1; j < 2; ++j) {
528 iy = iy0 + j;
529 if (iy < 0 || iy > ny - 1)
530 continue;
531 if (i == 0 && j == 0)
532 continue;
533 Float_t r20 = matrf[ix][iy];
534 if (matri[ix][iy] && r20 < 0.1)
535 continue; // active cell
536 Float_t r2 = (ix - ixOK) * (ix - ixOK) + (iy - iyOK) * (iy - iyOK);
537 if (r20 < 0.1 || r2 < r20) {
538 matrf[ix][iy] = r2;
539 matri[ix][iy] = matri[ix0][iy0];
540 if (r20 < 0.1)
541 que.push(fhInd->GetBin(ix + 1, iy + 1));
542 }
543 }
544 }
545 }
546
547 for (Int_t i = 1; i <= nx; ++i) {
548 for (Int_t j = 1; j <= ny; ++j) {
549 fhInd->SetBinContent(i, j, matri[i - 1][j - 1]);
550 }
551 }
552 } // if (first)
553
554 Int_t ix = fhInd->GetXaxis()->FindBin(theX);
555 if (ix < 1)
556 ix = 1;
557 else if (ix > fhInd->GetNbinsX())
558 ix = fhInd->GetNbinsX();
559 Int_t iy = fhInd->GetYaxis()->FindBin(theY);
560 if (iy < 1)
561 iy = 1;
562 else if (iy > fhInd->GetNbinsY())
563 iy = fhInd->GetNbinsY();
564
565 return fhInd->GetBinContent(ix, iy);
566}
567
568// -------------------------------------------------------------------------
569
570Double_t BmnGemFastDigitize::GetNumberOfClusters(Double_t beta,
571 Double_t gamma,
572 Double_t charge,
573 Double_t p0,
574 Double_t p1)
575{
576 Double_t beta2 = beta * beta;
577 Double_t gamma2 = gamma * gamma;
578 Double_t val = p0 / beta2 * (p1 + TMath::Log(beta2 * gamma2) - beta2);
579 return TMath::Min(val * charge * charge, 20000.);
580}
581
582// -------------------------------------------------------------------------
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
i
int i
Definition P4_F32vec4.h:22
Int_t
int Int_t
Float_t
float Float_t
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
gDirectory
#define gDirectory
gSystem
R__EXTERN TSystem * gSystem
BmnGemFastDigitize::~BmnGemFastDigitize
~BmnGemFastDigitize()
Destructor.
Definition BmnGemFastDigitize.cxx:50
BmnGemFastDigitize::BmnGemFastDigitize
BmnGemFastDigitize()
Default constructor.
Definition BmnGemFastDigitize.cxx:43
BmnGemFastDigitize::Init
virtual InitStatus Init()
Virtual method Init.
Definition BmnGemFastDigitize.cxx:54
BmnGemFastDigitize::FindDataSet
Int_t FindDataSet(Double_t theX, Double_t theY)
Definition BmnGemFastDigitize.cxx:483
BmnGemFastDigitize::ProduceHitResponseFast
void ProduceHitResponseFast(CbmStsSensor *sensor, std::set< Int_t > &pSet, Double_t *stripSignalF, Double_t *stripSignalB, std::map< Int_t, std::set< Int_t > > &chanPointMapF, std::map< Int_t, std::set< Int_t > > &chanPointMapB)
BmnGemFastDigitize::fStrRef
std::map< Int_t, Int_t > fStrRef[2]
Definition BmnGemFastDigitize.h:64
BmnGemFastDigitize::Interp
Double_t Interp(Double_t *yy, Int_t indx, Int_t ntot, Double_t dx)
Definition BmnGemFastDigitize.cxx:450
BmnGemFastDigitize::Exec
virtual void Exec(Option_t *opt)
Virtual method Exec.
Definition BmnGemFastDigitize.cxx:209
BmnGemFastDigitize::fvDig1
std::vector< Float_t > * fvDig1
Definition BmnGemFastDigitize.h:70
BmnGemFastDigitize::fvDig0
std::vector< Float_t > * fvDig0
Definition BmnGemFastDigitize.h:70
BmnGemFastDigitize::fPoints
TClonesArray * fPoints
Input arrays.
Definition BmnGemFastDigitize.h:57
BmnGemFastDigitize::fvDigis
std::vector< Float_t > fvDigis[2]
Definition BmnGemFastDigitize.h:69
BmnGemFastDigitize::fMCTracks
TClonesArray * fMCTracks
Definition BmnGemFastDigitize.h:58
BmnGemFastDigitize::GetNumberOfClusters
Double_t GetNumberOfClusters(Double_t beta, Double_t gamma, Double_t charge, Double_t p0, Double_t p1)
Definition BmnGemFastDigitize.cxx:570
BmnGemFastDigitize::fTreeMap
std::map< Int_t, TTree * > fTreeMap
Definition BmnGemFastDigitize.h:66
BmnGemFastDigitize::fXYmap
std::map< Int_t, std::multimap< Float_t, Int_t > > fXYmap[2]
Output array of CbmStsHits.
Definition BmnGemFastDigitize.h:63
BmnGemFastDigitize::Finish
virtual void Finish()
Virtual method Finish.
Definition BmnGemFastDigitize.cxx:213
BmnGemFastDigitize::fHistMap
std::map< Int_t, TH2I * > fHistMap
Definition BmnGemFastDigitize.h:65
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Accessors.
Definition CbmMCTrack.h:56
CbmStsPoint::PositionOut
void PositionOut(TVector3 &pos)
Definition CbmStsPoint.h:74
CbmStsPoint::PositionIn
void PositionIn(TVector3 &pos)
Definition CbmStsPoint.h:73
CbmStsSensor::GetNChannelsFront
Int_t GetNChannelsFront() const
Definition CbmStsSensor.h:123
CbmStsSensor::GetDy
Double_t GetDy() const
Definition CbmStsSensor.h:116
CbmStsSensor::GetNChannelsBack
Int_t GetNChannelsBack() const
Definition CbmStsSensor.h:124
CbmStsSensor::GetFrontChannel
Int_t GetFrontChannel(Double_t x, Double_t y, Double_t z, Double_t &dPitch)
Definition CbmStsSensor.cxx:570
CbmStsSensor::GetDx
Double_t GetDx() const
Definition CbmStsSensor.h:115
CbmStsSensor::GetBackChannel
Int_t GetBackChannel(Double_t x, Double_t y, Double_t z, Double_t &dPitch)
Definition CbmStsSensor.cxx:611
CbmStsSensor::GetStationNr
Int_t GetStationNr() const
Definition CbmStsSensor.h:96
TDatabasePDG::Instance
static TDatabasePDG * Instance()
TDatabasePDG::GetParticle
TParticlePDG * GetParticle(const char *name) const
TDirectoryFile::Get
T * Get(const char *namecycle)
TEntryList::GetEntry
virtual Long64_t GetEntry(Long64_t index)
TEntryList::GetN
virtual Long64_t GetN() const
TFile::Open
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
TH1::Clone
TObject * Clone(const char *newname="") const override
TH2::GetRandom2
virtual void GetRandom2(Double_t &x, Double_t &y, TRandom *rng=nullptr)
TLorentzVector::Beta
Double_t Beta() const
TLorentzVector::Gamma
Double_t Gamma() const
TLorentzVector::SetVectM
void SetVectM(const TVector3 &spatial, Double_t mass)
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
TParticlePDG::Charge
Double_t Charge() const
TParticlePDG::Mass
Double_t Mass() const
TString::Atoi
Int_t Atoi() const
TTree::SetBranchStatus
virtual void SetBranchStatus(const char *bname, bool status=true, UInt_t *found=nullptr)
TTree::GetEntry
virtual Int_t GetEntry(Long64_t entry, Int_t getall=0)
TTree::Draw
virtual Long64_t Draw(const char *varexp, const char *selection, Option_t *option="", Long64_t nentries=kMaxEntries, Long64_t firstentry=0)
TTree::SetEntryList
virtual void SetEntryList(TEntryList *list, Option_t *opt="")
TTree::SetBranchAddress
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
TVector3::Z
Double_t Z() const
TVector3::Y
Double_t Y() const
TVector3::X
Double_t X() const
beta
double beta(double x, double y)
TMath::Min
Double_t Min(Double_t a, Double_t b)
TMath::ATan2
Double_t ATan2(Double_t y, Double_t x)
TMath::Log
Double_t Log(Double_t x)
TMath::Abs
Double_t Abs(Double_t d)
TMath::Max
Double_t Max(Double_t a, Double_t b)
TMath::RadToDeg
constexpr Double_t RadToDeg()
nlohmann::find
iterator find(const Key &key)
Definition json.hpp:16128
nlohmann::end
*brief create JSON pointer *sa reference_tokens end()
nlohmann::begin
*brief append another JSON pointer at the end of this JSON pointer *sa ptr reference_tokens begin()