CbmRoot
Loading...
Searching...
No Matches
CbmTrdHitRateFastQa.cxx
Go to the documentation of this file.
1/* Copyright (C) 2013-2021 Institut fuer Kernphysik, Westfaelische Wilhelms-Universitaet Muenster, Muenster
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Cyrano Bergmann [committer], David Emschermann */
4
5#include "CbmTrdHitRateFastQa.h"
6
7#include "CbmTrdAddress.h"
8#include "CbmTrdGeoHandler.h"
9#include "CbmTrdParAsic.h"
10#include "CbmTrdParModAsic.h"
11#include "CbmTrdParModDigi.h"
12#include "CbmTrdParModGeo.h"
13#include "CbmTrdParSetAsic.h"
14#include "CbmTrdParSetDigi.h"
15#include "CbmTrdParSetGeo.h"
16#include "CbmTrdRadiator.h"
17#include "CbmTrdUtils.h"
18
19#include "FairRootManager.h"
20#include "FairRunAna.h"
21#include "FairRuntimeDb.h"
22#include <Logger.h>
23
24#include "TBox.h"
25#include "TCanvas.h"
26#include "TClonesArray.h"
27#include "TColor.h"
28#include "TGeoManager.h"
29#include "TH2.h"
30#include "TImage.h"
31#include "TLine.h"
32#include "TMath.h"
33#include "TPRegexp.h"
34#include "TPolyMarker.h"
35#include "TStyle.h"
36#include <TFile.h>
37
38#include <fstream>
39#include <iostream>
40#include <vector>
41
42#include <cmath>
43
44#include "Riostream.h"
45
46#define winsize 6000 // 1500 // 5500 // 6000
47
48using std::cin;
49using std::cout;
50using std::endl;
51using std::flush;
52using std::ios;
53using std::pair;
54using std::setfill;
55using std::setiosflags;
56using std::setprecision;
57using std::setw;
58using std::vector;
59
60// ---- Default constructor -------------------------------------------
61CbmTrdHitRateFastQa::CbmTrdHitRateFastQa() : CbmTrdHitRateFastQa("TrdHitRateFastQa", "") {}
62// --------------------------------------------------------------------
63
64// ---- Constructor ----------------------------------------------------
65CbmTrdHitRateFastQa::CbmTrdHitRateFastQa(const char* name, const char*)
66 : FairTask(name)
67 , myfile()
68 , nTotalAsics(0)
69 , nTotalOptLinks(0)
70 , trdTotalDataRate(0)
71 , h1DataModule(NULL)
72 , h1OptLinksModule(NULL)
73 , Digicounter(-1)
74 , tFile(NULL)
75 , fDraw(kFALSE)
76 , fPlane(-1)
77 , fStation(-1)
78 , fLayer(-1)
79 , fCol_mean(-1)
80 , fCol_in(-1)
81 , fCol_out(-1)
82 , fRow_mean(-1)
83 , fRow_in(-1)
84 , fRow_out(-1)
85 , fModuleID(-1)
86 , fMCindex(-1)
87 , local_meanLL()
88 , local_meanC()
89 , global_meanLL()
90 , global_meanC()
91 , local_inLL()
92 , local_inC()
93 , global_inLL()
94 , global_inC()
95 , local_outLL()
96 , local_outC()
97 , global_outLL()
98 , global_outC()
99 , fx_in(0.)
100 , fx_out(0.)
101 , fy_in(0.)
102 , fy_out(0.)
103 , fz_in(0.)
104 , fz_out(0.)
105 , fx_mean(0.)
106 , fy_mean(0.)
107 , fz_mean(0.)
108 , fSector(-1)
109 , padsize()
110 , modulesize()
111 , fELoss(-1.)
112 , fELossdEdX(-1.)
113 , fELossTR(-1.)
114 , fPosXLL(0.)
115 , fPosYLL(0.)
116 , fPadPosxLL(0.)
117 , fPadPosyLL(0.)
118 , fPadPosxC(0.)
119 , fPadPosyC(0.)
120 , fDeltax(0.)
121 , fDeltay(0.)
122 , fPadCharge()
123 , fPRFHitPositionLL(0.)
124 , fPRFHitPositionC(0.)
125 , fEfficiency(1.)
126 , fTrdPoints(NULL)
127 , fDigiCollection(NULL)
128 , fDigiMatchCollection(NULL)
129 , fMCStacks(NULL)
130 , fAsicPar(NULL)
131 , fDigiPar(NULL)
132 , fGeoPar(NULL)
133 , fGeoHandler(new CbmTrdGeoHandler())
134 , fColors()
135 , fZLevel()
136 , fBitPerHit(0)
137 , fDigiMap()
138 , fDigiMapIt()
139{
140}
141// --------------------------------------------------------------------
142
143// ---- Destructor ----------------------------------------------------
144CbmTrdHitRateFastQa::~CbmTrdHitRateFastQa()
145{
146 // FairRootManager *ioman =FairRootManager::Instance();
147 //ioman->Write();
148 //fDigiCollection->Clear("C");
149 //delete fDigiCollection;
150 //fDigiMatchCollection->Clear("C");
151 //delete fDigiMatchCollection;
152}
153// --------------------------------------------------------------------
154
155// ---- Initialisation ----------------------------------------------
156void CbmTrdHitRateFastQa::SetParContainers()
157{
158 cout << " * CbmTrdHitRateFastQa * :: SetParContainers() " << endl;
159
160
161 // Get Base Container
162 FairRunAna* ana = FairRunAna::Instance();
163 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
164
165 fAsicPar = (CbmTrdParSetAsic*) (rtdb->getContainer("CbmTrdParSetAsic"));
166 fDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
167 fGeoPar = (CbmTrdParSetGeo*) (rtdb->getContainer("CbmTrdParSetGeo"));
168}
169// --------------------------------------------------------------------
170
171// ---- ReInit -------------------------------------------------------
172InitStatus CbmTrdHitRateFastQa::ReInit()
173{
174
175 cout << " * CbmTrdHitRateFastQa * :: ReInit() " << endl;
176
177
178 FairRunAna* ana = FairRunAna::Instance();
179 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
180
181 fAsicPar = (CbmTrdParSetAsic*) (rtdb->getContainer("CbmTrdParSetAsic"));
182 fDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
183 fGeoPar = (CbmTrdParSetGeo*) (rtdb->getContainer("CbmTrdParSetGeo"));
184
185 return kSUCCESS;
186}
187// --------------------------------------------------------------------
188
189// ---- Init ----------------------------------------------------------
190InitStatus CbmTrdHitRateFastQa::Init()
191{
192
193 cout << " * CbmTrdHitRate * :: Init() " << endl;
194
195 FairRootManager* ioman = FairRootManager::Instance();
196 if (!ioman) Fatal("Init", "No FairRootManager");
197
198 fTrdPoints = (TClonesArray*) ioman->GetObject("TrdPoint");
199 if (!fTrdPoints) {
200 cout << "-W CbmTrdCluster::Init: No TrdPoints array!" << endl;
201 cout << " Task will be inactive" << endl;
202 return kERROR;
203 }
204
205 fMCStacks = (TClonesArray*) ioman->GetObject("MCTrack");
206 if (!fMCStacks) {
207 cout << "-W CbmTrdCluster::Init: No MCTrack array!" << endl;
208 cout << " Task will be inactive" << endl;
209 return kERROR;
210 }
211 /*
212 fDigiCollection = new TClonesArray("CbmTrdDigi", 100);
213 ioman->Register("TrdDigi","TRD Digis",fDigiCollection,IsOutputBranchPersistent("TrdDigi"));
214
215 fDigiMatchCollection = new TClonesArray("CbmTrdDigiMatch", 100);
216 ioman->Register("TrdDigiMatch","TRD Digis",fDigiMatchCollection,IsOutputBranchPersistent("TrdDigiMatch"));
217 */
218 fGeoHandler->Init();
219
220 //fRadiators->Init();
221
222 return kSUCCESS;
223}
224// --------------------------------------------------------------------
225
226
227// ---- Exec ----------------------------------------------------------
228void CbmTrdHitRateFastQa::Exec(Option_t*)
229{
230 myfile.open("hits_per_asic.txt", std::fstream::out);
231 myfile << "#" << endl;
232 myfile << "## TRD data per ASIC" << endl;
233 myfile << "#" << endl;
234 myfile << "#" << endl;
235 myfile << "#--------------------------" << endl;
236
237 fBitPerHit = 112 * 1.5; // (112 + 4*16) * 10 / 8.; // 6 time samples, 8b/10b encoding, CBMnet header
238 // fBitPerHit = 220; // (112 + 4*16) * 10 / 8.; // 6 time samples, 8b/10b encoding, CBMnet header
239 // fBitPerHit = 112; // 6 time samples 3 + (9 * 6 + 3) / 15 = 7 words = 7 * 16 bit = 112 bits
240
241 //TStyle::SetNumberContours(99);
242 gStyle->SetNumberContours(99);
243 //TH2::SetContour(99);
244 /*
245 const Int_t Number = 11;
246 Double_t r[Number] = {1.00, 0.50, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00, 0.50, 1.00};
247 Double_t g[Number] = {0.00, 0.00, 0.00, 0.50, 1.00, 1.00, 1.00, 1.00, 1.00, 0.50, 1.00};
248 Double_t b[Number] = {1.00, 1.00, 1.00, 1.00, 1.00, 0.50, 0.00, 0.00, 0.00, 0.50, 1.00};
249 Double_t length[Number] = {0.00, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00};
250 Int_t nb = 99;
251 Int_t a = TColor::CreateGradientColorTable(Number,length,r,g,b,nb);
252 */
253 printf("Introduction:\n");
254 // HitRateGeoPara2 *GeoPara = new HitRateGeoPara2;
255 // Bool_t Lines;
256 // Bool_t Fast = true; // false; // true; // false; // will not fill the root file (in Histo)!!
257 // Bool_t firstLayer = false;
258 fDraw = true; // false;
259 Double_t ZRangeL = 1;
260 Double_t ZRangeU = 1e05;
261 Double_t mm2bin = 2.5; // 20.0; // 10.0; // 5.0; // 2.5;
262
263 fStation = 0;
264 fLayer = 0;
265
267 TFile* oldFile = gFile;
268 TDirectory* oldDir = gDirectory;
269
270 tFile = new TFile("CbmTrdHitRateFastQa.root", "RECREATE", " ROOT file with histograms");
271
273 gFile = oldFile;
274 gDirectory = oldDir;
275
276 size_t buf_size = 50;
277 Char_t name[buf_size];
278 Char_t title[buf_size];
279
280
281 snprintf(name, buf_size - 1, "HA_S%d_L%d", fStation, fLayer);
282 // sprintf(title,"DataAsic_Station %d, Layer %d",fStation,fLayer);
283 snprintf(title, buf_size - 1, "Data_per_Asic");
284 // TH1F* h1HitAsic = new TH1F(name,title,50*fBitPerHit,1,10*fBitPerHit);
285 // TH1F* h1HitAsic = new TH1F(name,title,1000,1,2000); // Mbit
286
287 // set even bin size on logx scale
288 const Int_t anbins = 100;
289 Double_t axmin = 1; // Mbit
290 Double_t axmax = 2000; // Mbit
291 Double_t alogxmin = TMath::Log10(axmin);
292 Double_t alogxmax = TMath::Log10(axmax);
293 Double_t abinwidth = (alogxmax - alogxmin) / anbins;
294 Double_t axbins[anbins + 1];
295 axbins[0] = axmin;
296 for (Int_t i = 1; i <= anbins; i++) {
297 axbins[i] = axmin + TMath::Power(10, alogxmin + i * abinwidth);
298 }
299 TH1F* h1HitAsic = new TH1F(name, title, anbins, axbins); // Mbit
300
301 if (fBitPerHit == 1.) h1HitAsic->SetXTitle("Hits/Asic [Hz]");
302 else
303 h1HitAsic->SetXTitle("Data/32ch Asic [Mbit/s]");
304 h1HitAsic->SetYTitle("count");
305 // h1HitAsic->GetYaxis()->SetRangeUser(0,20);
306
307 snprintf(name, buf_size - 1, "HM_S%d_L%d", fStation, fLayer);
308 // sprintf(title, buf_size-1, "DataModule_Station %d, Layer %d",fStation,fLayer);
309 snprintf(title, buf_size - 1, "Data_per_Module");
310 // h1DataModule = new TH1F(name,title,50*fBitPerHit,10,100*10*fBitPerHit);
311 // h1DataModule = new TH1F(name,title,1000,10,100*2000);
312
313 // set even bin size on logx scale
314 const Int_t bnbins = 100;
315 Double_t bxmin = 10; // Mbit
316 Double_t bxmax = 100 * 2000; // Mbit
317 Double_t blogxmin = TMath::Log10(bxmin);
318 Double_t blogxmax = TMath::Log10(bxmax);
319 Double_t bbinwidth = (blogxmax - blogxmin) / bnbins;
320 Double_t bxbins[bnbins + 1];
321 bxbins[0] = bxmin;
322 for (Int_t i = 1; i <= bnbins; i++) {
323 bxbins[i] = bxmin + TMath::Power(10, blogxmin + i * bbinwidth);
324 }
325 h1DataModule = new TH1F(name, title, bnbins, bxbins); // Mbit
326
327 if (fBitPerHit == 1.) h1DataModule->SetXTitle("Hits/Module [Hz]");
328 else
329 h1DataModule->SetXTitle("Data/Module [Mbit/s]");
330 h1DataModule->SetYTitle("count");
331 // h1DataModule->GetYaxis()->SetRangeUser(0,20);
332
333 snprintf(name, buf_size - 1, "HO_S%d_L%d", fStation, fLayer);
334 // sprintf(title,"OptLinksModule_Station %d, Layer %d",fStation,fLayer);
335 snprintf(title, buf_size - 1, "5_Gbps_optical_links_per_Module");
336 h1OptLinksModule = new TH1F(name, title, 20, 0.5, 20.5);
337 h1OptLinksModule->SetXTitle("optical links");
338 h1OptLinksModule->SetYTitle("count");
339 // h1OptLinksModule->GetYaxis()->SetRangeUser(0,20);
340
341 TH1F* h1HitPad = NULL;
342 TH2F* h2Layer = NULL;
343 TH2F* h2Topview[3] = {NULL, NULL, NULL};
344 TCanvas* c1 = NULL;
345 TCanvas* c3 = NULL;
346 TCanvas* c2 = NULL;
347 TCanvas* c0 = NULL;
348 if (fDraw) {
349 c0 = new TCanvas("c0", "c0", 1800, 450);
350 c0->Divide(3, 1);
351 }
352 const Int_t z_cave = 12000;
353 h2Topview[0] =
354 new TH2F("TopView", "TopView", int(2 * winsize / mm2bin), -winsize, winsize, int(z_cave / mm2bin), 0, z_cave);
355 h2Topview[0]->SetXTitle("x-Coordinate [mm]");
356 h2Topview[0]->SetYTitle("z-Coordinate [mm]");
357 h2Topview[0]->SetZTitle("#");
358 h2Topview[1] = new TH2F("FrontView", "FrontView", int(2 * winsize / mm2bin), -winsize, winsize,
359 int(2 * winsize / mm2bin), -winsize, winsize);
360 h2Topview[1]->SetXTitle("x-Coordinate [mm]");
361 h2Topview[1]->SetYTitle("y-Coordinate [mm]");
362 h2Topview[1]->SetZTitle("#");
363 h2Topview[2] =
364 new TH2F("SideView", "SideView", int(z_cave / mm2bin), 0, z_cave, int(2 * winsize / mm2bin), -winsize, winsize);
365 h2Topview[2]->SetXTitle("z-Coordinate [mm]");
366 h2Topview[2]->SetYTitle("y-Coordinate [mm]");
367 h2Topview[2]->SetZTitle("#");
368 for (Int_t i = 0; i < 3; i++) {
369 h2Topview[i]->SetContour(99);
370 h2Topview[i]->SetStats(kFALSE);
371 h2Topview[i]->GetXaxis()->SetLabelSize(0.02);
372 h2Topview[i]->GetYaxis()->SetLabelSize(0.02);
373 h2Topview[i]->GetZaxis()->SetLabelSize(0.02);
374 h2Topview[i]->GetXaxis()->SetTitleSize(0.02);
375 h2Topview[i]->GetXaxis()->SetTitleOffset(1.5);
376 h2Topview[i]->GetYaxis()->SetTitleSize(0.02);
377 h2Topview[i]->GetYaxis()->SetTitleOffset(2);
378 h2Topview[i]->GetZaxis()->SetTitleSize(0.02);
379 h2Topview[i]->GetZaxis()->SetTitleOffset(-2);
380 if (fDraw) {
381 c0->cd(i + 1);
382 h2Topview[i]->Draw("colz");
383 }
384 }
385 /*
386 for (Int_t i = 0; i < 8; i++){
387 fColors.push_back(19-i);
388 fZLevel.push_back(ZRangeL + (i+1) * ((ZRangeU - ZRangeL) / 9));
389 printf("%2i %3i %f\n",i,fColors[i],fZLevel[i]);
390 }
391 fColors.push_back(1);
392 fZLevel.push_back(ZRangeL + (9) * ((ZRangeU - ZRangeL) / 9));
393 */
394 for (Int_t i = 0; i < TColor::GetNumberOfColors(); i++) {
395 fColors.push_back(TColor::GetColorPalette(i));
396 fZLevel.push_back(ZRangeL + TMath::Power(10, TMath::Log10(ZRangeU) / TColor::GetNumberOfColors() * i));
397 //fZLevel.push_back(ZRangeL + (i+1) * ((ZRangeU - ZRangeL) / TColor::GetNumberOfColors()));
398 //printf("%2i %3i %f\n",i,TColor::GetColorPalette(i),ZRangeL + (i+1) * ((ZRangeU - ZRangeL) / TColor::GetNumberOfColors()));
399 }
400 //h2Topview->GetZaxis()->SetRangeUser(ZRangeL,ZRangeU);
401 // return;
402
403 vector<int> Plane01;
404 vector<int> Plane02;
405 vector<int> Plane03;
406 vector<int> Plane04;
407 vector<int> Plane05;
408 vector<int> Plane06;
409 vector<int> Plane07;
410 vector<int> Plane08;
411 vector<int> Plane09;
412 vector<int> Plane10;
413
414 Int_t ModuleID;
415 // Int_t Sector = 0;
416 const Int_t NofModules = fDigiPar->GetNrOfModules();
417
418 //const Int_t NofModules = 10; // DE
419 //TH2F *Module[NofModules];
420
421 for (Int_t i = 0; i < NofModules; i++) {
422
423 ModuleID = fDigiPar->GetModuleId(i);
424 // CbmTrdParSetAsic *fModuleAsic = (CbmTrdParSetAsic*)fAsicPar->GetModuleSet(ModuleID);
425 // CbmTrdParModDigi *fModuleInfo = (CbmTrdParModDigi*)fDigiPar->GetModulePar(ModuleID);
426 // CbmTrdParModGeo *fModuleGeo = (CbmTrdParModGeo*)fGeoPar->GetModulePar(ModuleID);
427 //printf("%d:\n %.1f %.1f %.1f\n",ModuleID,fModuleInfo->GetX(), fModuleInfo->GetY(),fModuleInfo->GetZ());
428
429 fPlane = CbmTrdAddress::GetLayerId(ModuleID);
430 fStation = fPlane / 4 + 1; // OK for SIS100 and SIS300
431 fLayer = fPlane % 4 + 1; // OK for SIS100 and SIS300
432
433 //cout << "module " << i+1 << ": ID " << ModuleID << " - P" << fPlane << " S" << fStation << " L" << fLayer << endl;
434 //printf("module %4i: ID $6i - P%2i S%2i L%2i\n",i+1,ModuleID,fPlane,fStation,fLayer);
435 // fill plane vectors
436 if (fPlane == 0) Plane01.push_back(ModuleID);
437 if (fPlane == 1) Plane02.push_back(ModuleID);
438 if (fPlane == 2) Plane03.push_back(ModuleID);
439 if (fPlane == 3) Plane04.push_back(ModuleID);
440 if (fPlane == 4) Plane05.push_back(ModuleID);
441 if (fPlane == 5) Plane06.push_back(ModuleID);
442 if (fPlane == 6) Plane07.push_back(ModuleID);
443 if (fPlane == 7) Plane08.push_back(ModuleID);
444 if (fPlane == 8) Plane09.push_back(ModuleID);
445 if (fPlane == 9) Plane10.push_back(ModuleID);
446 }
447
448 vector<vector<int>> LiSi;
449 //LiSi.push_back(Plane01);
450 //LiSi.push_back(Plane02);
451 //LiSi.push_back(Plane05);
452 //LiSi.push_back(Plane06);
453 //LiSi.push_back(Plane09);
454 //LiSi.push_back(Plane10);
455 LiSi.push_back(Plane01);
456 LiSi.push_back(Plane02);
457 LiSi.push_back(Plane03);
458 LiSi.push_back(Plane04);
459 LiSi.push_back(Plane05);
460 LiSi.push_back(Plane06);
461 LiSi.push_back(Plane07);
462 LiSi.push_back(Plane08);
463 LiSi.push_back(Plane09);
464 LiSi.push_back(Plane10);
465
466 TString OutFile1;
467 TString OutFile2;
468
469 TImage* Outimage1;
470 TImage* Outimage2;
471
472 TLine* MaxHitRatePerPad = new TLine(1e05, 0, 1e05, 1e04); // 100.000 hits per pad
473 MaxHitRatePerPad->SetLineColor(2); // make it red
474 MaxHitRatePerPad->SetLineWidth(8); // make it thick
475
476 /*
477 h2Layer = new TH2F("dummy1","dummy1",1,-0.5,0.5,1,-0.5,0.5);
478 h1HitPad = new TH1F("dummy2","dummy2",1,-0.5,0.5);
479 c1 = new TCanvas("dummy3","dummy3",10,10);
480 c2 = new TCanvas("dummy4","dummy4",10,10);
481 */
482 //**************************************************************************************
483
484 if ((Int_t) LiSi.size() > 0)
485 for (vector<vector<Int_t>>::size_type j = 0; j < LiSi.size(); j++) {
486 // fix layer number
487 UInt_t nModulesInThisLayer = LiSi[j].size();
488 if (nModulesInThisLayer == 0) break;
489 fPlane = CbmTrdAddress::GetLayerId(LiSi[j][0]);
490 fStation = fPlane / 4 + 1; // OK for SIS100 and SIS300
491 fLayer = fPlane % 4 + 1; // OK for SIS100 and SIS300
492
493 printf("Station: %i Layer: %i nModules: %i\n", fStation, fLayer, nModulesInThisLayer);
494
495 h2Topview[0]->Reset();
496 h2Topview[1]->Reset();
497 h2Topview[2]->Reset();
498
499 OutFile1.Form("pics/HitRateLayerPadView_S%d_L%d.png", fStation, fLayer);
500 OutFile2.Form("pics/HitRateLayerSpectrum_S%d_L%d.png", fStation, fLayer);
501
502 HistoInit(c1, c2, c3, h2Layer, h1HitPad, ZRangeL, ZRangeU, mm2bin);
503
504 if (fDraw) {
505 c2->cd(1);
506 MaxHitRatePerPad->Draw("same"); // draw red line
507 }
508 // generate HitPerPadSpectra
509 // generate png files
510 // Lines = false;
511 //c3 = (TCanvas*)c1->Clone("c3");
512 if (nModulesInThisLayer > 0)
513 for (vector<int>::size_type i = 0; i < nModulesInThisLayer; i++) {
514 printf(" ModuleAddress: %i\n", LiSi[j][i]);
515 myfile << "# ModuleAddress: " << LiSi[j][i] << endl;
516 ScanModulePlane(LiSi[j][i], c1, c3, h1HitPad, h1HitAsic);
517 }
518 /*
519 GetModuleInformationFromDigiPar(GeoPara, Fast, Lines, LiSi[j][i], h2Layer ,c1, h1HitPad, c2, h2Topview, c0, mm2bin);
520
521 if(fDraw)
522 {
523 c1->cd(1)->SetLogz(1);
524 h2Layer->Draw("colz"); // draw layer view
525 }
526
527 Lines = true;
528 for (vector<int>::size_type i = 0; i < nModulesInThisLayer; i++)
529 ScanModulePlane();
530 GetModuleInformationFromDigiPar(GeoPara, Fast, Lines, LiSi[j][i], h2Layer ,c1, h1HitPad, c2, h2Topview, c0, mm2bin);
531 */
532 if (fDraw) // dump png file for this layer
533 {
534 Outimage1 = TImage::Create();
535 Outimage1->FromPad(c1);
536 Outimage1->WriteImage(OutFile1);
537 OutFile1.ReplaceAll("png", "pdf");
538 c1->SaveAs(OutFile1);
539 OutFile1.ReplaceAll("Pad", "ASIC");
540 if (c1) delete c1;
541 c3->SaveAs(OutFile1);
542 OutFile1.ReplaceAll("pdf", "png");
543 c3->SaveAs(OutFile1);
544 if (c3) delete c3;
545
546 c2->cd(1);
547 h1HitPad->Draw("same");
548
549 c2->cd(2);
550 h1HitAsic->Draw();
551 h1HitAsic->Write("", TObject::kOverwrite);
552 c2->Update();
553 TLine* MaxDataRatePerUplink = new TLine(500, 0.7 * gPad->GetUymax(), 500,
554 gPad->GetUymax()); // 500 Mbit per Uplink
555 // TLine* MaxDataRatePerUplink = new TLine(500,60,500,80); // 500 Mbit per Uplink
556 MaxDataRatePerUplink->SetLineColor(2); // make it red
557 MaxDataRatePerUplink->SetLineWidth(8); // make it thick
558 MaxDataRatePerUplink->Draw("same"); // draw red line
559
560 c2->cd(3);
561 h1DataModule->Draw();
562 h1DataModule->Write("", TObject::kOverwrite);
563 c2->Update();
564
565 TLine* MaxDataRatePerOptLink = new TLine(5000, 0.7 * gPad->GetUymax(), 5000,
566 gPad->GetUymax()); // 500 Mbit per Uplink
567 // TLine* MaxDataRatePerOptLink = new TLine(5000,5,5000,8); // 500 Mbit per Uplink
568 MaxDataRatePerOptLink->SetLineColor(2); // make it red
569 MaxDataRatePerOptLink->SetLineWidth(8); // make it thick
570 MaxDataRatePerOptLink->Draw("same"); // draw red line
571
572 c2->cd(4);
573 h1OptLinksModule->Draw();
574 h1OptLinksModule->Write("", TObject::kOverwrite);
575
576 Outimage2 = TImage::Create();
577 Outimage2->FromPad(c2);
578 Outimage2->WriteImage(OutFile2);
579 OutFile1.ReplaceAll("png", "pdf");
580 c2->SaveAs(OutFile2);
581 OutFile1.ReplaceAll("Pad", "ASIC");
582 c2->SaveAs(OutFile2);
583
584 h1HitAsic->Reset();
585 h1DataModule->Reset();
586 h1OptLinksModule->Reset();
587
588 if (c2) delete c2;
589 if (Outimage1) delete Outimage1;
590 if (Outimage2) delete Outimage2;
591 }
592 if (h2Layer) delete h2Layer;
593 if (h1HitPad) delete h1HitPad;
594 if (fDraw) c0->Update();
595 }
596 if (fDraw) c0->Update();
597 h1HitAsic->Write("", TObject::kOverwrite);
598
599 printf(" --------------------------\n");
600 printf(" total number of ASICs : %d\n", nTotalAsics);
601 printf(" total number of optical links: %d\n", nTotalOptLinks);
602 printf(" total TRD data rate : %.2f (Gbit/s)\n", trdTotalDataRate * 1e-3);
603 printf(" --------------------------\n");
604 printf(" --------------------------\n");
605
606 myfile << "# total number of ASICs: " << nTotalAsics << endl;
607 myfile << "#--------------------------" << endl;
608 myfile << "#--------------------------" << endl;
609 myfile.close();
610}
611
612
613void CbmTrdHitRateFastQa::ScanModulePlane(const Int_t moduleAddress, TCanvas*& c1, TCanvas*& c2, TH1F*& HitPad,
614 TH1F*& HitAsic)
615{
616 c1->cd(1);
617 Double_t ratePerModule = 0; // sum of data from this module
618 TVector3 padPos;
619 TVector3 padSize;
620 CbmTrdParModAsic* fModuleAsic = (CbmTrdParModAsic*) fAsicPar->GetModulePar(moduleAddress);
621 CbmTrdParModDigi* fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(moduleAddress);
622 CbmTrdParModGeo* fModuleGeo = (CbmTrdParModGeo*) fGeoPar->GetModulePar(moduleAddress);
623 gGeoManager->FindNode(fModuleGeo->GetX(), fModuleGeo->GetY(), fModuleGeo->GetZ());
624 /*gGeoManager->FindNode(fModuleInfo->GetX(),
625 fModuleInfo->GetY(),
626 fModuleInfo->GetZ());
627 printf("%s\n",TString(gGeoManager->GetPath()).Data());
628 */
629 std::vector<Int_t> AsicAddresses;
630 fModuleAsic->GetAsicAddresses(&AsicAddresses);
631 Int_t nofAsics = fModuleAsic->GetNofAsics();
632 printf(" NofAsics: %3i\n", nofAsics);
633 myfile << "# NofAsics : " << nofAsics << endl;
634 myfile << "# moduleAddress / Asic ID / hits per 32ch Asic per second" << endl;
635 nTotalAsics += nofAsics;
636 std::map<Int_t, Double_t> ratePerAsicMap;
637 for (Int_t iAsic = 0; iAsic < nofAsics; iAsic++) {
638 ratePerAsicMap[AsicAddresses[iAsic]] = 0.;
639 }
640 const Int_t nSec = fModuleInfo->GetNofSectors();
641 for (Int_t s = 0; s < nSec; s++) {
642 const Int_t nRow = fModuleInfo->GetNofRowsInSector(s);
643 const Int_t nCol = fModuleInfo->GetNofColumnsInSector(s);
644 for (Int_t r = 0; r < nRow; r++) {
645 for (Int_t c = 0; c < nCol; c++) {
646 fModuleInfo->GetPosition(/*moduleAddress, */ s, c, r, padPos,
647 padSize); // padPos local or global???
648 Int_t channelAddress = CbmTrdAddress::GetAddress(CbmTrdAddress::GetLayerId(moduleAddress),
649 CbmTrdAddress::GetModuleId(moduleAddress), s, r, c);
650 Double_t local_min[3] = {padPos[0] - 0.5 * padSize[0], padPos[1] - 0.5 * padSize[1], padPos[2]};
651 Double_t local_max[3] = {padPos[0] + 0.5 * padSize[0], padPos[1] + 0.5 * padSize[1], padPos[2]};
652 //Double_t master_min[3];
653 //Double_t master_max[3];
654 //gGeoManager->LocalToMaster(local_min, master_min);
655 //gGeoManager->LocalToMaster(local_max, master_max);
656 if (fModuleInfo->GetOrientation() == 1
657 || fModuleInfo->GetOrientation()
658 == 3) { // Pad size is in local coordinate system where position is in global coordinate system
659 local_min[0] = padPos[0] - 0.5 * padSize[1];
660 local_min[1] = padPos[1] - 0.5 * padSize[0];
661 local_max[0] = padPos[0] + 0.5 * padSize[1];
662 local_max[1] = padPos[1] + 0.5 * padSize[0];
663 }
664 //Double_t global_min[3];
665 //Double_t global_max[3];
666 //Double_t x_min(padPos[0]-0.5*padSize[0]), x_max(padPos[0]+0.5*padSize[0]), y_min(padPos[1]-0.5*padSize[1]), y_max(padPos[1]+0.5*padSize[1]);
667 //gGeoManager->LocalToMaster(local_min, global_min);
668 //gGeoManager->LocalToMaster(local_max, global_max);
669 //TBox *pad = new TBox(x_min, y_min, x_max, y_max);
670 for (Int_t i = 0; i < 3; i++) {
671 //master_min[i] *= 10.;
672 //master_max[i] *= 10.;
673 local_min[i] *= 10.;
674 local_max[i] *= 10.;
675 }
676 //Double_t rate = CalcHitRatePad(global_min[0], global_max[0], global_min[1], global_max[1], global_max[2]);
677 Double_t rate = CalcHitRatePad(local_min[0], local_max[0], local_min[1], local_max[1], local_max[2]);
678 HitPad->Fill(rate);
679 ratePerModule += rate; // increase sum by rate of this asic
680
681 //TBox *pad = new TBox(global_min[0], global_min[1], global_max[0], global_max[1]);
682 TBox* pad = new TBox(local_min[0], local_min[1], local_max[0], local_max[1]);
683 //printf(" %i %i %i (%f, %f) (%f, %f) %f\n",s,r,c,local_min[0],local_min[1],global_min[0],global_min[1],rate);
684 pad->SetLineColor(0);
685 pad->SetLineWidth(0);
686 // Int_t color(0), j(0);
687 UInt_t j(0);
688 while ((rate > fZLevel[j]) && (j < fZLevel.size())) {
689 //printf (" %i<%i %i %E <= %E\n",j,(Int_t)fZLevel.size(),fColors[j], rate, fZLevel[j]);
690 j++;
691 }
692 //printf ("%i<%i %i %E <= %E\n\n",j,(Int_t)fZLevel.size(),fColors[j], rate, fZLevel[j]);
693 pad->SetFillColor(fColors[j]);
694 if (j >= fZLevel.size() || rate > 1E05) pad->SetFillColor(12);
695 pad->Draw("same");
696
697 //delete pad;
698 Int_t AsicAddress = fModuleAsic->GetAsicAddress(channelAddress);
699 if (AsicAddress < 0)
700 LOG(error) << "CbmTrdHitRateFastQa::ScanModulePlane: Channel address:" << channelAddress
701 << " is not initialized in module " << moduleAddress << "(s:" << s << ", r:" << r << ", c:" << c
702 << ")";
703 ratePerAsicMap[AsicAddress] += rate;
704 }
705 }
706 //TLine *sec = new TLine();
707 //sec->SetLineStyle(2);
708 //sec->Draw("same");
709 }
710 TBox* module = new TBox(
711 fModuleGeo->GetX() * 10 - fModuleInfo->GetSizeX() * 10, fModuleGeo->GetY() * 10 - fModuleInfo->GetSizeY() * 10,
712 fModuleGeo->GetX() * 10 + fModuleInfo->GetSizeX() * 10, fModuleGeo->GetY() * 10 + fModuleInfo->GetSizeY() * 10);
713 module->SetFillStyle(0);
714 module->Draw("same");
715 c1->Update();
716 c2->cd(1);
717 module->Draw("same");
718 c2->Update();
719 CbmTrdParAsic* Asic = NULL;
720 TBox* asic = NULL;
721 CbmTrdUtils* utils = new CbmTrdUtils();
722 utils->InitColorVector(true, 0., 4e6);
723 for (Int_t iAsic = 0; iAsic < nofAsics; iAsic++) {
724 Asic = fModuleAsic->GetAsicPar(AsicAddresses[iAsic]);
725
726 Double_t local_min[3];
727 Double_t local_max[3];
728 Double_t master_min[3];
729 Double_t master_max[3];
730 local_min[0] = Asic->GetX() - 0.5 * Asic->GetSizeX();
731 local_min[1] = Asic->GetY() - 0.5 * Asic->GetSizeY();
732 local_min[2] = fModuleGeo->GetZ();
733 local_max[0] = Asic->GetX() + 0.5 * Asic->GetSizeX();
734 local_max[1] = Asic->GetY() + 0.5 * Asic->GetSizeY();
735 local_max[2] = fModuleGeo->GetZ();
736 gGeoManager->LocalToMaster(local_min, master_min);
737 gGeoManager->LocalToMaster(local_max, master_max);
738 //printf("moduleAddress:%10i O:%i Asic:%3i rate:%E\n",moduleAddress,fModuleInfo->GetOrientation(),AsicAddresses[iAsic],ratePerAsicMap[AsicAddresses[iAsic]]);
739 asic = new TBox(10. * master_min[0], 10. * master_min[1], 10. * master_max[0], 10. * master_max[1]);
740 c2->cd(1);
741 if (ratePerAsicMap[AsicAddresses[iAsic]] > 4e6) asic->SetFillColor(kBlack);
742 else
743 asic->SetFillColor(utils->GetColorCode(ratePerAsicMap[AsicAddresses[iAsic]]));
744
745 // myfile << iAsic << " " << AsicAddresses[iAsic] << " " << ratePerAsicMap[AsicAddresses[iAsic]] << endl;
746 myfile << moduleAddress << " " << setfill('0') << setw(2) << iAsic << " " << setiosflags(ios::fixed)
747 << setprecision(0) << setfill(' ') << setw(8) << ratePerAsicMap[AsicAddresses[iAsic]] << endl;
748
749 Double_t dataPerAsic = ratePerAsicMap[AsicAddresses[iAsic]] * 3 * 1e-6 * fBitPerHit; // Mbit, incl. neighbor
750 HitAsic->Fill(dataPerAsic);
751 asic->Draw("same");
752 //c2->Update();
753 }
754
755 Double_t dataPerModule = ratePerModule * 3 * 1e-6 * fBitPerHit; // Mbit, incl. neighbor
756 Int_t nOptLinks = 1 + dataPerModule / 4000.; // 5000.; // 1 link plus 1 for each 4 Gbps (fill links to 80% max)
757 h1DataModule->Fill(dataPerModule);
758 h1OptLinksModule->Fill(nOptLinks);
759
760 // global statistics
761 nTotalOptLinks += nOptLinks;
762 trdTotalDataRate += dataPerModule;
763
764 printf(" data rate: %11.4f (Gbit/s)\n", dataPerModule * 1e-3);
765 printf(" opt links: %6i\n", nOptLinks);
766 printf(" --------------------------\n");
767
768 myfile << "#--------------------------" << endl;
769}
770
771
772void CbmTrdHitRateFastQa::HistoInit(TCanvas*& c1, TCanvas*& c2, TCanvas*& c3, TH2F*& Layer, TH1F*& HitPad,
773 Double_t ZRangeL, Double_t ZRangeU, Double_t mm2bin)
774{
775 size_t buf_size = 50;
776 Char_t name[buf_size];
777 Char_t title[buf_size];
778
779
780 snprintf(name, buf_size - 1, "HP_S%d_L%d", fStation, fLayer);
781 snprintf(title, buf_size - 1, "HitPad_Station %d, Layer %d", fStation, fLayer);
782
783 // set even bin size on logx scale
784 const Int_t cnbins = 200;
785 Double_t cxmin = 1e2;
786 Double_t cxmax = 1e6;
787 Double_t clogxmin = TMath::Log10(cxmin);
788 Double_t clogxmax = TMath::Log10(cxmax);
789 Double_t cbinwidth = (clogxmax - clogxmin) / cnbins;
790 Double_t cxbins[cnbins + 1];
791 cxbins[0] = cxmin;
792 for (Int_t i = 1; i <= cnbins; i++) {
793 cxbins[i] = cxmin + TMath::Power(10, clogxmin + i * cbinwidth);
794 }
795 HitPad = new TH1F(name, title, cnbins, cxbins);
796
797 // HitPad = new TH1F(name,title,10000,1e02,1e06);
798 HitPad->SetXTitle("Hits/Pad [Hz]");
799 HitPad->SetYTitle("count");
800 HitPad->GetYaxis()->SetRangeUser(1, 1e04);
801
802 snprintf(name, buf_size - 1, "c2_S%d_L%d", fStation, fLayer);
803 snprintf(title, buf_size - 1, "c2 Station %d, Layer %d", fStation, fLayer);
804 if (fDraw) {
805 c2 = new TCanvas(name, title, 1600, 900);
806 c2->Divide(2, 2); // (3,1);
807 c2->cd(1)->SetLogx(1);
808 c2->cd(1)->SetLogy(1);
809 c2->cd(1)->SetGridx(1);
810 c2->cd(1)->SetGridy(1);
811 HitPad->Draw();
812
813 // h1HitAsic
814 c2->cd(2)->SetLogx(1); // Data per 32ch Asic lin/log scale
815 c2->cd(2)->SetLogy(0); // will be overwritten later
816 c2->cd(2)->SetGridx(1);
817 c2->cd(2)->SetGridy(1);
818
819 // h1DataModule
820 c2->cd(3)->SetLogx(1); // Data per Module lin/log scale
821 c2->cd(3)->SetLogy(0); // will be overwritten later
822 c2->cd(3)->SetGridx(1);
823 c2->cd(3)->SetGridy(1);
824
825 // h1OptLinksModule
826 c2->cd(4)->SetLogx(0); // Data per Module lin/log scale
827 c2->cd(4)->SetLogy(0); // will be overwritten later
828 c2->cd(4)->SetGridx(1);
829 c2->cd(4)->SetGridy(1);
830 }
831
832 snprintf(name, buf_size - 1, "S%d_L%d", fStation, fLayer);
833 snprintf(title, buf_size - 1, "Station %d, Layer %d", fStation, fLayer);
834 printf("%s\n", title);
835
836 Layer =
837 new TH2F(name, title, int(2 * winsize / mm2bin), -winsize, winsize, int(2 * winsize / mm2bin), -winsize, winsize);
838 Layer->SetContour(99);
839 Layer->SetXTitle("x-Coordinate [mm]");
840 Layer->SetYTitle("y-Coordinate [mm]");
841 Layer->SetZTitle("Hits/Pad [Hz]");
842 Layer->SetStats(kFALSE);
843 Layer->GetXaxis()->SetLabelSize(0.02);
844 Layer->GetYaxis()->SetLabelSize(0.02);
845 Layer->GetZaxis()->SetLabelSize(0.02);
846 Layer->GetXaxis()->SetTitleSize(0.02);
847 Layer->GetXaxis()->SetTitleOffset(1.5);
848 Layer->GetYaxis()->SetTitleSize(0.02);
849 Layer->GetYaxis()->SetTitleOffset(2);
850 Layer->GetZaxis()->SetTitleSize(0.02);
851 Layer->GetZaxis()->SetTitleOffset(-2);
852 Layer->GetZaxis()->SetRangeUser(ZRangeL, ZRangeU);
853 Layer->Fill(0., 0., 0.);
854
855 snprintf(name, buf_size - 1, "c1_S%d_L%d", fStation, fLayer);
856 snprintf(title, buf_size - 1, "c1 Station %d, Layer %d", fStation, fLayer);
857 if (fDraw) {
858 c1 = new TCanvas(name, title, 1000, 900);
859 c1->Divide(1, 1);
860 c1->cd(1)->SetLogz(1);
861 Layer->DrawCopy("colz");
862 snprintf(name, buf_size - 1, "c3_S%d_L%d", fStation, fLayer);
863 snprintf(title, buf_size - 1, "c3 Station %d, Layer %d", fStation, fLayer);
864 c3 = new TCanvas(name, title, 1000, 900);
865 c3->Divide(1, 1);
866 c3->cd(1)->SetLogz(1);
867 Layer->GetZaxis()->SetRangeUser(1., 4e6);
868 Layer->SetZTitle("Hits/Asic [Hz]");
869 Layer->DrawCopy("colz");
870 }
871}
872
873// --------------------------------------------------------------------
874
875// ---- FinishTask-----------------------------------------------------
882
883// --------------------------------------------------------------------
884// ----GetModuleInformation ------------------------------------------
885void CbmTrdHitRateFastQa::GetModuleInformation()
886{
887 // Extract the information about station, layer, module type
888 // and cpoy number of the module from the full path to the
889 // node.
890 // The full path is tokenized at the "/" which diide the different
891 // levels of the geometry.
892 // Knowing the nameing scheme of the volumes one gets the required
893 // information with simple string manipulation.
894 // This is probably not the fastes way, but the speed has to be checked.
895 // The methode works only for versions of Root > 5.20.0, before the
896 // class TStringTocken is not implemented
897
898 TString path = gGeoManager->GetPath();
899 cout << "Path: " << path << endl;
900 TStringToken* bla = new TStringToken(path, "/");
901
902 while (bla->NextToken()) {
903 if (bla->Contains("layer")) {
904 TString bla3 = (TString) *bla;
905 Ssiz_t pos = bla3.Last('_');
906 Ssiz_t substringLength = bla3.Length() - pos - 1;
907 TString bla2 = bla3((bla3.Last('_') + 1), substringLength);
908 TString bla1 = bla3(3, 1);
909 fStation = bla1.Atoi();
910 fLayer = bla2.Atoi();
911 }
912 if (bla->Contains("mod")) {
913 TString bla3 = (TString) *bla;
914 Ssiz_t pos = bla3.Last('_');
915 Ssiz_t substringLength = bla3.Length() - pos - 1;
916 TString bla2 = bla3(pos + 1, substringLength);
917 substringLength = pos - 7;
918 TString bla1 = bla3(7, substringLength);
919 break;
920 }
921 }
922}
923
924// --------------------------------------------------------------------
925// ----GetModuleInformationFromDigiPar ------------------------------------------
926void CbmTrdHitRateFastQa::GetModuleInformationFromDigiPar(HitRateGeoPara2* GeoPara, Bool_t Fast, Bool_t Lines,
927 Int_t VolumeID, TH2F* Layer, TCanvas* c1, TH1F* HitPad,
928 TCanvas* c2, TH2F* Topview[3], TCanvas* c0, Double_t mm2bin)
929{
930 // fPos is >0 for x and y and not rotated
931 // origin of the local coordinate system in
932 // the lower left corner of the chamber,
933 // x to the right (small side of the pads), y up
934 //cout << "GetModuleInformationFromDigiPar" << endl;
935
936 CbmTrdParModDigi* fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(VolumeID);
937 CbmTrdParModGeo* fModuleGeo = (CbmTrdParModGeo*) fGeoPar->GetModulePar(VolumeID);
938 if (fModuleInfo != NULL) {
939
940 // fill GeoPara
941 TVector3 padPos;
942 TVector3 padSize;
943
944 GeoPara->nSec = fModuleInfo->GetNofSectors(); // is always 3
945 GeoPara->moduleId = VolumeID;
946 GeoPara->layerId = fLayer;
947 GeoPara->stationId = fStation;
948
949 GeoPara->mPos[0] = fModuleGeo->GetX() * 10;
950 GeoPara->mPos[1] = fModuleGeo->GetY() * 10;
951 GeoPara->mPos[2] = fModuleGeo->GetZ() * 10; // == z(station) ??
952
953 GeoPara->mSize[0] = fModuleInfo->GetSizeX() * 10;
954 GeoPara->mSize[1] = fModuleInfo->GetSizeY() * 10;
955 GeoPara->mSize[2] = 0;
956
957 GeoPara->nCol = 0; // reset total number of columns
958 GeoPara->nRow = 0; // reset total number of rows
959
960 for (Int_t s = 0; s < GeoPara->nSec; s++) // for all (3) sectors
961 {
962 GeoPara->sCol[s] = 0; // reset number of columns in sector
963 GeoPara->sRow[s] = 0; // reset number of rows in sector
964
965 fModuleInfo->GetPosition(/*VolumeID, */ s, 0, 0, padPos, padSize);
966 GeoPara->pSize[s][2] = 0;
967
968 for (Int_t i = 0; i < 2; i++) // for x and y
969 {
970 GeoPara->pSize[s][i] = padSize[i] * 10; // convert to mm
971
972 if (i == 0) // x direction
973 {
974 GeoPara->sSize[s][i] = fModuleInfo->GetSectorSizeX(s) * 10;
975 if (GeoPara->sSize[s][i] < 2 * GeoPara->mSize[i]) // if sector smaller than module
976 {
977 GeoPara->sCol[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]); // calculate number of pads
978 }
979 else // only one sector
980 {
981 if (s == 0)
982 GeoPara->sCol[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]); // calculate number of pads
983 else
984 GeoPara->sCol[s] = 0; // set other sectors 0
985 }
986 GeoPara->nCol += GeoPara->sCol[s]; // sum up columns
987 }
988
989 if (i == 1) // y direction
990 {
991 GeoPara->sSize[s][i] = fModuleInfo->GetSectorSizeY(s) * 10;
992 if (GeoPara->sSize[s][i] < 2 * GeoPara->mSize[i]) // if sector smaller than module
993 {
994 GeoPara->sRow[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]);
995 }
996 else // only one sector
997 {
998 if (s == 0)
999 GeoPara->sRow[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]); // calculate number of pads
1000 else
1001 GeoPara->sRow[s] = 0; // set other sectors 0
1002 }
1003 GeoPara->nRow += GeoPara->sRow[s]; // sum up columns
1004 }
1005 }
1006 }
1007
1008 // swap, if required
1009 if (GeoPara->nRow > GeoPara->nCol) // if mofre rows than columns
1010 {
1011 Int_t itemp;
1012 Double_t dtemp;
1013
1014 itemp = GeoPara->nRow;
1015 GeoPara->nRow = GeoPara->nCol;
1016 GeoPara->nCol = itemp;
1017
1018 for (Int_t s = 0; s < GeoPara->nSec; s++) // for all (3) sectors
1019 {
1020 dtemp = GeoPara->sRow[s];
1021 GeoPara->sRow[s] = GeoPara->sCol[s];
1022 GeoPara->sCol[s] = dtemp;
1023
1024 dtemp = GeoPara->sSize[s][0];
1025 GeoPara->sSize[s][0] = GeoPara->sSize[s][1];
1026 GeoPara->sSize[s][1] = dtemp;
1027
1028 dtemp = GeoPara->pSize[s][0];
1029 GeoPara->pSize[s][0] = GeoPara->pSize[s][1];
1030 GeoPara->pSize[s][1] = dtemp;
1031 }
1032 cout << "swapped x and y" << endl;
1033 }
1034
1035 // get origin
1036 fModuleInfo->GetPosition(/*VolumeID, */ 0, 0, 0, padPos, padSize);
1037 GeoPara->vOrigin[0] = padPos[0] * 10;
1038 GeoPara->vOrigin[1] = padPos[1] * 10;
1039 GeoPara->vOrigin[2] = padPos[2] * 10;
1040
1041 // get col offset
1042 fModuleInfo->GetPosition(/*VolumeID, */ 0, 1, 0, padPos, padSize);
1043 GeoPara->vX[0] = padPos[0] * 10 - GeoPara->vOrigin[0];
1044 GeoPara->vX[1] = padPos[1] * 10 - GeoPara->vOrigin[1];
1045 GeoPara->vX[2] = padPos[2] * 10 - GeoPara->vOrigin[2];
1046
1047 // get row offset
1048 // fModuleInfo->GetPosition(VolumeID, 0, 0, 1, padPos, padSize); // rather take next of 3 sectors
1049 fModuleInfo->GetPosition(/*VolumeID, */ 1, 0, 0, padPos, padSize);
1050 GeoPara->vY[0] = padPos[0] * 10 - GeoPara->vOrigin[0];
1051 GeoPara->vY[1] = padPos[1] * 10 - GeoPara->vOrigin[1];
1052 GeoPara->vY[2] = padPos[2] * 10 - GeoPara->vOrigin[2];
1053
1054 // normal vector
1055 GeoPara->vN[0] = GeoPara->vX[1] * GeoPara->vY[2] - GeoPara->vX[2] * GeoPara->vY[1];
1056 GeoPara->vN[1] = GeoPara->vX[2] * GeoPara->vY[0] - GeoPara->vX[0] * GeoPara->vY[2];
1057 GeoPara->vN[2] = GeoPara->vX[0] * GeoPara->vY[1] - GeoPara->vX[1] * GeoPara->vY[0];
1058
1059 // inclination angle
1060 GeoPara->lambda = (GeoPara->vOrigin[0] * GeoPara->vN[0] + GeoPara->vOrigin[1] * GeoPara->vN[1]
1061 + GeoPara->vOrigin[2] * GeoPara->vN[2]);
1062
1063 /*x-direction*/
1064 GeoPara->cosX = (GeoPara->vX[0] * 1 + GeoPara->vX[1] * 0 + GeoPara->vX[2] * 0)
1065 / sqrt(pow(GeoPara->vX[0], 2) + pow(GeoPara->vX[1], 2) + pow(GeoPara->vX[2], 2));
1066 /*y-direction*/
1067 GeoPara->cosY = (GeoPara->vY[0] * 0 + GeoPara->vY[1] * 1 + GeoPara->vY[2] * 0)
1068 / sqrt(pow(GeoPara->vY[0], 2) + pow(GeoPara->vY[1], 2) + pow(GeoPara->vY[2], 2));
1069
1070 // check and fix cos values (were larger than 1)
1071 // cout << "a cosX " << GeoPara->cosX << " cosY " << GeoPara->cosY << endl;
1072
1073 if (GeoPara->cosX > 1) GeoPara->cosX = 1;
1074 if (GeoPara->cosX < -1) GeoPara->cosX = -1;
1075 if (GeoPara->cosY > 1) GeoPara->cosY = 1;
1076 if (GeoPara->cosY < -1) GeoPara->cosY = -1;
1077
1078 // cout << "b cosX " << GeoPara->cosX << " cosY " << GeoPara->cosY << endl;
1079
1080 if (GeoPara->vX[0] != 0)
1081 GeoPara->stepDirection[0] =
1082 fabs(GeoPara->vX[0]) / (GeoPara->vX[0]); // is the next pad (1,0,0) on the left or right side?
1083 else
1084 GeoPara->stepDirection[0] =
1085 -fabs(GeoPara->vX[1]) / (GeoPara->vX[1]); // is the next pad (1,0,0) on the left or right side?
1086
1087 if (GeoPara->vY[1] != 0)
1088 GeoPara->stepDirection[1] =
1089 fabs(GeoPara->vY[1]) / (GeoPara->vY[1]); // is the next pad (0,1,0) on the upper or lower side?
1090 else
1091 GeoPara->stepDirection[1] =
1092 -fabs(GeoPara->vY[0]) / (GeoPara->vY[0]); // is the next pad (0,1,0) on the upper or lower side?
1093
1094 // transfer values to old variables
1095 Double_t Mpos[3];
1096 Double_t Msize[3];
1097 for (Int_t i = 0; i < 3; i++) // for x, y and z
1098 {
1099 Mpos[i] = GeoPara->mPos[i]; // set module position
1100 Msize[i] = GeoPara->mSize[i]; // set module size
1101 }
1102
1103 cout.setf(ios::fixed);
1104 cout.precision(0);
1105 if (Lines)
1106 cout << "pos " << setw(6) << Mpos[0] << setw(6) << Mpos[1] << setw(6) << Mpos[2] << " size " << setw(6)
1107 << Msize[0] << setw(6) << Msize[1] << setw(6) << Msize[2] << endl;
1108
1109
1110 const Int_t nSec = fModuleInfo->GetNofSectors(); // is always 3
1111 // cout << nSec << " nSec" << endl;
1112 Double_t Ssize[3 * nSec]; // sector size
1113 Double_t Psize[3 * nSec]; // pad size
1114 for (Int_t iSec = 0; iSec < nSec; iSec++) // iSec = Sector
1115 {
1116 Ssize[0 + iSec * nSec] = GeoPara->sSize[iSec][0];
1117 Ssize[1 + iSec * nSec] = GeoPara->sSize[iSec][1];
1118 Ssize[2 + iSec * nSec] = 0;
1119
1120 Psize[0 + iSec * nSec] = GeoPara->pSize[iSec][0];
1121 Psize[1 + iSec * nSec] = GeoPara->pSize[iSec][1];
1122 Psize[2 + iSec * nSec] = 0;
1123 }
1124
1125 if (Lines)
1126 cout << "sec " << setw(10) << Ssize[0] << setw(10) << Ssize[1] << setw(10) << Ssize[3] << setw(10) << Ssize[4]
1127 << setw(10) << Ssize[6] << setw(10) << Ssize[7] << endl;
1128
1129 cout.precision(2);
1130 if (Lines)
1131 cout << "pad " << setw(10) << Psize[0] << setw(10) << Psize[1] << setw(10) << Psize[3] << setw(10) << Psize[4]
1132 << setw(10) << Psize[6] << setw(10) << Psize[7] << endl;
1133
1134 // nCol and nRow
1135 Int_t nCol = GeoPara->nCol; // columns per module - or - pads per row
1136 Int_t nRow = GeoPara->nRow; // rows per module
1137
1138 const Int_t mCol = fModuleInfo->GetNofColumns();
1139 const Int_t mRow = fModuleInfo->GetNofRows();
1140
1141 if (Lines) {
1142 cout << "col0 " << setw(10) << GeoPara->sCol[0] << " row0 " << setw(10) << GeoPara->sRow[0] << endl;
1143 cout << "col1 " << setw(10) << GeoPara->sCol[1] << " row1 " << setw(10) << GeoPara->sRow[1] << endl;
1144 cout << "col2 " << setw(10) << GeoPara->sCol[2] << " row2 " << setw(10) << GeoPara->sRow[2] << endl;
1145 cout << "col " << setw(10) << nCol << " row " << setw(10) << nRow << endl;
1146
1147 cout << "mcol " << setw(10) << mCol << " mrow " << setw(10) << mRow << endl;
1148 cout << "0col " << setw(10) << fModuleInfo->GetNofColumnsInSector(0) << " 0row " << setw(10)
1149 << fModuleInfo->GetNofRowsInSector(0) << endl;
1150 cout << "1col " << setw(10) << fModuleInfo->GetNofColumnsInSector(1) << " 1row " << setw(10)
1151 << fModuleInfo->GetNofRowsInSector(1) << endl;
1152 cout << "2col " << setw(10) << fModuleInfo->GetNofColumnsInSector(2) << " 2row " << setw(10)
1153 << fModuleInfo->GetNofRowsInSector(2) << endl;
1154 }
1155
1156
1158 // Int_t nCol = 0; // columns per module - or - pads per row
1159 // Int_t nRow = 0; // rows per module
1160 //
1161 // for (Int_t i = 0; i < fnSec; i++)
1162 // {
1163 // if (Ssize[0+i*nSec] < 2 * Msize[0] && Ssize[0+i*nSec] > 0) // sector size within module size
1164 // {
1165 // nCol += int(Ssize[0+i*nSec]/Psize[0+i*nSec]);
1166 // }
1167 // else
1168 // {
1169 // nCol = int(Ssize[0+i*nSec]/Psize[0+i*nSec]);
1170 // }
1171 //
1172 // if (Ssize[1+i*nSec] < 2 * Msize[1] && Ssize[1+i*nSec] > 0) // sector size within module size
1173 // {
1174 // nRow += int(Ssize[1+i*nSec]/Psize[1+i*nSec]);
1175 // }
1176 // else
1177 // {
1178 // nRow = int(Ssize[1+i*nSec]/Psize[1+i*nSec]);
1179 // }
1180 // }
1181 //
1182 // if (Lines)
1183 // cout << "nCol " << setw(10) << nCol << " nRow "<< setw(10) << nRow << endl;
1184
1185 if (Lines) {
1186 cout << "vX0 " << setw(10) << GeoPara->vX[0] << endl;
1187 cout << "vX1 " << setw(10) << GeoPara->vX[1] << endl;
1188 cout << "vY0 " << setw(10) << GeoPara->vY[0] << endl;
1189 cout << "vY1 " << setw(10) << GeoPara->vY[1] << endl;
1190 cout << "step " << setw(10) << GeoPara->stepDirection[0] << setw(10) << GeoPara->stepDirection[1] << endl;
1191 cout << "cos " << setw(10) << GeoPara->cosX << setw(10) << GeoPara->cosY << endl;
1192 cout << "v00xy " << setw(10) << GeoPara->mPos[0] - GeoPara->vOrigin[0] << setw(10)
1193 << GeoPara->mPos[1] - GeoPara->vOrigin[1] << endl;
1194
1195 // cout << "------------------------------------------------------" << endl;
1196 }
1197
1198 // angles are relative to y direction
1199 if ((GeoPara->stepDirection[0] == 1) && (GeoPara->stepDirection[1] == 1)) // OK - vert up
1200 GeoPara->rot_angle = 0;
1201 if ((GeoPara->stepDirection[0] == -1) && (GeoPara->stepDirection[1] == 1)) // OK - hori left
1202 GeoPara->rot_angle = 1;
1203 if ((GeoPara->stepDirection[0] == -1) && (GeoPara->stepDirection[1] == -1)) // OK - vert down
1204 GeoPara->rot_angle = 2;
1205 if ((GeoPara->stepDirection[0] == 1) && (GeoPara->stepDirection[1] == -1)) // OK - hori right
1206 GeoPara->rot_angle = 3;
1207
1208 // GeoPara->ori = fModuleInfo->GetOrientation(); // is 0,1,2,3
1209 if (Lines)
1210 cout << "ori " << setw(10) << fModuleInfo->GetOrientation() << " rot " << setw(10) << GeoPara->rot_angle
1211 << endl;
1212
1213 // if (GeoPara->rot_angle == 3)
1214 // {
1215 // fModuleInfo->GetPosition(VolumeID, 2, GeoPara->sCol[0], GeoPara->sRow[2], padPos, padSize);
1216 // if (Lines)
1217 // {
1224 // cout << "v11xy " << setw(10) << padPos[0] * 10 - GeoPara->vOrigin[0]
1225 // << setw(10) << padPos[1] * 10 - GeoPara->vOrigin[1] << endl;
1226 // }
1227 // GeoPara->vOrigin[0] = padPos[0] * 10;
1228 // GeoPara->vOrigin[1] = padPos[1] * 10;
1229 // }
1230
1231 if (Lines) cout << "------------------------------------------------------" << endl;
1232
1233 Topview[0]->Fill(GeoPara->mPos[0], GeoPara->mPos[2]);
1234 Topview[1]->Fill(GeoPara->mPos[0], GeoPara->mPos[1]);
1235 Topview[2]->Fill(GeoPara->mPos[2], GeoPara->mPos[1]);
1236
1237 if (Lines) {
1238 DrawPads(GeoPara, Layer, c1);
1239 DrawBorders(GeoPara, Layer, c1);
1240 }
1241 else {
1242 Histo(GeoPara, Fast, Layer, c1, HitPad, c2, Topview, c0, mm2bin);
1243 }
1244
1245 if (Lines)
1246 for (Int_t s = 0; s < GeoPara->nSec; s++) // for all (3) sectors
1247 {
1248 fModuleInfo->GetPosition(/*VolumeID,*/ s, 0, 0, padPos, padSize);
1249
1250 TPolyMarker* start = new TPolyMarker(1);
1251 start->SetPoint(0, padPos(0) * 10, padPos(1) * 10);
1252 start->SetMarkerStyle(22);
1253 start->SetMarkerSize(.8);
1254
1255 if (fDraw) {
1256 c1->cd(1);
1257 start->Draw("same");
1258 }
1259
1260 fModuleInfo->GetPosition(/*VolumeID,*/ s, GeoPara->sCol[0] - 1, GeoPara->sRow[s] - 1, padPos, padSize);
1261
1262 TPolyMarker* end = new TPolyMarker(1);
1263 end->SetPoint(0, padPos(0) * 10, padPos(1) * 10);
1264 end->SetMarkerStyle(20);
1265 end->SetMarkerSize(.8);
1266
1267 if (fDraw) {
1268 c1->cd(1);
1269 end->Draw("same");
1270 }
1271 }
1272 }
1273
1274 else
1275 printf("fModuleInfo == NULL\n");
1276}
1277Double_t CbmTrdHitRateFastQa::CalcHitRatePad(const Double_t x_min, const Double_t x_max, const Double_t y_min,
1278 const Double_t y_max, const Double_t z)
1279{ //[mm]
1280 // Double_t a[3] = { 7.66582e+00, 6.97712e+00, 6.53780e+00};
1281 // Double_t b[3] = {-2.72375e-03, -1.85168e-03, -1.42673e-03};
1282 Double_t r(0), alpha(0), HitRate(0);
1283 Double_t xStepWidth(1), yStepWidth(1); // [mm]
1284 const Int_t xSteps = Int_t(fabs(x_max - x_min) / xStepWidth) - 1;
1285 const Int_t ySteps = Int_t(fabs(y_max - y_min) / yStepWidth) - 1;
1286 Double_t x(x_min + 0.5), y(y_min + 0.5);
1287 // const Double_t padArea = fabs(x_max - x_min) * fabs(y_max - y_min);
1288 for (Int_t yStep = 0; yStep < ySteps; yStep++) {
1289 x = x_min + 0.5;
1290 for (Int_t xStep = 0; xStep < xSteps; xStep++) {
1291 r = sqrt(pow(x, 2) + pow(y, 2));
1292 alpha = atan(r / z) * 1000.; //[mrad]
1293 HitRate += // fit including errors
1294 exp(4.536355e00 + -8.393716e-03 * alpha) + exp(2.400547e01 + -1.208306e-02 * alpha) / (z * z);
1295 x += xStepWidth;
1296 }
1297 y += yStepWidth;
1298 }
1299 return (HitRate); // [Hz/mm^2]
1300}
1301
1302Double_t CbmTrdHitRateFastQa::CalcHitRate(HitRateGeoPara2* GeoPara, Double_t StartX, Double_t /*StopX*/, Int_t xSteps,
1303 Double_t StartY, Double_t /*StopY*/, Int_t ySteps, Double_t* /*Mpos*/,
1304 TH2F** /*Topview[3]*/, TCanvas* /*c0*/)
1305{
1306 //cout << "CalcHitRate" << endl;
1307 Double_t HitRate = 0; //1. / sqrt( pow( StartX,2) + pow( StartY,2));
1308 Double_t r = 0;
1309 Double_t alpha = 0;
1310 // Int_t counter = 0;
1311 // Double_t a[3] = { 7.66582e+00, 6.97712e+00, 6.53780e+00};
1312 // Double_t b[3] = {-2.72375e-03, -1.85168e-03, -1.42673e-03};
1313
1314 Double_t xStepWidth = GeoPara->cosX;
1315 Double_t yStepWidth = GeoPara->cosY;
1316 Double_t x = StartX + 0.5 * GeoPara->cosX; // don't start on the pad border line
1317 Double_t y = StartY + 0.5 * GeoPara->cosY;
1318
1319 // Double_t xStepWidth = fabs(GeoPara->cosX);
1320 // Double_t yStepWidth = fabs(GeoPara->cosY);
1321 // Double_t x = StartX + 0.5 * fabs(GeoPara->cosX); // don't start on the pad border line
1322 // Double_t y = StartY + 0.5 * fabs(GeoPara->cosY);
1323 Double_t z = (GeoPara->lambda - (x * GeoPara->vN[0] + y * GeoPara->vN[1])) / GeoPara->vN[2];
1324 for (Int_t yStep = 0; yStep < ySteps; yStep++) {
1325 x = StartX + 0.5 * GeoPara->cosX;
1326 // x = StartX + 0.5 * fabs(GeoPara->cosX);
1327 for (Int_t xStep = 0; xStep < xSteps; xStep++) {
1328 z = (GeoPara->lambda - (x * GeoPara->vN[0] + y * GeoPara->vN[1])) / GeoPara->vN[2];
1329 r = sqrt(pow(x / 1.5, 2) + pow(y, 2));
1330 alpha = atan(r / z) * 1000.;
1331 /* //Fit without errors
1332 HitRate +=
1333 exp(4.54156e00 + -8.47377e-03 * alpha) +
1334 exp(2.40005e01 + -1.19541e-02 * alpha) /
1335 (z * z)
1336 ;
1337 */
1338 HitRate += // fit including errors
1339 exp(4.536355e00 + -8.393716e-03 * alpha) + exp(2.400547e01 + -1.208306e-02 * alpha) / (z * z);
1340 /*
1341 Topview[0]->Fill(x,z);
1342 Topview[1]->Fill(x,y);
1343 Topview[2]->Fill(z,y);
1344 */
1345 x += xStepWidth;
1346 }
1347 y += yStepWidth;
1348 }
1349 return (HitRate /*/(counter/100.)*/); /*Convertes Hits/Pad -> Hits/cm² on each Pad*/
1350}
1351
1352
1353void CbmTrdHitRateFastQa::Histo(HitRateGeoPara2* GeoPara, Bool_t Fast, TH2F* h2Layer, TCanvas* /*c1*/, TH1F* h1HitPad,
1354 TCanvas* /*c2*/, TH2F* Topview[3], TCanvas* c0, Double_t mm2bin)
1355{
1356 Double_t ZRangeL = 1e00; //1e05;
1357 Double_t ZRangeU = 1e05; //1e06;
1358 TString name;
1359
1360 // show only current module name
1361 name.Form("ModuleID %5d", GeoPara->moduleId);
1362 // cout << " " << name << "\r" << flush;
1363 cout << " " << name << "\n" << flush;
1364
1365 name.Form("Module%05d", GeoPara->moduleId);
1366 TH2F* h2Module;
1367 if (GeoPara->mSize[0] < 500) // small module
1368 h2Module = new TH2F(name, name, 600 / mm2bin + 1, -300.5, 300.5, 600 / mm2bin + 1, -300.5,
1369 300.5); // 501 x 501 bins
1370 else // large module
1371 h2Module = new TH2F(name, name, 1000 / mm2bin + 1, -500.5, 500.5, 1000 / mm2bin + 1, -500.5,
1372 500.5); // 1001 x 1001 bins
1373 // TH2F *h2Module = new TH2F(name,name,1500/mm2bin+1,-750.5,750.5,1500/mm2bin+1,-750.5,750.5);
1374 h2Module->GetZaxis()->SetRangeUser(ZRangeL, ZRangeU);
1375
1376 name.Form("Module%05dHitPad", GeoPara->moduleId);
1377 TH1F* h1HitPadModule = new TH1F(name, name, 1200, 0, 120000);
1378 // TH1F* h1HitPadModule = new TH1F(name,name,10000,1e00,1e06);
1379
1380 //cout << "Histo" << endl;
1381
1382 // printf(" Xp %6d Xs %6d Yp %6d Ys %6d\n", (Int_t)GeoPara->mPos[0], (Int_t)GeoPara->mSize[0], (Int_t)GeoPara->mPos[1], (Int_t)GeoPara->mSize[1]);
1383 // cout << endl << "Xp " << GeoPara->mPos[0] << " Xs " << GeoPara->mSize[0] << " Yp " << GeoPara->mPos[1] << " Ys " << GeoPara->mSize[1] << endl;
1384
1385 Double_t HiteRate = 0;
1386
1387 Int_t iSecX = 0;
1388 Int_t iSecY = 0;
1389
1390 Double_t planeStartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1391 Double_t planeStopX = planeStartX + GeoPara->pSize[iSecX][0];
1392
1393 Double_t planeStartY = GeoPara->mPos[1] - GeoPara->mSize[1];
1394 Double_t planeStopY = planeStartY + GeoPara->pSize[iSecY][1];
1395
1396 Double_t StartX = GeoPara->vOrigin[0]
1397 + 0.5 * GeoPara->pSize[iSecX][0] * GeoPara->cosX; // vOrigin points to the center of pad (0,0,0)
1398 Double_t StopX = StartX + GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1399
1400 Double_t StartY = GeoPara->vOrigin[1] + 0.5 * GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1401 Double_t StopY = StartY + GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1402
1403 Int_t xSteps = 0;
1404 Int_t ySteps = 0;
1405
1406 // cout << "x " << GeoPara->pSize[iSecX][0] << " y " << GeoPara->pSize[iSecX][1] << endl;
1407
1408 // printf("stX %6d spX %6d stY %6d spY %6d\n", (Int_t)StartX, (Int_t)StopX, (Int_t)StartY, (Int_t)StopY);
1409
1410 //Int_t xStepDirection = GeoPara->vX[0] / fabs(GeoPara->vX[0]); // is the next pad (1,0,0) on the left or right side?
1411 //Int_t yStepDirection = GeoPara->vY[1] / fabs(GeoPara->vY[1]); // is the next pad (0,1,0) on the upper or lower side?
1412
1413 //----------------------------------------------------------------------------------------
1414
1415 //if (GeoPara->rot_angle == 1)
1416 //if (GeoPara->rot_angle != 3)
1417 for (Int_t iR = 0; iR < GeoPara->nRow; iR++) // rows
1418 {
1419 StartX = GeoPara->vOrigin[0] + 0.5 * GeoPara->pSize[iSecX][0] * GeoPara->cosX; // why to the left?
1420 StopX = StartX + GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1421 planeStartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1422 planeStopX = planeStartX + GeoPara->pSize[iSecX][0];
1423
1424 for (Int_t iC = 0; iC < GeoPara->nCol; iC++) // for all columns
1425 {
1426 xSteps = GeoPara->pSize[iSecX][0] - 1; // 1 Step / mm
1427 ySteps = GeoPara->pSize[iSecY][1] - 1;
1428
1429 if (Fast) {
1430 //DE
1431 // if (GeoPara->mPos[0] > -1 && GeoPara->mPos[1] > -1)
1432 if (GeoPara->mPos[0] < 1 && GeoPara->mPos[1] < 1) // copy the bottom left quadrant
1433 {
1434 HiteRate = CalcHitRate(GeoPara, StartX, StopX, xSteps, StartY, StopY, ySteps, GeoPara->mPos, Topview, c0);
1435 h1HitPad->Fill(HiteRate, 4); // weight 4 for each quadrant
1436 }
1437 }
1438
1439 else // slow, for each module
1440 {
1441 HiteRate = CalcHitRate(GeoPara, StartX, StopX, xSteps, StartY, StopY, ySteps, GeoPara->mPos, Topview, c0);
1442 h1HitPad->Fill(HiteRate);
1443 h1HitPadModule->Fill(HiteRate); // single module, goes into root file
1444 }
1445
1446 // fill HitRate into h2Layer histo
1447
1448 Int_t mStepX = Int_t((planeStartX - GeoPara->mPos[0]));
1449 Int_t mStepY = Int_t((planeStartY - GeoPara->mPos[1]));
1450
1451 for (Int_t stepY = int(planeStartY / mm2bin); stepY < int(planeStopY / mm2bin); stepY++) {
1452 // mStepX = Int_t((planeStartX-GeoPara->mPos[0]));
1453
1454 for (Int_t stepX = int(planeStartX / mm2bin); stepX < int(planeStopX / mm2bin); stepX++) {
1455
1456 if (Fast) {
1457 //DE
1458 // if (GeoPara->mPos[0]/*GeoPara->mPos[0]*/ > -1 && GeoPara->mPos[1]/*GeoPara->mPos[1]*/ > -1)
1459 if (GeoPara->mPos[0] /*GeoPara->mPos[0]*/ < 1 && GeoPara->mPos[1] /*GeoPara->mPos[1]*/ < 1) {
1460
1461 h2Layer->SetBinContent(stepX + int(winsize / mm2bin), stepY + int(winsize / mm2bin),
1462 HiteRate); // bin coordiantes
1463 h2Layer->SetBinContent(-1 * stepX + int(winsize / mm2bin), stepY + int(winsize / mm2bin),
1464 HiteRate); // bin coordiantes
1465 h2Layer->SetBinContent(stepX + int(winsize / mm2bin), -1 * stepY + int(winsize / mm2bin),
1466 HiteRate); // bin coordiantes
1467 h2Layer->SetBinContent(-1 * stepX + int(winsize / mm2bin), -1 * stepY + int(winsize / mm2bin),
1468 HiteRate); // bin coordiantes
1469
1470 if (GeoPara->moduleId == 5) {
1471 // cout << "filling 5 " << mStepX << " " << mStepY << " " << HiteRate << endl;
1472 h2Module->Fill(mStepX, mStepY,
1473 HiteRate); // single module, goes into root file
1474 }
1475 }
1476 }
1477
1478 else // slow
1479 {
1480
1481 h2Layer->SetBinContent(stepX + int(winsize / mm2bin), stepY + int(winsize / mm2bin),
1482 HiteRate); // bin coordinates
1483
1484 // cout << "StepX " << stepX << " StepY " << stepY << endl;
1485
1486 if (GeoPara->moduleId == 5) {
1487 cout << "filling 5 " << mStepX << " " << mStepY << " " << HiteRate << endl;
1488 h2Module->Fill(mStepX, mStepY,
1489 HiteRate); // single module, goes into root file
1490 }
1491 // cout << "mStepX " << mStepX << " mStepY " << mStepY << endl;
1492 }
1493
1494 mStepX += mm2bin;
1495 }
1496
1497 mStepY += mm2bin;
1498 }
1499
1500 //printf("Sx Sy (%d,%d) nC nR (%d,%d) iC iR (%d,%d) SC SR (%d,%d) Px Py (%.1f,%.1f) StartX StartY (%.1f,%.1f)\n",iSecX,iSecY,nC,nR,iC,iR,SecCol,SecRow,Psize[0+iSecX*nSec],Psize[1+iSecY*nSec],StartX,StartY);
1501
1502 // if (iC == GeoPara->sCol[iSecX]-1) // why -1 ??
1503 if (iC == GeoPara->sCol[iSecX]) {
1504 // cout << "iC " << iC << " " << GeoPara->nCol << " " << iSecX << endl;
1505 iSecX++;
1506 }
1507 StartX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1508 StopX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1509 planeStartX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0];
1510 planeStopX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0];
1511 // StartX += GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1512 // StopX += GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1513 // planeStartX += GeoPara->pSize[iSecX][0];
1514 // planeStopX += GeoPara->pSize[iSecX][0];
1515 }
1516
1517 iSecX = 0;
1518
1519 // if (iR == GeoPara->sRow[iSecY]-1) // why -1 ??
1520 if (iR == GeoPara->sRow[iSecY]) {
1521 // cout << "iR " << iR << " " << GeoPara->nRow << " " << iSecY << endl;
1522 iSecY++;
1523 }
1524
1525 StartY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1526 StopY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1527 planeStartY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1];
1528 planeStopY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1];
1529 // StartY += GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1530 // StopY += GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1531 // planeStartY += GeoPara->pSize[iSecY][1];
1532 // planeStopY += GeoPara->pSize[iSecY][1];
1533 }
1534
1535 //----------------------------------------------------------------------------------------
1536
1537 for (Int_t i = 0; i < 3; i++)
1538 Topview[i]->Write("", TObject::kOverwrite);
1539
1540 h2Layer->Write("", TObject::kOverwrite);
1541 h1HitPad->Write("", TObject::kOverwrite);
1542
1543 h2Module->Write("", TObject::kOverwrite);
1544 h1HitPadModule->Write("", TObject::kOverwrite);
1545 if (h2Module) delete h2Module;
1546 if (h1HitPadModule) delete h1HitPadModule;
1547}
1548
1549
1550void CbmTrdHitRateFastQa::DrawBorders(HitRateGeoPara2* GeoPara, TH2F* /*Layer*/, TCanvas* c1)
1551{
1552 //----------------------Sector--------------------------------------
1553 Double_t SecXStart = 0.0;
1554 Double_t SecYStart = 0.0;
1555 Double_t SecXStop = 0.0;
1556 Double_t SecYStop = 0.0;
1557
1558 for (Int_t iSec = 0; iSec < GeoPara->nSec - 1; iSec++) // would be enough to iterate up to nSec-1
1559 {
1560
1561 if (GeoPara->sSize[iSec][0] < 2 * GeoPara->mSize[0] && GeoPara->sSize[iSec][0] > 0) // in x direction
1562 {
1563 SecXStart += GeoPara->sSize[iSec][0];
1564 SecXStop = SecXStart;
1565 }
1566 else {
1567 SecXStart = 0.0;
1568 SecXStop = GeoPara->sSize[iSec][0];
1569 }
1570
1571 if (GeoPara->sSize[iSec][1] < 2 * GeoPara->mSize[1] && GeoPara->sSize[iSec][1] > 0) // in x direction
1572 {
1573 SecYStart += GeoPara->sSize[iSec][1];
1574 SecYStop = SecYStart;
1575 }
1576 else {
1577 SecYStart = 0.0;
1578 SecYStop = GeoPara->sSize[iSec][1];
1579 }
1580
1581
1582 TLine* sector =
1583 new TLine(GeoPara->mPos[0] - GeoPara->mSize[0] + SecXStart, GeoPara->mPos[1] - GeoPara->mSize[1] + SecYStart,
1584 GeoPara->mPos[0] - GeoPara->mSize[0] + SecXStop, GeoPara->mPos[1] - GeoPara->mSize[1] + SecYStop);
1585 sector->SetLineColor(15);
1586 sector->SetLineStyle(2);
1587
1588 TPolyMarker* corner = new TPolyMarker(1);
1589 corner->SetPoint(0, GeoPara->mPos[0] - GeoPara->stepDirection[0] * GeoPara->mSize[0],
1590 GeoPara->mPos[1] - GeoPara->stepDirection[1] * GeoPara->mSize[1]);
1591 corner->SetMarkerStyle(21);
1592 corner->SetMarkerSize(.8);
1593
1594 if (fDraw) {
1595 c1->cd(1);
1596 sector->Draw("same");
1597 corner->Draw("same");
1598 }
1599 }
1600
1601
1602 //----------------------Module--------------------------------------
1603 TBox* module = new TBox(GeoPara->mPos[0] - GeoPara->mSize[0], GeoPara->mPos[1] - GeoPara->mSize[1],
1604 GeoPara->mPos[0] + GeoPara->mSize[0], GeoPara->mPos[1] + GeoPara->mSize[1]);
1605 module->SetFillStyle(0);
1606 module->SetLineColor(1);
1607
1608 TBox* M_inner =
1609 new TBox(GeoPara->mPos[0] - 100, GeoPara->mPos[1] - 100, GeoPara->mPos[0] + 100, GeoPara->mPos[1] + 100);
1610 // M_inner->SetUniqueID(ModuleId);
1611 M_inner->SetFillColor(kWhite);
1612
1613 if (fDraw) {
1614 c1->cd(1);
1615 module->Draw("same"); // black module frame
1616 // M_inner->Draw("same"); // white box in module center // do not draw fot the time being
1617 c1->Write("", TObject::kOverwrite);
1618 }
1619}
1620// --------------------------------------------------------------------
1621
1622
1623void CbmTrdHitRateFastQa::DrawPads(HitRateGeoPara2* GeoPara, TH2F* /*Layer*/, TCanvas* c1)
1624{
1625 //----------------------Pad--------------------------------------
1626 // const Int_t nR = GeoPara->nRow;
1627 // const Int_t nC = GeoPara->nCol;
1628
1629 Int_t iSecX = 0;
1630 Int_t iSecY = 0;
1631
1632 Double_t StartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1633 Double_t StartY = GeoPara->mPos[1] - GeoPara->mSize[1];
1634
1635 Int_t SecCol = GeoPara->sCol[iSecX];
1636 Int_t SecRow = GeoPara->sRow[iSecY];
1637
1638 for (Int_t iR = 0; iR < GeoPara->nRow; iR++) // rows
1639 {
1640 StartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1641 for (Int_t iC = 0; iC < GeoPara->nCol; iC++) // cols
1642 {
1643 TLine* Pa = new TLine(StartX, // - low
1644 StartY, StartX + GeoPara->pSize[iSecX][0], StartY);
1645 TLine* Pb = new TLine(StartX, // - high
1646 StartY + GeoPara->pSize[iSecY][1], StartX + GeoPara->pSize[iSecX][0],
1647 StartY + GeoPara->pSize[iSecY][1]);
1648 TLine* Pc = new TLine(StartX, // | left
1649 StartY, StartX, StartY + GeoPara->pSize[iSecY][1]);
1650 TLine* Pd = new TLine(StartX + GeoPara->pSize[iSecX][0], // | right
1651 StartY, StartX + GeoPara->pSize[iSecX][0], StartY + GeoPara->pSize[iSecY][1]);
1652 c1->cd(1);
1653 Pa->SetLineColor(17);
1654 //Pa->SetLineStyle(3);
1655 Pa->Draw("same");
1656
1657 Pb->SetLineColor(17);
1658 //Pb->SetLineStyle(3);
1659 Pb->Draw("same");
1660
1661 Pc->SetLineColor(17);
1662 //Pc->SetLineStyle(3);
1663 Pc->Draw("same");
1664
1665 Pd->SetLineColor(17);
1666 //Pd->SetLineStyle(3);
1667 Pd->Draw("same");
1668
1669 //printf("Sx Sy (%d,%d) nC nR (%d,%d) iC iR (%d,%d) SC SR (%d,%d) Px Py (%.1f,%.1f) StartX StartY (%.1f,%.1f)\n",iSecX,iSecY,GeoPara->nCol,GeoPara->nRow,iC,iR,SecCol,SecRow,GeoPara->pSize[iSecX][0],GeoPara->pSize[iSecY][1],StartX,StartY);
1670
1671 if (iC == SecCol - 1) {
1672 iSecX++;
1673 SecCol += GeoPara->sCol[iSecX];
1674 }
1675 StartX += GeoPara->pSize[iSecX][0];
1676 }
1677 iSecX = 0;
1678 if (iR == SecRow - 1) {
1679 iSecY++;
1680 SecRow += GeoPara->sRow[iSecY];
1681 }
1682 StartY += GeoPara->pSize[iSecY][1];
1683 }
1684}
1685// --------------------------------------------------------------------
1686
1687
1688// ------DrawDigi------------------------------------------------------
1689void CbmTrdHitRateFastQa::DrawDigi() {}
1690// --------------------------------------------------------------------
1691
1692
1693// ---- Register ------------------------------------------------------
1694void CbmTrdHitRateFastQa::Register()
1695{
1696 //FairRootManager::Instance()->Register("TrdDigi","Trd Digi", fDigiCollection, IsOutputBranchPersistent("TrdDigi"));
1697 //FairRootManager::Instance()->Register("TrdDigiMatch","Trd Digi Match", fDigiMatchCollection, IsOutputBranchPersistent("TrdDigiMatch"));
1698}
1699// --------------------------------------------------------------------
1700
1701ClassImp(CbmTrdHitRateFastQa)
ClassImp
ClassImp(CbmConverterManager)
CbmTrdAddress.h
Helper class to convert unique channel ID back and forth.
CbmTrdGeoHandler.h
Helper class to extract information from the GeoManager.
winsize
#define winsize
Definition CbmTrdHitRateFastQa.cxx:46
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
exp
friend fvec exp(const fvec &a)
Definition KfSimdPseudo.h:121
CbmTrdAddress::GetModuleId
static uint32_t GetModuleId(uint32_t address)
Return module ID from address.
Definition CbmTrdAddress.h:79
CbmTrdAddress::GetLayerId
static uint32_t GetLayerId(uint32_t address)
Return layer ID from address.
Definition CbmTrdAddress.h:69
CbmTrdAddress::GetAddress
static uint32_t GetAddress(int32_t layerId, int32_t moduleId, int32_t sectorId, int32_t rowId, int32_t columnId)
Return address from system ID, layer, module, sector, column and row IDs.
Definition CbmTrdAddress.h:44
CbmTrdGeoHandler::Init
void Init(Bool_t isSimulation=kFALSE)
Definition CbmTrdGeoHandler.cxx:52
CbmTrdHitRateFastQa::CalcHitRate
Double_t CalcHitRate(HitRateGeoPara2 *GeoPara, Double_t StartX, Double_t StopX, Int_t xSteps, Double_t StartY, Double_t StopY, Int_t ySteps, Double_t *Mpos, TH2F *Topview[3], TCanvas *c0)
Definition CbmTrdHitRateFastQa.cxx:1302
CbmTrdHitRateFastQa::fAsicPar
CbmTrdParSetAsic * fAsicPar
MC Track information.
Definition CbmTrdHitRateFastQa.h:206
CbmTrdHitRateFastQa::Init
virtual InitStatus Init()
Definition CbmTrdHitRateFastQa.cxx:190
CbmTrdHitRateFastQa::fMCStacks
TClonesArray * fMCStacks
Corresponding MCPoints to TRD digis.
Definition CbmTrdHitRateFastQa.h:204
CbmTrdHitRateFastQa::CalcHitRatePad
Double_t CalcHitRatePad(const Double_t x_min, const Double_t x_max, const Double_t y_min, const Double_t y_max, const Double_t z)
Definition CbmTrdHitRateFastQa.cxx:1277
CbmTrdHitRateFastQa::fColors
std::vector< Int_t > fColors
Definition CbmTrdHitRateFastQa.h:212
CbmTrdHitRateFastQa::fDigiPar
CbmTrdParSetDigi * fDigiPar
Definition CbmTrdHitRateFastQa.h:207
CbmTrdHitRateFastQa::GetModuleInformationFromDigiPar
void GetModuleInformationFromDigiPar(HitRateGeoPara2 *GeoPara, Bool_t Fast, Bool_t Lines, Int_t VolumeID, TH2F *Layer, TCanvas *c1, TH1F *HitPad, TCanvas *c2, TH2F *Topview[3], TCanvas *c0, Double_t mm2bin)
Definition CbmTrdHitRateFastQa.cxx:926
CbmTrdHitRateFastQa::ScanModulePlane
void ScanModulePlane(const Int_t moduleId, TCanvas *&c1, TCanvas *&c2, TH1F *&HitPad, TH1F *&HitAsic)
Definition CbmTrdHitRateFastQa.cxx:613
CbmTrdHitRateFastQa::HistoInit
void HistoInit(TCanvas *&c1, TCanvas *&c2, TCanvas *&c3, TH2F *&Layer, TH1F *&HitPad, Double_t ZRangeL, Double_t ZRangeU, Double_t mm2bin)
Definition CbmTrdHitRateFastQa.cxx:772
CbmTrdHitRateFastQa::Histo
void Histo(HitRateGeoPara2 *GeoPara, Bool_t Fast, TH2F *Layer, TCanvas *c1, TH1F *HitPad, TCanvas *c2, TH2F *Topview[3], TCanvas *c0, Double_t mm2bin)
Definition CbmTrdHitRateFastQa.cxx:1353
CbmTrdHitRateFastQa::ReInit
virtual InitStatus ReInit()
Definition CbmTrdHitRateFastQa.cxx:172
CbmTrdHitRateFastQa::fDigiCollection
TClonesArray * fDigiCollection
Trd MC points.
Definition CbmTrdHitRateFastQa.h:202
CbmTrdHitRateFastQa::fGeoPar
CbmTrdParSetGeo * fGeoPar
Definition CbmTrdHitRateFastQa.h:208
CbmTrdHitRateFastQa::fGeoHandler
CbmTrdGeoHandler * fGeoHandler
Definition CbmTrdHitRateFastQa.h:210
CbmTrdHitRateFastQa::fDigiMap
std::map< std::pair< Int_t, std::pair< Int_t, Int_t > >, CbmTrdDigi * > fDigiMap
Definition CbmTrdHitRateFastQa.h:217
CbmTrdHitRateFastQa::fDigiMatchCollection
TClonesArray * fDigiMatchCollection
TRD digis.
Definition CbmTrdHitRateFastQa.h:203
CbmTrdHitRateFastQa::DrawBorders
void DrawBorders(HitRateGeoPara2 *GeoPara, TH2F *Layer, TCanvas *c1)
Definition CbmTrdHitRateFastQa.cxx:1550
CbmTrdHitRateFastQa::DrawPads
void DrawPads(HitRateGeoPara2 *GeoPara, TH2F *Layer, TCanvas *c1)
Definition CbmTrdHitRateFastQa.cxx:1623
CbmTrdHitRateFastQa::Exec
virtual void Exec(Option_t *option)
Definition CbmTrdHitRateFastQa.cxx:228
CbmTrdHitRateFastQa::fZLevel
std::vector< Double_t > fZLevel
Definition CbmTrdHitRateFastQa.h:213
CbmTrdParAsic
Definition of ASIC parameters.
Definition CbmTrdParAsic.h:20
CbmTrdParAsic::GetSizeX
virtual Double_t GetSizeX() const =0
CbmTrdParAsic::GetSizeY
virtual Double_t GetSizeY() const =0
CbmTrdParAsic::GetY
virtual Double_t GetY() const
Definition CbmTrdParAsic.h:40
CbmTrdParAsic::GetX
virtual Double_t GetX() const
Definition CbmTrdParAsic.h:39
CbmTrdParModAsic
Describe TRD module ASIC settings (electronic gain, delays, etc)
Definition CbmTrdParModAsic.h:30
CbmTrdParModAsic::GetAsicAddresses
virtual void GetAsicAddresses(std::vector< Int_t > *a) const
Query the ASICs in the module for their DAQ address. It applies to the list of ASICs....
Definition CbmTrdParModAsic.cxx:81
CbmTrdParModAsic::GetNofAsics
virtual int GetNofAsics() const
Returns the DEFAULT number of ASICs for the current module.
Definition CbmTrdParModAsic.cxx:105
CbmTrdParModAsic::GetAsicAddress
virtual Int_t GetAsicAddress(Int_t chAddress) const
Look for the ASIC which operates on a specific channel It applies to the list of ASICs.
Definition CbmTrdParModAsic.cxx:68
CbmTrdParModAsic::GetAsicPar
virtual const CbmTrdParAsic * GetAsicPar(Int_t address) const
Look for the ASIC parameters of a given DAQ id It applies to the list of ASICs.
Definition CbmTrdParModAsic.cxx:91
CbmTrdParModDigi
Definition of chamber gain conversion for one TRD module.
Definition CbmTrdParModDigi.h:19
CbmTrdParModDigi::GetSectorSizeX
Double_t GetSectorSizeX(Int_t i) const
Definition CbmTrdParModDigi.h:36
CbmTrdParModDigi::GetNofSectors
Int_t GetNofSectors() const
Definition CbmTrdParModDigi.h:43
CbmTrdParModDigi::GetPosition
void GetPosition(Int_t sectorId, Int_t columnId, Int_t rowId, TVector3 &padPos, TVector3 &padSize) const
Definition CbmTrdParModDigi.cxx:752
CbmTrdParModDigi::GetSizeY
Double_t GetSizeY() const
Definition CbmTrdParModDigi.h:84
CbmTrdParModDigi::GetNofColumnsInSector
Int_t GetNofColumnsInSector(Int_t i) const
Definition CbmTrdParModDigi.cxx:336
CbmTrdParModDigi::GetOrientation
Int_t GetOrientation() const
Definition CbmTrdParModDigi.h:31
CbmTrdParModDigi::GetNofRows
Int_t GetNofRows() const
Definition CbmTrdParModDigi.cxx:317
CbmTrdParModDigi::GetSectorSizeY
Double_t GetSectorSizeY(Int_t i) const
Definition CbmTrdParModDigi.h:37
CbmTrdParModDigi::GetNofRowsInSector
Int_t GetNofRowsInSector(Int_t i) const
Definition CbmTrdParModDigi.cxx:343
CbmTrdParModDigi::GetNofColumns
Int_t GetNofColumns() const
Definition CbmTrdParModDigi.cxx:298
CbmTrdParModDigi::GetSizeX
Double_t GetSizeX() const
Definition CbmTrdParModDigi.h:83
CbmTrdParModGeo
Definition of geometry for one TRD module.
CbmTrdParModGeo::GetZ
virtual Double_t GetZ() const
Definition CbmTrdParModGeo.h:37
CbmTrdParModGeo::GetY
virtual Double_t GetY() const
Definition CbmTrdParModGeo.h:31
CbmTrdParModGeo::GetX
virtual Double_t GetX() const
Definition CbmTrdParModGeo.h:25
CbmTrdParSetAsic
Describe TRD module ASIC settings (electronic gain, delays, etc)
Definition CbmTrdParSetAsic.h:28
CbmTrdParSet::GetNrOfModules
virtual Int_t GetNrOfModules() const
Definition CbmTrdParSet.h:38
CbmTrdParSet::GetModuleId
virtual Int_t GetModuleId(Int_t i) const
Definition CbmTrdParSet.cxx:37
CbmTrdParSet::GetModulePar
virtual const CbmTrdParMod * GetModulePar(Int_t detId) const
Definition CbmTrdParSet.cxx:49
CbmTrdUtils::GetColorCode
Int_t GetColorCode(Double_t value)
Definition CbmTrdUtils.cxx:56
CbmTrdUtils::InitColorVector
void InitColorVector(Bool_t logScale, Double_t min, Double_t max)
Definition CbmTrdUtils.cxx:43
HitRateGeoPara2::pSize
Double_t pSize[3][3]
Definition CbmTrdHitRateFastQa.h:37
HitRateGeoPara2::sSize
Double_t sSize[3][3]
Definition CbmTrdHitRateFastQa.h:36