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