BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnRescale.h
Go to the documentation of this file.
1
2#ifndef BMNRESCALE_H
3#define BMNRESCALE_H
4
5#include <DigiRunHeader.h>
6#include <FairMCEventHeader.h>
7#include <TBranch.h>
8#include <TCanvas.h>
9#include <TChain.h>
10#include <TClass.h>
11#include <TClonesArray.h>
12#include <TFile.h>
13#include <TGeoTrack.h>
14#include <TH1D.h>
15#include <TLeaf.h>
16#include <TStopwatch.h>
17#include <TStyle.h>
18#include <TTree.h>
19#include <math.h>
20#include <regex>
21#include <vector>
22
23// #include <CbmMCTrack.h>
24// #include <CbmStsPoint.h>
25// #include <BmnSiliconPoint.h>
26// #include <BmnCSCPoint.h>
27#include "BmnCSCStationSet.h"
28#include "BmnGemStripConfiguration.h"
29#include "BmnGemStripStationSet.h"
30#include "BmnGemStripStationSet_RunSpring2017.h"
31#include "BmnGemStripStationSet_RunSummer2016.h"
32#include "BmnGemStripStationSet_RunWinter2016.h"
33#include "BmnSiliconStationSet.h"
34
35#include <BmnCSCDigit.h>
36#include <BmnCSCHit.h>
37#include <BmnEnums.h>
38#include <BmnEventHeader.h>
39#include <BmnGemStripDigit.h>
40#include <BmnGemStripHit.h>
41#include <BmnGlobalTrack.h>
42#include <BmnMath.h>
43#include <BmnSiliconDigit.h>
44#include <BmnSiliconHit.h>
45#include <DstEventHeader.h>
46
47#define TRACK_HITS 1
48
49using namespace std;
50using namespace TMath;
51
52class BmnSigInfo
53{
54 public:
55 BmnSigInfo()
56 {
57 minVal = DBL_MAX;
58 maxVal = DBL_MIN;
59 hSig = NULL;
60 hIntSig = NULL;
61 func = NULL;
62 }
63
64 ~BmnSigInfo()
65 {
66 if (hSig)
67 delete hSig;
68 if (hIntSig)
69 delete hIntSig;
70 if (func)
71 delete func;
72 hSig = NULL;
73 hIntSig = NULL;
74 func = NULL;
75 }
76 Double_t minVal;
77 Double_t maxVal;
78 TH1D* hSig;
79 TH1D* hIntSig;
80 TF1* func;
81
82 private:
83};
84
85/**
86 * BmnRescale - class implementing ADC signal rescaling
87 * specific for each strip plane
88 *
89 * \author Ilnur Gabdrakhmanov <ilnur@jinr.ru>
90 */
91class BmnRescale
92{
93 public:
94 /**
95 *
96 * @param periodId
97 * @param runId
98 * @param lowThr lower cluster signal threshold
99 * @param ClusterSizeThr lower cluster size threshold (in strip count)
100 * @param nBins Histograms' binning
101 */
102 BmnRescale(UInt_t periodId, UInt_t runId, Double_t lowThr = 0, Double_t ClusterSizeThr = 0, Int_t nBins = 1e6);
103 virtual ~BmnRescale();
104
105 BmnStatus CreateRescales(TString fileNameMC, TString fileNameEx);
106 TF1* GetRescaleFunc(TString name, TF1* mc, TF1* ex);
107
108 vector<vector<vector<TF1*>>> GetGemRescaleVector() { return fStripGemRescale; }
109
110 vector<vector<vector<TF1*>>> GetSilRescaleVector() { return fStripSilRescale; }
111
112 vector<vector<vector<TF1*>>> GetCSCRescaleVector() { return fStripCSCRescale; }
113
114 vector<vector<vector<BmnSigInfo*>>> GetGemInfoVectorExp() { return fStripGemInfoEx; }
115
116 vector<vector<vector<BmnSigInfo*>>> GetSilInfoVectorExp() { return fStripSilInfoEx; }
117
118 vector<vector<vector<BmnSigInfo*>>> GetCSCInfoVectorExp() { return fStripCSCInfoEx; }
119
120 TCanvas* GetCanvasGem();
121 TCanvas* GetCanvasSil();
122 TCanvas* GetCanvasCSC();
123
124 private:
125 static const UInt_t Pad_Width = 640;
126 static const UInt_t Pad_Height = 360;
127 BmnSetup fSetup;
128 UInt_t fPeriodId;
129 UInt_t fRunId;
130 Double_t fLowThr;
131 Double_t fClusterSizeThr;
132 Int_t fNBins;
133
134 BmnGemStripStationSet* fGemStationSet;
135 BmnSiliconStationSet* fSilStationSet;
136 BmnCSCStationSet* fCscStationSet;
137
138 vector<vector<vector<BmnSigInfo*>>> fStripGemInfoMC;
139 vector<vector<vector<BmnSigInfo*>>> fStripSilInfoMC;
140 vector<vector<vector<BmnSigInfo*>>> fStripCSCInfoMC;
141 vector<vector<vector<BmnSigInfo*>>> fStripGemInfoEx;
142 vector<vector<vector<BmnSigInfo*>>> fStripSilInfoEx;
143 vector<vector<vector<BmnSigInfo*>>> fStripCSCInfoEx;
144 vector<vector<vector<TF1*>>> fStripGemRescale;
145 vector<vector<vector<TF1*>>> fStripSilRescale;
146 vector<vector<vector<TF1*>>> fStripCSCRescale;
147
148 template<class ClStationSet, class ClSt, class ClMod>
149 TCanvas* GetCanvas(ClStationSet* ss)
150 {
151 Int_t sumMods = 0;
152 Int_t maxLayers = 0;
153 for (Int_t iStation = 0; iStation < ss->GetNStations(); iStation++) {
154 ClSt* st = ss->GetStation(iStation);
155 sumMods += st->GetNModules();
156 for (Int_t iModule = 0; iModule < st->GetNModules(); iModule++) {
157 ClMod* mod = st->GetModule(iModule);
158 if (maxLayers < mod->GetNStripLayers())
159 maxLayers = mod->GetNStripLayers();
160 }
161 }
162 TString name = Form("Strip_Digit_%d", rand());
163 TString title = "Canvas";
164 TCanvas* can = new TCanvas(name, title, Pad_Width * maxLayers, Pad_Height * sumMods);
165 can->Divide(maxLayers, sumMods);
166 return can;
167 }
168
169 BmnStatus LoadDigitDistributionsFromFile(TString name,
170 vector<TString> digiNames,
171 vector<vector<vector<vector<BmnSigInfo*>>>> digiVecs);
172 BmnStatus LoadClusterDistributionsFromFile(TString name,
173 vector<TString> clusterNames,
174 vector<vector<vector<vector<BmnSigInfo*>>>> clusterVecs);
175 BmnStatus LoadClustersOfTrackDistributionsFromFile(TString fileName,
176 TString globalTrackNames,
177 vector<TString> subTrackNames,
178 vector<TString> hitNames,
179 vector<vector<vector<vector<BmnSigInfo*>>>> clusterVecs);
180 void CreateStripVectors();
181
182 template<class cl>
183 void FreeVector3D(vector<vector<vector<cl*>>> vvv)
184 {
185 for (auto& row : vvv)
186 for (auto& col : row)
187 for (auto& el : col)
188 delete el;
189 }
190
191 /**
192 * Template function for generating 3D vector of BmnSigInfo objects
193 * @param ss station set for the specific detector subsystem
194 * @param vvv empty 3D vector to fill
195 */
196 template<class ClStationSet, class ClSt, class ClMod>
197 void FillInfoVector(ClStationSet* ss, vector<vector<vector<BmnSigInfo*>>>& vvv)
198 {
199 for (Int_t iStation = 0; iStation < ss->GetNStations(); iStation++) {
200 vector<vector<BmnSigInfo*>> rowGEM;
201 ClSt* st = ss->GetStation(iStation);
202 for (Int_t iModule = 0; iModule < st->GetNModules(); iModule++) {
203 vector<BmnSigInfo*> colGEM;
204 ClMod* mod = st->GetModule(iModule);
205 for (Int_t iLayer = 0; iLayer < mod->GetNStripLayers(); iLayer++) {
206 BmnSigInfo* r = new BmnSigInfo();
207 colGEM.push_back(r);
208 }
209 rowGEM.push_back(colGEM);
210 }
211 vvv.push_back(rowGEM);
212 }
213 }
214
215 template<class ClStationSet, class ClSt, class ClMod>
216 void FillRescaleVector(ClStationSet* ss,
217 vector<vector<vector<BmnSigInfo*>>>& vMC,
218 vector<vector<vector<BmnSigInfo*>>>& vEx,
219 vector<vector<vector<TF1*>>>& vRescale)
220 {
221 for (Int_t iStation = 0; iStation < ss->GetNStations(); iStation++) {
222 vector<vector<TF1*>> rowGEM;
223 ClSt* st = ss->GetStation(iStation);
224 for (Int_t iModule = 0; iModule < st->GetNModules(); iModule++) {
225 vector<TF1*> colGEM;
226 ClMod* mod = st->GetModule(iModule);
227 for (Int_t iLayer = 0; iLayer < mod->GetNStripLayers(); iLayer++) {
228 TF1* r = nullptr;
229 if ((vMC[iStation][iModule][iLayer]->func) && (vEx[iStation][iModule][iLayer]->func))
230 r = GetRescaleFunc(Form("rescale_%s", vEx[iStation][iModule][iLayer]->func->GetName()),
231 vMC[iStation][iModule][iLayer]->func, vEx[iStation][iModule][iLayer]->func);
232 colGEM.push_back(r);
233 printf("st %d mod %d l %d func %08X of MC %08X Ex %08X \n", iStation, iModule, iLayer,
234 (UInt_t)(size_t)r, (UInt_t)(size_t)vMC[iStation][iModule][iLayer]->func,
235 (UInt_t)(size_t)vEx[iStation][iModule][iLayer]->func);
236 }
237 rowGEM.push_back(colGEM);
238 }
239 vRescale.push_back(rowGEM);
240 }
241 }
242 void CreateDigitDistributions(TTree* treeDig, TClonesArray* digits, vector<vector<vector<BmnSigInfo*>>>& infoVec);
243
244 /**
245 * Template function aimed to create exp cluster signal distribution function (Hits from tracks variant)
246 * @param treeDST an experimental DST tree
247 * @param globTracks
248 * @param subTracks
249 * @param gemHits inner tracker subdetector hit branch
250 * @param ss subdetector station set
251 * @param infoVec 3-d vector to store function info
252 */
253 template<class ClHit, class ClStationSet, class ClSt, class ClMod, class ClLay>
254 void CreateClusterDistributions(TTree* treeDST,
255 TClonesArray* globTracks,
256 TClonesArray* subTracks,
257 TClonesArray* gemHits,
258 ClStationSet* ss,
259 vector<vector<vector<BmnSigInfo*>>>& infoVec)
260 {
261 printf("CreateClusterDistributions started for %s\n", gemHits->GetClass()->GetName());
262 TStopwatch timer;
263 timer.Start();
264
265 Bool_t isGem = (typeid(ss) == typeid(BmnGemStripStationSet*));
266 Bool_t isSil = (typeid(ss) == typeid(BmnSiliconStationSet*));
267 Bool_t isCSC = (typeid(ss) == typeid(BmnCSCStationSet*));
268
269 // Find bounds
270 for (Long64_t iEv = 0; iEv < treeDST->GetEntries(); iEv++) {
271 treeDST->GetEntry(iEv);
272
273 for (Int_t iTrack = 0; iTrack < globTracks->GetEntriesFast(); iTrack++) {
274 BmnGlobalTrack* track = (BmnGlobalTrack*)globTracks->UncheckedAt(iTrack);
275 BmnTrack* subTrack = nullptr;
276 if (isGem && (track->GetGemTrackIndex() != -1))
277 subTrack = static_cast<BmnTrack*>(subTracks->UncheckedAt(track->GetGemTrackIndex()));
278 // if (isSil && (track->GetSilTrackIndex() != -1))
279 // subTrack = static_cast<BmnTrack*> (subTracks->UncheckedAt(track->GetSilTrackIndex()));
280 if ((!subTrack) && (!isCSC))
281 continue;
282 vector<Int_t> vrange; // range of hit indices
283 // vrange.resize(
284 // isCSC ? 1 : subTrack->GetNHits(),
285 // isCSC ? track->GetCscHitIndex(0) : 0);
286 if (!isCSC)
287 for (size_t i = 0; i < vrange.size(); i++)
288 vrange[i] = subTrack->GetHitIndex(i);
289 else
290 // if (track->GetCscHitIndex(0) == -1)
291 continue;
292 for (Int_t iHit : vrange) {
293 ClHit* hit = static_cast<ClHit*>(gemHits->UncheckedAt(iHit));
294 Int_t iSt = hit->GetStation();
295 Int_t iMod = hit->GetModule();
296 // if (iMod % 2 == 0)
297 // printf("ist %d imod %d ( %f : %f )\n", iSt, iMod, -hit->GetX(), hit->GetY());
298 ClMod* mod = ss->GetStation(iSt)->GetModule(iMod);
299 // printf("iev %d\n\t iHit %d iSt %d iMod %d\n",
300 // i, iHit, iSt, iMod);
301 for (size_t iLay = 0; iLay < mod->GetStripLayers().size(); iLay++) {
302 ClLay lay = mod->GetStripLayer(iLay);
303 Double_t x = hit->GetX();
304 Double_t y = hit->GetY();
305 // if (iMod % 2 == 0)
306 // printf("iLay %d xmin %f xmax %f ymin %f ymax %f\n", iLay,
307 // lay.GetXMinLayer(), lay.GetXMaxLayer(), lay.GetYMinLayer(),
308 // lay.GetYMaxLayer());
309 if ((lay.IsPointInsideStripLayer(x * (-1), y))) {
310 // printf("\t\tinside %d layer\n", iLay);
311 Int_t iLayCorr = (isSil && (iSt > 0)) ? (iLay / 2) : iLay; // For divided silicon layers
312 BmnSigInfo* info = infoVec[iSt][iMod][iLayCorr];
313 Double_t val = (lay.GetType() == UpperStripLayer) ? hit->GetStripTotalSignalInUpperLayer()
314 : hit->GetStripTotalSignalInLowerLayer();
315 Double_t cs = (lay.GetType() == UpperStripLayer) ? hit->GetClusterSizeInUpperLayer()
316 : hit->GetClusterSizeInLowerLayer();
317 // printf("\tiLay %d \n\t\tval %f fLowThr %f cs %f fClusterSizeThr
318 // %f\n",
319 // iLay, val, fLowThr, cs, fClusterSizeThr);
320 if ((val > fLowThr) && (cs > fClusterSizeThr)) {
321 if (info->maxVal < val)
322 info->maxVal = val;
323 if (info->minVal > val)
324 info->minVal = val;
325 }
326 }
327 }
328 }
329 }
330 }
331 // Create histograms
332 for (size_t iStation = 0; iStation < infoVec.size(); iStation++) {
333 printf("iStation %zu\n", iStation);
334 for (size_t iModule = 0; iModule < infoVec[iStation].size(); iModule++) {
335 printf("\tiModule %zu\n", iModule);
336 for (size_t iLayer = 0; iLayer < infoVec[iStation][iModule].size(); iLayer++) {
337 printf("\t\tiLayer %zu\n", iLayer);
338 BmnSigInfo* info = infoVec[iStation][iModule][iLayer];
339 Double_t minVal = info->minVal;
340 Double_t maxVal = info->maxVal;
341 printf("\t\t\tMin = %f max = %f\n", minVal, maxVal);
342 if (minVal > maxVal)
343 continue;
344 TString name = Form("Sig_exp_%s_%zu_%zu_%zu", gemHits->GetName(), iStation, iModule, iLayer);
345 info->hSig = new TH1D(name, name, fNBins, minVal - 0.5 * (maxVal - minVal) / (Double_t)fNBins,
346 maxVal + 0.5 * (maxVal - minVal) / (Double_t)fNBins);
347 name = Form("IntSignal_%s_%zu_%zu_%zu", gemHits->GetName(), iStation, iModule, iLayer);
348 info->hIntSig = new TH1D(name, name, fNBins, minVal - 0.5 * (maxVal - minVal) / (Double_t)fNBins,
349 maxVal + 0.5 * (maxVal - minVal) / (Double_t)fNBins);
350 }
351 }
352 }
353 // Fill signal hists
354 for (Long64_t iEv = 0; iEv < treeDST->GetEntries(); iEv++) {
355 treeDST->GetEntry(iEv);
356 for (Int_t iTrack = 0; iTrack < globTracks->GetEntriesFast(); iTrack++) {
357 BmnGlobalTrack* track = (BmnGlobalTrack*)globTracks->UncheckedAt(iTrack);
358 BmnTrack* subTrack = nullptr;
359 if (isGem && (track->GetGemTrackIndex() != -1))
360 subTrack = static_cast<BmnTrack*>(subTracks->UncheckedAt(track->GetGemTrackIndex()));
361 // if (isSil && (track->GetSilTrackIndex() != -1))
362 // subTrack = static_cast<BmnTrack*> (subTracks->UncheckedAt(track->GetSilTrackIndex()));
363 if ((!subTrack) && (!isCSC))
364 continue;
365 vector<Int_t> vrange; // range of hit indices
366 // vrange.resize(
367 // isCSC ? 1 : subTrack->GetNHits(),
368 // isCSC ? track->GetCscHitIndex(0) : 0);
369 if (!isCSC)
370 for (size_t i = 0; i < vrange.size(); i++)
371 vrange[i] = subTrack->GetHitIndex(i);
372 // else
373 // if (track->GetCscHitIndex(0) == -1)
374 // continue;
375 for (Int_t iHit : vrange) {
376 ClHit* hit = static_cast<ClHit*>(gemHits->UncheckedAt(iHit));
377 // printf("iev %lld ihit %d flag %d\n ", i, iHit, hit->GetFlag());
378 if (hit->GetFlag() == kFALSE) ///<- prevent multiple count
379 {
380 // printf("iev %lld ihit %d flag %d\n ", i, iHit, hit->GetFlag());
381 continue;
382 }
383 Int_t iSt = hit->GetStation();
384 Int_t iMod = hit->GetModule();
385 ClMod* mod = ss->GetStation(iSt)->GetModule(iMod);
386 for (size_t iLay = 0; iLay < mod->GetStripLayers().size(); iLay++) {
387 ClLay lay = mod->GetStripLayer(iLay);
388 Double_t x = hit->GetX();
389 Double_t y = hit->GetY();
390 if (lay.IsPointInsideStripLayer(x * (-1), y)) {
391 hit->SetFlag(kFALSE);
392 // ClHit *hit2 = static_cast<ClHit*> (gemHits->UncheckedAt(hitIndex));
393 // printf("iev %lld itrack %d ihit %d hitIndex %d set flag %d\n ",
394 // i, iTrack, iHit, hitIndex, hit2->GetFlag());
395 Int_t iLayCorr = (isSil && (iSt > 0)) ? (iLay / 2) : iLay; // For divided silicon layers
396 BmnSigInfo* info = infoVec[iSt][iMod][iLayCorr];
397 if (!(info->hSig))
398 continue;
399 Double_t val = (lay.GetType() == UpperStripLayer) ? hit->GetStripTotalSignalInUpperLayer()
400 : hit->GetStripTotalSignalInLowerLayer();
401 Int_t cs = (lay.GetType() == UpperStripLayer) ? hit->GetClusterSizeInUpperLayer()
402 : hit->GetClusterSizeInLowerLayer();
403 if ((val > fLowThr) && (cs > fClusterSizeThr)) {
404 info->hSig->Fill(val);
405 }
406 }
407 }
408 }
409 }
410 // clear marks
411 for (Int_t iHit = 0; iHit < gemHits->GetEntriesFast(); iHit++) {
412 BmnHit* hit = static_cast<BmnHit*>(gemHits->UncheckedAt(iHit));
413 hit->SetFlag(kTRUE);
414 }
415 }
416 // Fill integral hists & create funcs
417 for (size_t iStation = 0; iStation < infoVec.size(); iStation++) {
418 for (size_t iModule = 0; iModule < infoVec[iStation].size(); iModule++) {
419 for (size_t iLayer = 0; iLayer < infoVec[iStation][iModule].size(); iLayer++) {
420 BmnSigInfo* info = infoVec[iStation][iModule][iLayer];
421 if (!(info->hSig))
422 continue;
423 Double_t bc = 0;
424 for (Int_t i = 0; i < fNBins; i++) {
425 bc += info->hSig->GetBinContent(i);
426 info->hIntSig->SetBinContent(i, bc);
427 }
428 if (bc == 0)
429 continue;
430 info->hIntSig->Scale(1 / bc);
431 info->hIntSig->SetLineColor(kRed);
432 info->func = new TF1(
433 Form("%s_%zu_%zu_%zu", gemHits->GetName(), iStation, iModule, iLayer),
434 [info, bc](Double_t* x, Double_t* p) {
435 // Double_t val = (bc - x[0])/bc;//log(x[0]);
436 if (x[0] >= info->maxVal)
437 return 1.0;
438 Int_t iBin = info->hIntSig->GetXaxis()->FindBin(x[0]);
439 Int_t jBin = iBin + ((x[0] > info->hIntSig->GetBinCenter(iBin)) ? 1 : -1);
440 Double_t v = info->hIntSig->GetBinContent(iBin)
441 + (info->hIntSig->GetBinContent(jBin) - info->hIntSig->GetBinContent(iBin)) * 2
442 * abs(x[0] - info->hIntSig->GetBinCenter(iBin))
443 / (info->hIntSig->GetBinWidth(iBin) + info->hIntSig->GetBinWidth(jBin));
444 return v;
445 },
446 info->minVal, info->maxVal, 0);
447 info->func->SetNpx(fNBins);
448 info->func->SetLineColor(kBlue);
449 }
450 }
451 }
452 timer.Stop();
453 Double_t rtime = timer.RealTime();
454 Double_t ctime = timer.CpuTime();
455 printf("Real time %f s,\n CPU time %f s\n", rtime, ctime);
456 }
457
458 /**
459 * Template function aimed to create exp cluster signal distribution function (all hits)
460 * @param treeDST an experimental DST tree
461 * @param gemHits inner tracker subdetector hit branch
462 * @param ss subdetector station set
463 * @param infoVec 3-d vector to store function info
464 */
465 template<class ClHit, class ClStationSet, class ClSt, class ClMod, class ClLay>
466 void CreateClusterDistributions(TTree* treeDST,
467 TClonesArray* gemHits,
468 ClStationSet* ss,
469 vector<vector<vector<BmnSigInfo*>>>& infoVec)
470 {
471 printf("CreateClusterDistributions started for %s\n", gemHits->GetClass()->GetName());
472
473 Bool_t isGem = (typeid(ss) == typeid(BmnGemStripStationSet*));
474 printf("isGem = %d \n", isGem);
475 Bool_t isSil = (typeid(ss) == typeid(BmnSiliconStationSet*));
476 printf("isSil = %d \n", isSil);
477
478 // Find bounds
479 for (Long64_t iEv = 0; iEv < treeDST->GetEntries(); iEv++) {
480 treeDST->GetEntry(iEv);
481 for (Int_t iHit = 0; iHit < gemHits->GetEntriesFast(); iHit++) {
482 // for (Int_t iHit = 0; iHit < gemHits->GetEntriesFast(); iHit++) {
483 ClHit* hit = static_cast<ClHit*>(gemHits->UncheckedAt(iHit));
484 Int_t iSt = hit->GetStation();
485 Int_t iMod = hit->GetModule();
486 // if (iMod % 2 == 0)
487 // printf("ist %d imod %d ( %f : %f )\n", iSt, iMod, -hit->GetX(), hit->GetY());
488 ClMod* mod = ss->GetStation(iSt)->GetModule(iMod);
489 // printf("iev %d\n\t iHit %d iSt %d iMod %d\n",
490 // i, iHit, iSt, iMod);
491 for (size_t iLay = 0; iLay < mod->GetStripLayers().size(); iLay++) {
492 ClLay lay = mod->GetStripLayer(iLay);
493 Double_t x = hit->GetX();
494 Double_t y = hit->GetY();
495 // if (iMod % 2 == 0)
496 // printf("iLay %d xmin %f xmax %f ymin %f ymax %f\n", iLay,
497 // lay.GetXMinLayer(), lay.GetXMaxLayer(), lay.GetYMinLayer(),
498 // lay.GetYMaxLayer());
499 if ((lay.IsPointInsideStripLayer(x * (-1), y))) {
500 // printf("\t\tinside %d layer\n", iLay);
501 Int_t iLayCorr = (isSil && (iSt > 0)) ? (iLay / 2) : iLay; // For divided silicon layers
502 BmnSigInfo* info = infoVec[iSt][iMod][iLayCorr];
503 Double_t val = (lay.GetType() == UpperStripLayer) ? hit->GetStripTotalSignalInUpperLayer()
504 : hit->GetStripTotalSignalInLowerLayer();
505 Double_t cs = (lay.GetType() == UpperStripLayer) ? hit->GetClusterSizeInUpperLayer()
506 : hit->GetClusterSizeInLowerLayer();
507 // printf("\tiLay %d \n\t\tval %f fLowThr %f cs %f fClusterSizeThr %f\n",
508 // iLay, val, fLowThr, cs, fClusterSizeThr);
509 if ((val > fLowThr) && (cs > fClusterSizeThr)) {
510 if (info->maxVal < val)
511 info->maxVal = val;
512 if (info->minVal > val)
513 info->minVal = val;
514 }
515 }
516 }
517 }
518 }
519 // Create histograms
520 for (size_t iStation = 0; iStation < infoVec.size(); iStation++) {
521 printf("iStation %zu\n", iStation);
522 for (size_t iModule = 0; iModule < infoVec[iStation].size(); iModule++) {
523 printf("\tiModule %zu\n", iModule);
524 for (size_t iLayer = 0; iLayer < infoVec[iStation][iModule].size(); iLayer++) {
525 printf("\t\tiLayer %zu\n", iLayer);
526 BmnSigInfo* info = infoVec[iStation][iModule][iLayer];
527 Double_t minVal = info->minVal;
528 Double_t maxVal = info->maxVal;
529 printf("\t\t\tMin = %f max = %f\n", minVal, maxVal);
530 if (minVal > maxVal)
531 continue;
532 TString name = Form("Sig_exp_%s_%zu_%zu_%zu", gemHits->GetName(), iStation, iModule, iLayer);
533 info->hSig = new TH1D(name, name, fNBins, minVal - 0.5 * (maxVal - minVal) / (Double_t)fNBins,
534 maxVal + 0.5 * (maxVal - minVal) / (Double_t)fNBins);
535 name = Form("IntSignal_%s_%zu_%zu_%zu", gemHits->GetName(), iStation, iModule, iLayer);
536 info->hIntSig = new TH1D(name, name, fNBins, minVal - 0.5 * (maxVal - minVal) / (Double_t)fNBins,
537 maxVal + 0.5 * (maxVal - minVal) / (Double_t)fNBins);
538 }
539 }
540 }
541 // Fill signal hists
542 for (Long64_t iEv = 0; iEv < treeDST->GetEntries(); iEv++) {
543 treeDST->GetEntry(iEv);
544 for (Int_t iHit = 0; iHit < gemHits->GetEntriesFast(); iHit++) {
545 ClHit* hit = static_cast<ClHit*>(gemHits->UncheckedAt(iHit));
546 Int_t iSt = hit->GetStation();
547 Int_t iMod = hit->GetModule();
548 ClMod* mod = ss->GetStation(iSt)->GetModule(iMod);
549 for (size_t iLay = 0; iLay < mod->GetStripLayers().size(); iLay++) {
550 ClLay lay = mod->GetStripLayer(iLay);
551 Double_t x = hit->GetX();
552 Double_t y = hit->GetY();
553 if (lay.IsPointInsideStripLayer(x * (-1), y)) {
554 Int_t iLayCorr = (isSil && (iSt > 0)) ? (iLay / 2) : iLay; // For divided silicon layers
555 BmnSigInfo* info = infoVec[iSt][iMod][iLayCorr];
556 if (!(info->hSig))
557 continue;
558 Double_t val = (lay.GetType() == UpperStripLayer) ? hit->GetStripTotalSignalInUpperLayer()
559 : hit->GetStripTotalSignalInLowerLayer();
560 Int_t cs = (lay.GetType() == UpperStripLayer) ? hit->GetClusterSizeInUpperLayer()
561 : hit->GetClusterSizeInLowerLayer();
562 if ((val > fLowThr) && (cs > fClusterSizeThr)) {
563 info->hSig->Fill(val);
564 }
565 }
566 }
567 }
568 }
569 // Fill integral hists & create funcs
570 for (size_t iStation = 0; iStation < infoVec.size(); iStation++) {
571 for (size_t iModule = 0; iModule < infoVec[iStation].size(); iModule++) {
572 for (size_t iLayer = 0; iLayer < infoVec[iStation][iModule].size(); iLayer++) {
573 BmnSigInfo* info = infoVec[iStation][iModule][iLayer];
574 if (!(info->hSig))
575 continue;
576 Double_t bc = 0;
577 for (Int_t i = 0; i < fNBins; i++) {
578 bc += info->hSig->GetBinContent(i);
579 info->hIntSig->SetBinContent(i, bc);
580 }
581 if (bc == 0)
582 continue;
583 info->hIntSig->Scale(1 / bc);
584 info->hIntSig->SetLineColor(kRed);
585 info->func = new TF1(
586 Form("%s_%zu_%zu_%zu", gemHits->GetName(), iStation, iModule, iLayer),
587 [info, bc](Double_t* x, Double_t* p) {
588 // Double_t val = (bc - x[0])/bc;//log(x[0]);
589 if (x[0] >= info->maxVal)
590 return 1.0;
591 Int_t iBin = info->hIntSig->GetXaxis()->FindBin(x[0]);
592 Int_t jBin = iBin + ((x[0] > info->hIntSig->GetBinCenter(iBin)) ? 1 : -1);
593 Double_t v = info->hIntSig->GetBinContent(iBin)
594 + (info->hIntSig->GetBinContent(jBin) - info->hIntSig->GetBinContent(iBin)) * 2
595 * abs(x[0] - info->hIntSig->GetBinCenter(iBin))
596 / (info->hIntSig->GetBinWidth(iBin) + info->hIntSig->GetBinWidth(jBin));
597 return v;
598 },
599 info->minVal, info->maxVal, 0);
600 info->func->SetNpx(fNBins);
601 info->func->SetLineColor(kBlue);
602 }
603 }
604 }
605 }
606 // ClassDef(BmnRescale, 1);
607};
608
609#endif /* BMNRESCALE_H */
can
TCanvas * can(nullptr)
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
Long64_t
long long Long64_t
kTRUE
constexpr Bool_t kTRUE
kRed
kRed
kBlue
kBlue
Form
char * Form(const char *fmt,...)
BmnStatus
BmnStatus
Define status code.
Definition BmnEnums.h:93
UpperStripLayer
@ UpperStripLayer
Definition BmnStripData.h:13
BmnGlobalTrack
BmnGlobalTrack.h.
Definition BmnGlobalTrack.h:14
BmnGlobalTrack::GetGemTrackIndex
Int_t GetGemTrackIndex() const
Accessors.
Definition BmnGlobalTrack.h:34
BmnHit::SetFlag
void SetFlag(Bool_t fl)
Definition BmnHit.h:48
BmnRescale
BmnRescale - class implementing ADC signal rescaling specific for each strip plane.
Definition BmnRescale.h:92
BmnRescale::CreateRescales
BmnStatus CreateRescales(TString fileNameMC, TString fileNameEx)
BmnRescale::~BmnRescale
virtual ~BmnRescale()
BmnRescale::CreateStripVectors
void CreateStripVectors()
BmnRescale::GetSilRescaleVector
vector< vector< vector< TF1 * > > > GetSilRescaleVector()
Definition BmnRescale.h:110
BmnRescale::Pad_Width
static const UInt_t Pad_Width
Definition BmnRescale.h:125
BmnRescale::FreeVector3D
void FreeVector3D(vector< vector< vector< cl * > > > vvv)
Definition BmnRescale.h:183
BmnRescale::fStripGemInfoEx
vector< vector< vector< BmnSigInfo * > > > fStripGemInfoEx
Definition BmnRescale.h:141
BmnRescale::GetRescaleFunc
TF1 * GetRescaleFunc(TString name, TF1 *mc, TF1 *ex)
BmnRescale::CreateClusterDistributions
void CreateClusterDistributions(TTree *treeDST, TClonesArray *globTracks, TClonesArray *subTracks, TClonesArray *gemHits, ClStationSet *ss, vector< vector< vector< BmnSigInfo * > > > &infoVec)
Template function aimed to create exp cluster signal distribution function (Hits from tracks variant)...
Definition BmnRescale.h:254
BmnRescale::fRunId
UInt_t fRunId
Definition BmnRescale.h:129
BmnRescale::GetSilInfoVectorExp
vector< vector< vector< BmnSigInfo * > > > GetSilInfoVectorExp()
Definition BmnRescale.h:116
BmnRescale::fPeriodId
UInt_t fPeriodId
Definition BmnRescale.h:128
BmnRescale::fStripGemRescale
vector< vector< vector< TF1 * > > > fStripGemRescale
Definition BmnRescale.h:144
BmnRescale::LoadDigitDistributionsFromFile
BmnStatus LoadDigitDistributionsFromFile(TString name, vector< TString > digiNames, vector< vector< vector< vector< BmnSigInfo * > > > > digiVecs)
BmnRescale::GetGemInfoVectorExp
vector< vector< vector< BmnSigInfo * > > > GetGemInfoVectorExp()
Definition BmnRescale.h:114
BmnRescale::LoadClusterDistributionsFromFile
BmnStatus LoadClusterDistributionsFromFile(TString name, vector< TString > clusterNames, vector< vector< vector< vector< BmnSigInfo * > > > > clusterVecs)
BmnRescale::fStripCSCRescale
vector< vector< vector< TF1 * > > > fStripCSCRescale
Definition BmnRescale.h:146
BmnRescale::CreateDigitDistributions
void CreateDigitDistributions(TTree *treeDig, TClonesArray *digits, vector< vector< vector< BmnSigInfo * > > > &infoVec)
BmnRescale::fCscStationSet
BmnCSCStationSet * fCscStationSet
Definition BmnRescale.h:136
BmnRescale::fNBins
Int_t fNBins
Definition BmnRescale.h:132
BmnRescale::FillInfoVector
void FillInfoVector(ClStationSet *ss, vector< vector< vector< BmnSigInfo * > > > &vvv)
Template function for generating 3D vector of BmnSigInfo objects.
Definition BmnRescale.h:197
BmnRescale::GetGemRescaleVector
vector< vector< vector< TF1 * > > > GetGemRescaleVector()
Definition BmnRescale.h:108
BmnRescale::GetCanvas
TCanvas * GetCanvas(ClStationSet *ss)
Definition BmnRescale.h:149
BmnRescale::CreateClusterDistributions
void CreateClusterDistributions(TTree *treeDST, TClonesArray *gemHits, ClStationSet *ss, vector< vector< vector< BmnSigInfo * > > > &infoVec)
Template function aimed to create exp cluster signal distribution function (all hits).
Definition BmnRescale.h:466
BmnRescale::GetCanvasGem
TCanvas * GetCanvasGem()
BmnRescale::fStripCSCInfoEx
vector< vector< vector< BmnSigInfo * > > > fStripCSCInfoEx
Definition BmnRescale.h:143
BmnRescale::Pad_Height
static const UInt_t Pad_Height
Definition BmnRescale.h:126
BmnRescale::GetCSCInfoVectorExp
vector< vector< vector< BmnSigInfo * > > > GetCSCInfoVectorExp()
Definition BmnRescale.h:118
BmnRescale::GetCSCRescaleVector
vector< vector< vector< TF1 * > > > GetCSCRescaleVector()
Definition BmnRescale.h:112
BmnRescale::GetCanvasCSC
TCanvas * GetCanvasCSC()
BmnRescale::FillRescaleVector
void FillRescaleVector(ClStationSet *ss, vector< vector< vector< BmnSigInfo * > > > &vMC, vector< vector< vector< BmnSigInfo * > > > &vEx, vector< vector< vector< TF1 * > > > &vRescale)
Definition BmnRescale.h:216
BmnRescale::fLowThr
Double_t fLowThr
Definition BmnRescale.h:130
BmnRescale::fStripSilInfoMC
vector< vector< vector< BmnSigInfo * > > > fStripSilInfoMC
Definition BmnRescale.h:139
BmnRescale::LoadClustersOfTrackDistributionsFromFile
BmnStatus LoadClustersOfTrackDistributionsFromFile(TString fileName, TString globalTrackNames, vector< TString > subTrackNames, vector< TString > hitNames, vector< vector< vector< vector< BmnSigInfo * > > > > clusterVecs)
BmnRescale::fStripGemInfoMC
vector< vector< vector< BmnSigInfo * > > > fStripGemInfoMC
Definition BmnRescale.h:138
BmnRescale::GetCanvasSil
TCanvas * GetCanvasSil()
BmnRescale::fGemStationSet
BmnGemStripStationSet * fGemStationSet
Definition BmnRescale.h:134
BmnRescale::BmnRescale
BmnRescale(UInt_t periodId, UInt_t runId, Double_t lowThr=0, Double_t ClusterSizeThr=0, Int_t nBins=1e6)
BmnRescale::fStripCSCInfoMC
vector< vector< vector< BmnSigInfo * > > > fStripCSCInfoMC
Definition BmnRescale.h:140
BmnRescale::fSilStationSet
BmnSiliconStationSet * fSilStationSet
Definition BmnRescale.h:135
BmnRescale::fStripSilRescale
vector< vector< vector< TF1 * > > > fStripSilRescale
Definition BmnRescale.h:145
BmnRescale::fStripSilInfoEx
vector< vector< vector< BmnSigInfo * > > > fStripSilInfoEx
Definition BmnRescale.h:142
BmnRescale::fSetup
BmnSetup fSetup
Definition BmnRescale.h:127
BmnRescale::fClusterSizeThr
Double_t fClusterSizeThr
Definition BmnRescale.h:131
BmnSigInfo::maxVal
Double_t maxVal
Definition BmnRescale.h:77
BmnSigInfo::func
TF1 * func
Definition BmnRescale.h:80
BmnSigInfo::hIntSig
TH1D * hIntSig
Definition BmnRescale.h:79
BmnSigInfo::minVal
Double_t minVal
Definition BmnRescale.h:76
BmnSigInfo::hSig
TH1D * hSig
Definition BmnRescale.h:78
BmnTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition BmnTrack.h:51
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
TClonesArray::GetClass
TClass * GetClass() const
TCollection::GetName
const char * GetName() const override
TF1::SetNpx
virtual void SetNpx(Int_t npx=100)
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
TH1::Scale
virtual void Scale(Double_t c1=1, Option_t *option="")
TNamed::GetName
const char * GetName() const override
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
TStopwatch::RealTime
Double_t RealTime()
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
TStopwatch::CpuTime
Double_t CpuTime()
TStopwatch::Stop
void Stop()
TTree::GetEntry
virtual Int_t GetEntry(Long64_t entry, Int_t getall=0)
TTree::GetEntries
virtual Long64_t GetEntries() const
std
Definition json.hpp:4446