CbmRoot
Loading...
Searching...
No Matches
CbmTrdHitRateQa.cxx
Go to the documentation of this file.
1/* Copyright (C) 2010-2021 Institut fuer Kernphysik, Westfaelische Wilhelms-Universitaet Muenster, Muenster
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: David Emschermann, Cyrano Bergmann [committer] */
4
5#include "CbmTrdHitRateQa.h"
6
7#include "CbmTrdAddress.h"
8#include "CbmTrdGeoHandler.h"
9#include "CbmTrdParModDigi.h"
10#include "CbmTrdParModGeo.h"
11#include "CbmTrdParSetAsic.h"
12#include "CbmTrdParSetDigi.h"
13#include "CbmTrdParSetGeo.h"
14#include "CbmTrdRadiator.h"
15
16#include "FairRootManager.h"
17#include "FairRunAna.h"
18#include "FairRuntimeDb.h"
19
20#include "TBox.h"
21#include "TCanvas.h"
22#include "TClonesArray.h"
23#include "TColor.h"
24#include "TGeoManager.h"
25#include "TH2.h"
26#include "TImage.h"
27#include "TLine.h"
28#include "TPRegexp.h"
29#include "TPolyMarker.h"
30#include "TStyle.h"
31#include <TFile.h>
32
33#include <fstream>
34#include <iostream>
35#include <vector>
36
37#include <cmath>
38
39#include "Riostream.h"
40
41#define winsize 2000 // 1500 // 5500 // 6000
42
43using std::cin;
44using std::cout;
45using std::endl;
46using std::flush;
47using std::ios;
48using std::pair;
49using std::setw;
50using std::vector;
51
52// ---- Default constructor -------------------------------------------
53CbmTrdHitRateQa::CbmTrdHitRateQa() : CbmTrdHitRateQa("TrdHitRateQa", "") {}
54// --------------------------------------------------------------------
55
56// ---- Constructor ----------------------------------------------------
57CbmTrdHitRateQa::CbmTrdHitRateQa(const char* name, const char*)
58 : FairTask(name)
59 , Digicounter(-1)
60 , tFile(NULL)
61 , fDraw(kFALSE)
62 , fPlane(-1)
63 , fStation(-1)
64 , fLayer(-1)
65 , fCol_mean(-1)
66 , fCol_in(-1)
67 , fCol_out(-1)
68 , fRow_mean(-1)
69 , fRow_in(-1)
70 , fRow_out(-1)
71 , fModuleID(-1)
72 , fMCindex(-1)
73 , local_meanLL()
74 , local_meanC()
75 , global_meanLL()
76 , global_meanC()
77 , local_inLL()
78 , local_inC()
79 , global_inLL()
80 , global_inC()
81 , local_outLL()
82 , local_outC()
83 , global_outLL()
84 , global_outC()
85 , fx_in(0.)
86 , fx_out(0.)
87 , fy_in(0.)
88 , fy_out(0.)
89 , fz_in(0.)
90 , fz_out(0.)
91 , fx_mean(0.)
92 , fy_mean(0.)
93 , fz_mean(0.)
94 , fSector(-1)
95 , padsize()
96 , modulesize()
97 , fELoss(-1.)
98 , fELossdEdX(-1.)
99 , fELossTR(-1.)
100 , fPosXLL(0.)
101 , fPosYLL(0.)
102 , fPadPosxLL(0.)
103 , fPadPosyLL(0.)
104 , fPadPosxC(0.)
105 , fPadPosyC(0.)
106 , fDeltax(0.)
107 , fDeltay(0.)
108 , fPadCharge()
109 , fPRFHitPositionLL(0.)
110 , fPRFHitPositionC(0.)
111 , fEfficiency(1.)
112 , fTrdPoints(NULL)
113 , fDigiCollection(NULL)
114 , fDigiMatchCollection(NULL)
115 , fMCStacks(NULL)
116 , fAsicPar(NULL)
117 , fDigiPar(NULL)
118 , fGeoPar(NULL)
119 , fGeoHandler(new CbmTrdGeoHandler())
120 , fDigiMap()
121 , fDigiMapIt()
122{
123}
124// --------------------------------------------------------------------
125
126// ---- Destructor ----------------------------------------------------
127CbmTrdHitRateQa::~CbmTrdHitRateQa()
128{
129 // FairRootManager *ioman =FairRootManager::Instance();
130 //ioman->Write();
131 //fDigiCollection->Clear("C");
132 //delete fDigiCollection;
133 //fDigiMatchCollection->Clear("C");
134 //delete fDigiMatchCollection;
135}
136// --------------------------------------------------------------------
137
138// ---- Initialisation ----------------------------------------------
139void CbmTrdHitRateQa::SetParContainers()
140{
141 cout << " * CbmTrdHitRateQa * :: SetParContainers() " << endl;
142
143
144 // Get Base Container
145 FairRunAna* ana = FairRunAna::Instance();
146 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
147
148 fAsicPar = (CbmTrdParSetAsic*) (rtdb->getContainer("CbmTrdParSetAsic"));
149 fDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
150 fGeoPar = (CbmTrdParSetGeo*) (rtdb->getContainer("CbmTrdParSetGeo"));
151}
152// --------------------------------------------------------------------
153
154// ---- ReInit -------------------------------------------------------
155InitStatus CbmTrdHitRateQa::ReInit()
156{
157
158 cout << " * CbmTrdHitRateQa * :: ReInit() " << endl;
159
160
161 FairRunAna* ana = FairRunAna::Instance();
162 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
163
164 fAsicPar = (CbmTrdParSetAsic*) (rtdb->getContainer("CbmTrdParSetAsic"));
165 fDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
166 fGeoPar = (CbmTrdParSetGeo*) (rtdb->getContainer("CbmTrdParSetGeo"));
167
168 return kSUCCESS;
169}
170// --------------------------------------------------------------------
171
172// ---- Init ----------------------------------------------------------
173InitStatus CbmTrdHitRateQa::Init()
174{
175
176 cout << " * CbmTrdHitRate * :: Init() " << endl;
177
178 FairRootManager* ioman = FairRootManager::Instance();
179 if (!ioman) Fatal("Init", "No FairRootManager");
180
181 fTrdPoints = (TClonesArray*) ioman->GetObject("TrdPoint");
182 if (!fTrdPoints) {
183 cout << "-W CbmTrdCluster::Init: No TrdPoints array!" << endl;
184 cout << " Task will be inactive" << endl;
185 return kERROR;
186 }
187
188 fMCStacks = (TClonesArray*) ioman->GetObject("MCTrack");
189 if (!fMCStacks) {
190 cout << "-W CbmTrdCluster::Init: No MCTrack array!" << endl;
191 cout << " Task will be inactive" << endl;
192 return kERROR;
193 }
194 /*
195 fDigiCollection = new TClonesArray("CbmTrdDigi", 100);
196 ioman->Register("TrdDigi","TRD Digis",fDigiCollection,IsOutputBranchPersistent("TrdDigi"));
197
198 fDigiMatchCollection = new TClonesArray("CbmTrdDigiMatch", 100);
199 ioman->Register("TrdDigiMatch","TRD Digis",fDigiMatchCollection,IsOutputBranchPersistent("TrdDigiMatch"));
200 */
201 fGeoHandler->Init();
202
203 //fRadiators->Init();
204
205 return kSUCCESS;
206}
207// --------------------------------------------------------------------
208
209
210// ---- Exec ----------------------------------------------------------
211void CbmTrdHitRateQa::Exec(Option_t*)
212{
213 //TStyle::SetNumberContours(99);
214 gStyle->SetNumberContours(99);
215 //TH2::SetContour(99);
216 /*
217 const Int_t Number = 11;
218 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};
219 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};
220 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};
221 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};
222 Int_t nb = 99;
223 Int_t a = TColor::CreateGradientColorTable(Number,length,r,g,b,nb);
224 */
225 printf("Introduction:\n");
226 HitRateGeoPara* GeoPara = new HitRateGeoPara;
227 Bool_t Lines;
228 Bool_t Fast = true; // false; // true; // false; // will not fill the root file (in Histo)!!
229 // Bool_t firstLayer = false;
230 fDraw = true; // false;
231 Double_t ZRangeL = 1e00; //1e05;
232 Double_t ZRangeU = 1e05; //1e06;
233 Double_t mm2bin = 2.5; // 20.0; // 10.0; // 5.0; // 2.5;
234
235 fStation = 0;
236 fLayer = 0;
237
239 TFile* oldFile = gFile;
240 TDirectory* oldDir = gDirectory;
241
242 tFile = new TFile("CbmTrdHitRateQa.root", "RECREATE", " ROOT file with histograms");
243
245 gFile = oldFile;
246 gDirectory = oldDir;
247
248 TH1F* h1HitPad = NULL;
249 TH2F* h2Layer = NULL;
250 TH2F* h2Topview[3] = {NULL, NULL, NULL};
251 TCanvas* c1 = NULL;
252 TCanvas* c2 = NULL;
253 TCanvas* c0 = NULL;
254 if (fDraw) {
255 c0 = new TCanvas("c0", "c0", 1800, 450);
256 c0->Divide(3, 1);
257 }
258 const Int_t z_cave = 12000;
259 h2Topview[0] =
260 new TH2F("TopView", "TopView", int(2 * winsize / mm2bin), -winsize, winsize, int(z_cave / mm2bin), 0, z_cave);
261 h2Topview[0]->SetXTitle("x-Coordinate [mm]");
262 h2Topview[0]->SetYTitle("z-Coordinate [mm]");
263 h2Topview[0]->SetZTitle("#");
264 h2Topview[1] = new TH2F("FrontView", "FrontView", int(2 * winsize / mm2bin), -winsize, winsize,
265 int(2 * winsize / mm2bin), -winsize, winsize);
266 h2Topview[1]->SetXTitle("x-Coordinate [mm]");
267 h2Topview[1]->SetYTitle("y-Coordinate [mm]");
268 h2Topview[1]->SetZTitle("#");
269 h2Topview[2] =
270 new TH2F("SideView", "SideView", int(z_cave / mm2bin), 0, z_cave, int(2 * winsize / mm2bin), -winsize, winsize);
271 h2Topview[2]->SetXTitle("z-Coordinate [mm]");
272 h2Topview[2]->SetYTitle("y-Coordinate [mm]");
273 h2Topview[2]->SetZTitle("#");
274 for (Int_t i = 0; i < 3; i++) {
275 h2Topview[i]->SetContour(99);
276 h2Topview[i]->SetStats(kFALSE);
277 h2Topview[i]->GetXaxis()->SetLabelSize(0.02);
278 h2Topview[i]->GetYaxis()->SetLabelSize(0.02);
279 h2Topview[i]->GetZaxis()->SetLabelSize(0.02);
280 h2Topview[i]->GetXaxis()->SetTitleSize(0.02);
281 h2Topview[i]->GetXaxis()->SetTitleOffset(1.5);
282 h2Topview[i]->GetYaxis()->SetTitleSize(0.02);
283 h2Topview[i]->GetYaxis()->SetTitleOffset(2);
284 h2Topview[i]->GetZaxis()->SetTitleSize(0.02);
285 h2Topview[i]->GetZaxis()->SetTitleOffset(-2);
286 if (fDraw) {
287 c0->cd(i + 1);
288 h2Topview[i]->Draw("colz");
289 }
290 }
291 for (Int_t i = 0; i < TColor::GetNumberOfColors(); i++)
292 printf("%2i %3i\n", i, TColor::GetColorPalette(i));
293 //h2Topview->GetZaxis()->SetRangeUser(ZRangeL,ZRangeU);
294 // return;
295
296 vector<int> Plane01;
297 vector<int> Plane02;
298 vector<int> Plane03;
299 vector<int> Plane04;
300 vector<int> Plane05;
301 vector<int> Plane06;
302 vector<int> Plane07;
303 vector<int> Plane08;
304 vector<int> Plane09;
305 vector<int> Plane10;
306
307 Int_t ModuleID;
308 // Int_t Sector = 0;
309 const Int_t NofModules = fDigiPar->GetNrOfModules();
310
311 //const Int_t NofModules = 80; // DE
312 //TH2F *Module[NofModules];
313
314 for (Int_t i = 0; i < NofModules; i++) {
315
316 ModuleID = fDigiPar->GetModuleId(i);
317 // CbmTrdParModDigi *fModuleInfo = (CbmTrdParModDigi*)fDigiPar->GetModulePar(ModuleID);
318 //printf("%d:\n %.1f %.1f %.1f\n",ModuleID,fModuleInfo->GetX(), fModuleInfo->GetY(),fModuleInfo->GetZ());
319
320 fPlane = CbmTrdAddress::GetLayerId(ModuleID);
321 fStation = fPlane / 4 + 1; // OK for SIS100 and SIS300
322 fLayer = fPlane % 4 + 1; // OK for SIS100 and SIS300
323
324 cout << "module " << i + 1 << ": ID " << ModuleID << " - P" << fPlane << " S" << fStation << " L" << fLayer << endl;
325
326 // fill plane vectors
327 if (fPlane == 0) Plane01.push_back(ModuleID);
328 if (fPlane == 1) Plane02.push_back(ModuleID);
329 if (fPlane == 2) Plane03.push_back(ModuleID);
330 if (fPlane == 3) Plane04.push_back(ModuleID);
331 if (fPlane == 4) Plane05.push_back(ModuleID);
332 if (fPlane == 5) Plane06.push_back(ModuleID);
333 if (fPlane == 6) Plane07.push_back(ModuleID);
334 if (fPlane == 7) Plane08.push_back(ModuleID);
335 if (fPlane == 8) Plane09.push_back(ModuleID);
336 if (fPlane == 9) Plane10.push_back(ModuleID);
337 }
338
339 vector<vector<int>> LiSi;
340 //LiSi.push_back(Plane01);
341 //LiSi.push_back(Plane02);
342 //LiSi.push_back(Plane05);
343 //LiSi.push_back(Plane06);
344 //LiSi.push_back(Plane09);
345 //LiSi.push_back(Plane10);
346 LiSi.push_back(Plane01);
347 LiSi.push_back(Plane02);
348 LiSi.push_back(Plane03);
349 LiSi.push_back(Plane04);
350 LiSi.push_back(Plane05);
351 LiSi.push_back(Plane06);
352 LiSi.push_back(Plane07);
353 LiSi.push_back(Plane08);
354 LiSi.push_back(Plane09);
355 LiSi.push_back(Plane10);
356
357 size_t buf_size = 200;
358 Char_t OutFile1[buf_size];
359 Char_t OutFile2[buf_size];
360
361 TImage* Outimage1 = NULL;
362 TImage* Outimage2 = NULL;
363
364 TLine* MaxHitRatePerPad = new TLine(1e05, 0, 1e05, 1e04); // 100.000 hits per pad
365 MaxHitRatePerPad->SetLineColor(2); // make it red
366 MaxHitRatePerPad->SetLineWidth(8); // make it thick
367
368 /*
369 h2Layer = new TH2F("dummy1","dummy1",1,-0.5,0.5,1,-0.5,0.5);
370 h1HitPad = new TH1F("dummy2","dummy2",1,-0.5,0.5);
371 c1 = new TCanvas("dummy3","dummy3",10,10);
372 c2 = new TCanvas("dummy4","dummy4",10,10);
373 */
374 //**************************************************************************************
375
376 for (vector<vector<int>>::size_type j = 0; j < LiSi.size(); j++) {
377 // fix layer number
378 fPlane = CbmTrdAddress::GetLayerId(LiSi[j][0]);
379 fStation = fPlane / 4 + 1; // OK for SIS100 and SIS300
380 fLayer = fPlane % 4 + 1; // OK for SIS100 and SIS300
381 // cout << fStation << ", " << fLayer << endl;
382
383 UInt_t nModulesInThisLayer = LiSi[j].size();
384 // cout << nModulesInThisLayer << " modules in this layer" << endl;
385
386 h2Topview[0]->Reset();
387 h2Topview[1]->Reset();
388 h2Topview[2]->Reset();
389
390 snprintf(OutFile1, buf_size - 1, "pics/HitRateLayerPadView_S%d_L%d.png", fStation, fLayer);
391 snprintf(OutFile2, buf_size - 1, "pics/HitRateLayerSpectrum_S%d_L%d.png", fStation, fLayer);
392
393 HistoInit(c1, c2, h2Layer, h1HitPad, ZRangeL, ZRangeU, mm2bin);
394
395 if (fDraw) MaxHitRatePerPad->Draw("same"); // draw red line
396
397 // generate HitPerPadSpectra
398 // generate png files
399 Lines = false;
400 for (vector<int>::size_type i = 0; i < nModulesInThisLayer; i++)
401 GetModuleInformationFromDigiPar(GeoPara, Fast, Lines, LiSi[j][i], h2Layer, c1, h1HitPad, c2, h2Topview, c0,
402 mm2bin);
403
404 if (fDraw) {
405 c1->cd(1)->SetLogz(1);
406 h2Layer->Draw("colz"); // draw layer view
407 }
408
409 Lines = true;
410 for (vector<int>::size_type i = 0; i < nModulesInThisLayer; i++)
411 GetModuleInformationFromDigiPar(GeoPara, Fast, Lines, LiSi[j][i], h2Layer, c1, h1HitPad, c2, h2Topview, c0,
412 mm2bin);
413
414 if (fDraw) // dump png file for this layer
415 {
416 Outimage1 = TImage::Create();
417 Outimage1->FromPad(c1);
418 Outimage1->WriteImage(OutFile1);
419 }
420
421 if (fDraw) {
422 delete c1;
423 delete Outimage1;
424 h1HitPad->Draw("same");
425 Outimage2 = TImage::Create();
426 Outimage2->FromPad(c2);
427 Outimage2->WriteImage(OutFile2);
428 }
429
430 /*
431 snprintf(OutFile1,buf_size-1,"pics/%s_S%d_L%d.eps",trddigiparpath,fStation,fLayer);
432 c1->cd(1)->Print(OutFile1);
433 */
434 delete h2Layer;
435
436 /*
437 snprintf(OutFile2,buf_size-1,"pics/%s_HitPerPad_S%d_L%d.eps",trddigiparpath,fStation,fLayer);
438 c2->cd(1)->Print(OutFile2);
439 */
440 delete h1HitPad;
441
442 if (fDraw) {
443 delete c2;
444 delete Outimage2;
445 }
446
447 if (fDraw) c0->Update();
448 }
449
450 if (fDraw) c0->Update();
451
452 // // write modules histos to file
453 // for (Int_t i = 0; i < NofModules; i++)
454 // Module[i]->Write("", TObject::kOverwrite);
455
456 // tFile->Close(); // this causes a segfalut
457}
458
459
460void CbmTrdHitRateQa::HistoInit(TCanvas*& c1, TCanvas*& c2, TH2F*& Layer, TH1F*& HitPad, Double_t ZRangeL,
461 Double_t ZRangeU, Double_t mm2bin)
462{
463 size_t buf_size = 50;
464 Char_t name[buf_size];
465 Char_t title[buf_size];
466
467 snprintf(name, buf_size - 1, "S%d_L%d", fStation, fLayer);
468 snprintf(title, buf_size - 1, "Station %d, Layer %d", fStation, fLayer);
469 printf("%s\n", title);
470
471 Layer =
472 new TH2F(name, title, int(2 * winsize / mm2bin), -winsize, winsize, int(2 * winsize / mm2bin), -winsize, winsize);
473 Layer->SetContour(99);
474 Layer->SetXTitle("x-Coordinate [mm]");
475 Layer->SetYTitle("y-Coordinate [mm]");
476 Layer->SetZTitle("Hits/Pad [Hz]");
477 Layer->SetStats(kFALSE);
478 Layer->GetXaxis()->SetLabelSize(0.02);
479 Layer->GetYaxis()->SetLabelSize(0.02);
480 Layer->GetZaxis()->SetLabelSize(0.02);
481 Layer->GetXaxis()->SetTitleSize(0.02);
482 Layer->GetXaxis()->SetTitleOffset(1.5);
483 Layer->GetYaxis()->SetTitleSize(0.02);
484 Layer->GetYaxis()->SetTitleOffset(2);
485 Layer->GetZaxis()->SetTitleSize(0.02);
486 Layer->GetZaxis()->SetTitleOffset(-2);
487 Layer->GetZaxis()->SetRangeUser(ZRangeL, ZRangeU);
488
489 snprintf(name, buf_size - 1, "HP_S%d_L%d", fStation, fLayer);
490 snprintf(title, buf_size - 1, "HitPad_Station %d, Layer %d", fStation, fLayer);
491 HitPad = new TH1F(name, title, 10000, 1e00, 1e06);
492 HitPad->SetXTitle("Hits/Pad [Hz]");
493 HitPad->SetYTitle("count");
494 HitPad->GetYaxis()->SetRangeUser(1, 1e04);
495
496 snprintf(name, buf_size - 1, "c1_S%d_L%d", fStation, fLayer);
497 snprintf(title, buf_size - 1, "c1 Station %d, Layer %d", fStation, fLayer);
498 if (fDraw) {
499 c1 = new TCanvas(name, title, 1000, 900);
500 c1->Divide(1, 1);
501 c1->cd(1)->SetLogz(1);
502 Layer->Draw();
503 }
504 snprintf(name, buf_size - 1, "c2_S%d_L%d", fStation, fLayer);
505 snprintf(title, buf_size - 1, "c2 Station %d, Layer %d", fStation, fLayer);
506 if (fDraw) {
507 c2 = new TCanvas(name, title, 1000, 900 / 2);
508 c2->Divide(1, 1);
509 c2->cd(1)->SetLogx(1);
510 c2->cd(1)->SetLogy(1);
511 c2->cd(1)->SetGridx(1);
512 c2->cd(1)->SetGridy(1);
513 HitPad->Draw();
514 }
515}
516
517// --------------------------------------------------------------------
518
519// ---- FinishTask-----------------------------------------------------
520void CbmTrdHitRateQa::FinishEvent()
521{
522 fDigiMap.clear();
523 if (fDigiCollection) fDigiCollection->Clear();
524 if (fDigiMatchCollection) fDigiMatchCollection->Clear();
525}
526
527// --------------------------------------------------------------------
528// ----GetModuleInformation ------------------------------------------
529void CbmTrdHitRateQa::GetModuleInformation()
530{
531 // Extract the information about station, layer, module type
532 // and cpoy number of the module from the full path to the
533 // node.
534 // The full path is tokenized at the "/" which diide the different
535 // levels of the geometry.
536 // Knowing the nameing scheme of the volumes one gets the required
537 // information with simple string manipulation.
538 // This is probably not the fastes way, but the speed has to be checked.
539 // The methode works only for versions of Root > 5.20.0, before the
540 // class TStringTocken is not implemented
541
542 TString path = gGeoManager->GetPath();
543 cout << "Path: " << path << endl;
544 TStringToken* bla = new TStringToken(path, "/");
545
546 while (bla->NextToken()) {
547 if (bla->Contains("layer")) {
548 TString bla3 = (TString) *bla;
549 Ssiz_t pos = bla3.Last('_');
550 Ssiz_t substringLength = bla3.Length() - pos - 1;
551 TString bla2 = bla3((bla3.Last('_') + 1), substringLength);
552 TString bla1 = bla3(3, 1);
553 fStation = bla1.Atoi();
554 fLayer = bla2.Atoi();
555 }
556 if (bla->Contains("mod")) {
557 TString bla3 = (TString) *bla;
558 Ssiz_t pos = bla3.Last('_');
559 Ssiz_t substringLength = bla3.Length() - pos - 1;
560 TString bla2 = bla3(pos + 1, substringLength);
561 substringLength = pos - 7;
562 TString bla1 = bla3(7, substringLength);
563 break;
564 }
565 }
566}
567
568// --------------------------------------------------------------------
569// ----GetModuleInformationFromDigiPar ------------------------------------------
570void CbmTrdHitRateQa::GetModuleInformationFromDigiPar(HitRateGeoPara* GeoPara, Bool_t Fast, Bool_t Lines,
571 Int_t VolumeID, TH2F* Layer, TCanvas* c1, TH1F* HitPad,
572 TCanvas* c2, TH2F* Topview[3], TCanvas* c0, Double_t mm2bin)
573{
574 // fPos is >0 for x and y and not rotated
575 // origin of the local coordinate system in
576 // the lower left corner of the chamber,
577 // x to the right (small side of the pads), y up
578 //cout << "GetModuleInformationFromDigiPar" << endl;
579
580 CbmTrdParModDigi* fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(VolumeID);
581 CbmTrdParModGeo* fModuleGeo = (CbmTrdParModGeo*) fGeoPar->GetModulePar(VolumeID);
582 if (fModuleInfo != NULL) {
583
584 // fill GeoPara
585 TVector3 padPos;
586 TVector3 padSize;
587
588 GeoPara->nSec = fModuleInfo->GetNofSectors(); // is always 3
589 GeoPara->moduleId = VolumeID;
590 GeoPara->layerId = fLayer;
591 GeoPara->stationId = fStation;
592
593 GeoPara->mPos[0] = fModuleGeo->GetX() * 10;
594 GeoPara->mPos[1] = fModuleGeo->GetY() * 10;
595 GeoPara->mPos[2] = fModuleGeo->GetZ() * 10; // == z(station) ??
596
597 GeoPara->mSize[0] = fModuleInfo->GetSizeX() * 10;
598 GeoPara->mSize[1] = fModuleInfo->GetSizeY() * 10;
599 GeoPara->mSize[2] = 0;
600
601 GeoPara->nCol = 0; // reset total number of columns
602 GeoPara->nRow = 0; // reset total number of rows
603
604 for (Int_t s = 0; s < GeoPara->nSec; s++) // for all (3) sectors
605 {
606 GeoPara->sCol[s] = 0; // reset number of columns in sector
607 GeoPara->sRow[s] = 0; // reset number of rows in sector
608
609 fModuleInfo->GetPosition(/*VolumeID, */ s, 0, 0, padPos, padSize);
610 GeoPara->pSize[s][2] = 0;
611
612 for (Int_t i = 0; i < 2; i++) // for x and y
613 {
614 GeoPara->pSize[s][i] = padSize[i] * 10; // convert to mm
615
616 if (i == 0) // x direction
617 {
618 GeoPara->sSize[s][i] = fModuleInfo->GetSectorSizeX(s) * 10;
619 if (GeoPara->sSize[s][i] < 2 * GeoPara->mSize[i]) // if sector smaller than module
620 {
621 GeoPara->sCol[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]); // calculate number of pads
622 }
623 else // only one sector
624 {
625 if (s == 0)
626 GeoPara->sCol[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]); // calculate number of pads
627 else
628 GeoPara->sCol[s] = 0; // set other sectors 0
629 }
630 GeoPara->nCol += GeoPara->sCol[s]; // sum up columns
631 }
632
633 if (i == 1) // y direction
634 {
635 GeoPara->sSize[s][i] = fModuleInfo->GetSectorSizeY(s) * 10;
636 if (GeoPara->sSize[s][i] < 2 * GeoPara->mSize[i]) // if sector smaller than module
637 {
638 GeoPara->sRow[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]);
639 }
640 else // only one sector
641 {
642 if (s == 0)
643 GeoPara->sRow[s] = round(GeoPara->sSize[s][i] / GeoPara->pSize[s][i]); // calculate number of pads
644 else
645 GeoPara->sRow[s] = 0; // set other sectors 0
646 }
647 GeoPara->nRow += GeoPara->sRow[s]; // sum up columns
648 }
649 }
650 }
651
652 // swap, if required
653 if (GeoPara->nRow > GeoPara->nCol) // if mofre rows than columns
654 {
655 Int_t itemp;
656 Double_t dtemp;
657
658 itemp = GeoPara->nRow;
659 GeoPara->nRow = GeoPara->nCol;
660 GeoPara->nCol = itemp;
661
662 for (Int_t s = 0; s < GeoPara->nSec; s++) // for all (3) sectors
663 {
664 dtemp = GeoPara->sRow[s];
665 GeoPara->sRow[s] = GeoPara->sCol[s];
666 GeoPara->sCol[s] = dtemp;
667
668 dtemp = GeoPara->sSize[s][0];
669 GeoPara->sSize[s][0] = GeoPara->sSize[s][1];
670 GeoPara->sSize[s][1] = dtemp;
671
672 dtemp = GeoPara->pSize[s][0];
673 GeoPara->pSize[s][0] = GeoPara->pSize[s][1];
674 GeoPara->pSize[s][1] = dtemp;
675 }
676 cout << "swapped x and y" << endl;
677 }
678
679 // get origin
680 fModuleInfo->GetPosition(/*VolumeID, */ 0, 0, 0, padPos, padSize);
681 GeoPara->vOrigin[0] = padPos[0] * 10;
682 GeoPara->vOrigin[1] = padPos[1] * 10;
683 GeoPara->vOrigin[2] = padPos[2] * 10;
684
685 // get col offset
686 fModuleInfo->GetPosition(/*VolumeID, */ 0, 1, 0, padPos, padSize);
687 GeoPara->vX[0] = padPos[0] * 10 - GeoPara->vOrigin[0];
688 GeoPara->vX[1] = padPos[1] * 10 - GeoPara->vOrigin[1];
689 GeoPara->vX[2] = padPos[2] * 10 - GeoPara->vOrigin[2];
690
691 // get row offset
692 // fModuleInfo->GetPosition(VolumeID, 0, 0, 1, padPos, padSize); // rather take next of 3 sectors
693 fModuleInfo->GetPosition(/*VolumeID, */ 1, 0, 0, padPos, padSize);
694 GeoPara->vY[0] = padPos[0] * 10 - GeoPara->vOrigin[0];
695 GeoPara->vY[1] = padPos[1] * 10 - GeoPara->vOrigin[1];
696 GeoPara->vY[2] = padPos[2] * 10 - GeoPara->vOrigin[2];
697
698 // normal vector
699 GeoPara->vN[0] = GeoPara->vX[1] * GeoPara->vY[2] - GeoPara->vX[2] * GeoPara->vY[1];
700 GeoPara->vN[1] = GeoPara->vX[2] * GeoPara->vY[0] - GeoPara->vX[0] * GeoPara->vY[2];
701 GeoPara->vN[2] = GeoPara->vX[0] * GeoPara->vY[1] - GeoPara->vX[1] * GeoPara->vY[0];
702
703 // inclination angle
704 GeoPara->lambda = (GeoPara->vOrigin[0] * GeoPara->vN[0] + GeoPara->vOrigin[1] * GeoPara->vN[1]
705 + GeoPara->vOrigin[2] * GeoPara->vN[2]);
706
707 /*x-direction*/
708 GeoPara->cosX = (GeoPara->vX[0] * 1 + GeoPara->vX[1] * 0 + GeoPara->vX[2] * 0)
709 / sqrt(pow(GeoPara->vX[0], 2) + pow(GeoPara->vX[1], 2) + pow(GeoPara->vX[2], 2));
710 /*y-direction*/
711 GeoPara->cosY = (GeoPara->vY[0] * 0 + GeoPara->vY[1] * 1 + GeoPara->vY[2] * 0)
712 / sqrt(pow(GeoPara->vY[0], 2) + pow(GeoPara->vY[1], 2) + pow(GeoPara->vY[2], 2));
713
714 // check and fix cos values (were larger than 1)
715 // cout << "a cosX " << GeoPara->cosX << " cosY " << GeoPara->cosY << endl;
716
717 if (GeoPara->cosX > 1) GeoPara->cosX = 1;
718 if (GeoPara->cosX < -1) GeoPara->cosX = -1;
719 if (GeoPara->cosY > 1) GeoPara->cosY = 1;
720 if (GeoPara->cosY < -1) GeoPara->cosY = -1;
721
722 // cout << "b cosX " << GeoPara->cosX << " cosY " << GeoPara->cosY << endl;
723
724 if (GeoPara->vX[0] != 0)
725 GeoPara->stepDirection[0] =
726 fabs(GeoPara->vX[0]) / (GeoPara->vX[0]); // is the next pad (1,0,0) on the left or right side?
727 else
728 GeoPara->stepDirection[0] =
729 -fabs(GeoPara->vX[1]) / (GeoPara->vX[1]); // is the next pad (1,0,0) on the left or right side?
730
731 if (GeoPara->vY[1] != 0)
732 GeoPara->stepDirection[1] =
733 fabs(GeoPara->vY[1]) / (GeoPara->vY[1]); // is the next pad (0,1,0) on the upper or lower side?
734 else
735 GeoPara->stepDirection[1] =
736 -fabs(GeoPara->vY[0]) / (GeoPara->vY[0]); // is the next pad (0,1,0) on the upper or lower side?
737
738 // transfer values to old variables
739 Double_t Mpos[3];
740 Double_t Msize[3];
741 for (Int_t i = 0; i < 3; i++) // for x, y and z
742 {
743 Mpos[i] = GeoPara->mPos[i]; // set module position
744 Msize[i] = GeoPara->mSize[i]; // set module size
745 }
746
747 cout.setf(ios::fixed);
748 cout.precision(0);
749 if (Lines)
750 cout << "pos " << setw(6) << Mpos[0] << setw(6) << Mpos[1] << setw(6) << Mpos[2] << " size " << setw(6)
751 << Msize[0] << setw(6) << Msize[1] << setw(6) << Msize[2] << endl;
752
753
754 const Int_t nSec = fModuleInfo->GetNofSectors(); // is always 3
755 // cout << nSec << " nSec" << endl;
756 Double_t Ssize[3 * nSec]; // sector size
757 Double_t Psize[3 * nSec]; // pad size
758 for (Int_t iSec = 0; iSec < nSec; iSec++) // iSec = Sector
759 {
760 Ssize[0 + iSec * nSec] = GeoPara->sSize[iSec][0];
761 Ssize[1 + iSec * nSec] = GeoPara->sSize[iSec][1];
762 Ssize[2 + iSec * nSec] = 0;
763
764 Psize[0 + iSec * nSec] = GeoPara->pSize[iSec][0];
765 Psize[1 + iSec * nSec] = GeoPara->pSize[iSec][1];
766 Psize[2 + iSec * nSec] = 0;
767 }
768
769 if (Lines)
770 cout << "sec " << setw(10) << Ssize[0] << setw(10) << Ssize[1] << setw(10) << Ssize[3] << setw(10) << Ssize[4]
771 << setw(10) << Ssize[6] << setw(10) << Ssize[7] << endl;
772
773 cout.precision(2);
774 if (Lines)
775 cout << "pad " << setw(10) << Psize[0] << setw(10) << Psize[1] << setw(10) << Psize[3] << setw(10) << Psize[4]
776 << setw(10) << Psize[6] << setw(10) << Psize[7] << endl;
777
778 // nCol and nRow
779 Int_t nCol = GeoPara->nCol; // columns per module - or - pads per row
780 Int_t nRow = GeoPara->nRow; // rows per module
781
782 const Int_t mCol = fModuleInfo->GetNofColumns();
783 const Int_t mRow = fModuleInfo->GetNofRows();
784
785 if (Lines) {
786 cout << "col0 " << setw(10) << GeoPara->sCol[0] << " row0 " << setw(10) << GeoPara->sRow[0] << endl;
787 cout << "col1 " << setw(10) << GeoPara->sCol[1] << " row1 " << setw(10) << GeoPara->sRow[1] << endl;
788 cout << "col2 " << setw(10) << GeoPara->sCol[2] << " row2 " << setw(10) << GeoPara->sRow[2] << endl;
789 cout << "col " << setw(10) << nCol << " row " << setw(10) << nRow << endl;
790
791 cout << "mcol " << setw(10) << mCol << " mrow " << setw(10) << mRow << endl;
792 cout << "0col " << setw(10) << fModuleInfo->GetNofColumnsInSector(0) << " 0row " << setw(10)
793 << fModuleInfo->GetNofRowsInSector(0) << endl;
794 cout << "1col " << setw(10) << fModuleInfo->GetNofColumnsInSector(1) << " 1row " << setw(10)
795 << fModuleInfo->GetNofRowsInSector(1) << endl;
796 cout << "2col " << setw(10) << fModuleInfo->GetNofColumnsInSector(2) << " 2row " << setw(10)
797 << fModuleInfo->GetNofRowsInSector(2) << endl;
798 }
799
800
802 // Int_t nCol = 0; // columns per module - or - pads per row
803 // Int_t nRow = 0; // rows per module
804 //
805 // for (Int_t i = 0; i < fnSec; i++)
806 // {
807 // if (Ssize[0+i*nSec] < 2 * Msize[0] && Ssize[0+i*nSec] > 0) // sector size within module size
808 // {
809 // nCol += int(Ssize[0+i*nSec]/Psize[0+i*nSec]);
810 // }
811 // else
812 // {
813 // nCol = int(Ssize[0+i*nSec]/Psize[0+i*nSec]);
814 // }
815 //
816 // if (Ssize[1+i*nSec] < 2 * Msize[1] && Ssize[1+i*nSec] > 0) // sector size within module size
817 // {
818 // nRow += int(Ssize[1+i*nSec]/Psize[1+i*nSec]);
819 // }
820 // else
821 // {
822 // nRow = int(Ssize[1+i*nSec]/Psize[1+i*nSec]);
823 // }
824 // }
825 //
826 // if (Lines)
827 // cout << "nCol " << setw(10) << nCol << " nRow "<< setw(10) << nRow << endl;
828
829 if (Lines) {
830 cout << "vX0 " << setw(10) << GeoPara->vX[0] << endl;
831 cout << "vX1 " << setw(10) << GeoPara->vX[1] << endl;
832 cout << "vY0 " << setw(10) << GeoPara->vY[0] << endl;
833 cout << "vY1 " << setw(10) << GeoPara->vY[1] << endl;
834 cout << "step " << setw(10) << GeoPara->stepDirection[0] << setw(10) << GeoPara->stepDirection[1] << endl;
835 cout << "cos " << setw(10) << GeoPara->cosX << setw(10) << GeoPara->cosY << endl;
836 cout << "v00xy " << setw(10) << GeoPara->mPos[0] - GeoPara->vOrigin[0] << setw(10)
837 << GeoPara->mPos[1] - GeoPara->vOrigin[1] << endl;
838
839 // cout << "------------------------------------------------------" << endl;
840 }
841
842 // angles are relative to y direction
843 if ((GeoPara->stepDirection[0] == 1) && (GeoPara->stepDirection[1] == 1)) // OK - vert up
844 GeoPara->rot_angle = 0;
845 if ((GeoPara->stepDirection[0] == -1) && (GeoPara->stepDirection[1] == 1)) // OK - hori left
846 GeoPara->rot_angle = 1;
847 if ((GeoPara->stepDirection[0] == -1) && (GeoPara->stepDirection[1] == -1)) // OK - vert down
848 GeoPara->rot_angle = 2;
849 if ((GeoPara->stepDirection[0] == 1) && (GeoPara->stepDirection[1] == -1)) // OK - hori right
850 GeoPara->rot_angle = 3;
851
852 // GeoPara->ori = fModuleInfo->GetOrientation(); // is 0,1,2,3
853 if (Lines)
854 cout << "ori " << setw(10) << fModuleInfo->GetOrientation() << " rot " << setw(10) << GeoPara->rot_angle
855 << endl;
856
857 // if (GeoPara->rot_angle == 3)
858 // {
859 // fModuleInfo->GetPosition(VolumeID, 2, GeoPara->sCol[0], GeoPara->sRow[2], padPos, padSize);
860 // if (Lines)
861 // {
868 // cout << "v11xy " << setw(10) << padPos[0] * 10 - GeoPara->vOrigin[0]
869 // << setw(10) << padPos[1] * 10 - GeoPara->vOrigin[1] << endl;
870 // }
871 // GeoPara->vOrigin[0] = padPos[0] * 10;
872 // GeoPara->vOrigin[1] = padPos[1] * 10;
873 // }
874
875 if (Lines) cout << "------------------------------------------------------" << endl;
876
877 Topview[0]->Fill(GeoPara->mPos[0], GeoPara->mPos[2]);
878 Topview[1]->Fill(GeoPara->mPos[0], GeoPara->mPos[1]);
879 Topview[2]->Fill(GeoPara->mPos[2], GeoPara->mPos[1]);
880
881 if (Lines) {
882 DrawPads(GeoPara, Layer, c1);
883 DrawBorders(GeoPara, Layer, c1);
884 }
885 else {
886 Histo(GeoPara, Fast, Layer, c1, HitPad, c2, Topview, c0, mm2bin);
887 }
888
889 if (Lines)
890 for (Int_t s = 0; s < GeoPara->nSec; s++) // for all (3) sectors
891 {
892 fModuleInfo->GetPosition(/*VolumeID, */ s, 0, 0, padPos, padSize);
893
894 TPolyMarker* start = new TPolyMarker(1);
895 start->SetPoint(0, padPos(0) * 10, padPos(1) * 10);
896 start->SetMarkerStyle(22);
897 start->SetMarkerSize(.8);
898
899 if (fDraw) {
900 c1->cd(1);
901 start->Draw("same");
902 }
903
904 fModuleInfo->GetPosition(/*VolumeID, */ s, GeoPara->sCol[0] - 1, GeoPara->sRow[s] - 1, padPos, padSize);
905
906 TPolyMarker* end = new TPolyMarker(1);
907 end->SetPoint(0, padPos(0) * 10, padPos(1) * 10);
908 end->SetMarkerStyle(20);
909 end->SetMarkerSize(.8);
910
911 if (fDraw) {
912 c1->cd(1);
913 end->Draw("same");
914 }
915 }
916 }
917
918 else
919 printf("fModuleInfo == NULL\n");
920}
921
922
923Double_t CbmTrdHitRateQa::CalcHitRate(HitRateGeoPara* GeoPara, Double_t StartX, Double_t /*StopX*/, Int_t xSteps,
924 Double_t StartY, Double_t /*StopY*/, Int_t ySteps, Double_t* /*Mpos*/,
925 TH2F* Topview[3], TCanvas* /*c0*/)
926{
927 //cout << "CalcHitRate" << endl;
928 Double_t HitRate = 0; //1. / sqrt( pow( StartX,2) + pow( StartY,2));
929 Double_t r = 0;
930 Double_t alpha = 0;
931 // Int_t counter = 0;
932 // Double_t a[3] = { 7.66582e+00, 6.97712e+00, 6.53780e+00};
933 // Double_t b[3] = {-2.72375e-03, -1.85168e-03, -1.42673e-03};
934
935 Double_t xStepWidth = GeoPara->cosX;
936 Double_t yStepWidth = GeoPara->cosY;
937 Double_t x = StartX + 0.5 * GeoPara->cosX; // don't start on the pad border line
938 Double_t y = StartY + 0.5 * GeoPara->cosY;
939
940 // Double_t xStepWidth = fabs(GeoPara->cosX);
941 // Double_t yStepWidth = fabs(GeoPara->cosY);
942 // Double_t x = StartX + 0.5 * fabs(GeoPara->cosX); // don't start on the pad border line
943 // Double_t y = StartY + 0.5 * fabs(GeoPara->cosY);
944 Double_t z = (GeoPara->lambda - (x * GeoPara->vN[0] + y * GeoPara->vN[1])) / GeoPara->vN[2];
945 for (Int_t yStep = 0; yStep < ySteps; yStep++) {
946 x = StartX + 0.5 * GeoPara->cosX;
947 // x = StartX + 0.5 * fabs(GeoPara->cosX);
948 for (Int_t xStep = 0; xStep < xSteps; xStep++) {
949 z = (GeoPara->lambda - (x * GeoPara->vN[0] + y * GeoPara->vN[1])) / GeoPara->vN[2];
950 r = sqrt(pow(x, 2) + pow(y, 2));
951 alpha = atan(r / z) * 1000.;
952 /* //Fit without errors
953 HitRate +=
954 exp(4.54156e00 + -8.47377e-03 * alpha) +
955 exp(2.40005e01 + -1.19541e-02 * alpha) /
956 (z * z)
957 ;
958 */
959 HitRate += // fit including errors
960 exp(4.536355e00 + -8.393716e-03 * alpha) + exp(2.400547e01 + -1.208306e-02 * alpha) / (z * z);
961
962 Topview[0]->Fill(x, z);
963 Topview[1]->Fill(x, y);
964 Topview[2]->Fill(z, y);
965
966 x += xStepWidth;
967 }
968 y += yStepWidth;
969 }
970 return (HitRate /*/(counter/100.)*/); /*Convertes Hits/Pad -> Hits/cm² on each Pad*/
971}
972
973
974void CbmTrdHitRateQa::Histo(HitRateGeoPara* GeoPara, Bool_t Fast, TH2F* h2Layer, TCanvas* /*c1*/, TH1F* h1HitPad,
975 TCanvas* /*c2*/, TH2F* Topview[3], TCanvas* c0, Double_t mm2bin)
976{
977 Double_t ZRangeL = 1e00; //1e05;
978 Double_t ZRangeU = 1e05; //1e06;
979 TString name;
980
981 // show only current module name
982 name.Form("ModuleID %5d", GeoPara->moduleId);
983 // cout << " " << name << "\r" << flush;
984 cout << " " << name << "\n" << flush;
985
986 name.Form("Module%05d", GeoPara->moduleId);
987 TH2F* h2Module;
988 if (GeoPara->mSize[0] < 500) // small module
989 h2Module = new TH2F(name, name, 600 / mm2bin + 1, -300.5, 300.5, 600 / mm2bin + 1, -300.5,
990 300.5); // 501 x 501 bins
991 else // large module
992 h2Module = new TH2F(name, name, 1000 / mm2bin + 1, -500.5, 500.5, 1000 / mm2bin + 1, -500.5,
993 500.5); // 1001 x 1001 bins
994 // TH2F *h2Module = new TH2F(name,name,1500/mm2bin+1,-750.5,750.5,1500/mm2bin+1,-750.5,750.5);
995 h2Module->GetZaxis()->SetRangeUser(ZRangeL, ZRangeU);
996
997 name.Form("Module%05dHitPad", GeoPara->moduleId);
998 TH1F* h1HitPadModule = new TH1F(name, name, 1200, 0, 120000);
999 // TH1F* h1HitPadModule = new TH1F(name,name,10000,1e00,1e06);
1000
1001 //cout << "Histo" << endl;
1002
1003 // 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]);
1004 // cout << endl << "Xp " << GeoPara->mPos[0] << " Xs " << GeoPara->mSize[0] << " Yp " << GeoPara->mPos[1] << " Ys " << GeoPara->mSize[1] << endl;
1005
1006 Double_t HiteRate = 0;
1007
1008 Int_t iSecX = 0;
1009 Int_t iSecY = 0;
1010
1011 Double_t planeStartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1012 Double_t planeStopX = planeStartX + GeoPara->pSize[iSecX][0];
1013
1014 Double_t planeStartY = GeoPara->mPos[1] - GeoPara->mSize[1];
1015 Double_t planeStopY = planeStartY + GeoPara->pSize[iSecY][1];
1016
1017 Double_t StartX = GeoPara->vOrigin[0]
1018 + 0.5 * GeoPara->pSize[iSecX][0] * GeoPara->cosX; // vOrigin points to the center of pad (0,0,0)
1019 Double_t StopX = StartX + GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1020
1021 Double_t StartY = GeoPara->vOrigin[1] + 0.5 * GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1022 Double_t StopY = StartY + GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1023
1024 Int_t xSteps = 0;
1025 Int_t ySteps = 0;
1026
1027 // cout << "x " << GeoPara->pSize[iSecX][0] << " y " << GeoPara->pSize[iSecX][1] << endl;
1028
1029 // printf("stX %6d spX %6d stY %6d spY %6d\n", (Int_t)StartX, (Int_t)StopX, (Int_t)StartY, (Int_t)StopY);
1030
1031 //Int_t xStepDirection = GeoPara->vX[0] / fabs(GeoPara->vX[0]); // is the next pad (1,0,0) on the left or right side?
1032 //Int_t yStepDirection = GeoPara->vY[1] / fabs(GeoPara->vY[1]); // is the next pad (0,1,0) on the upper or lower side?
1033
1034 //----------------------------------------------------------------------------------------
1035
1036 //if (GeoPara->rot_angle == 1)
1037 //if (GeoPara->rot_angle != 3)
1038 for (Int_t iR = 0; iR < GeoPara->nRow; iR++) // rows
1039 {
1040 StartX = GeoPara->vOrigin[0] + 0.5 * GeoPara->pSize[iSecX][0] * GeoPara->cosX; // why to the left?
1041 StopX = StartX + GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1042 planeStartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1043 planeStopX = planeStartX + GeoPara->pSize[iSecX][0];
1044
1045 for (Int_t iC = 0; iC < GeoPara->nCol; iC++) // for all columns
1046 {
1047 xSteps = GeoPara->pSize[iSecX][0] - 1; // 1 Step / mm
1048 ySteps = GeoPara->pSize[iSecY][1] - 1;
1049
1050 if (Fast) {
1051 //DE
1052 // if (GeoPara->mPos[0] > -1 && GeoPara->mPos[1] > -1)
1053 if (GeoPara->mPos[0] < 1 && GeoPara->mPos[1] < 1) // copy the bottom left quadrant
1054 {
1055 HiteRate = CalcHitRate(GeoPara, StartX, StopX, xSteps, StartY, StopY, ySteps, GeoPara->mPos, Topview, c0);
1056 h1HitPad->Fill(HiteRate, 4); // weight 4 for each quadrant
1057 }
1058 }
1059
1060 else // slow, for each module
1061 {
1062 HiteRate = CalcHitRate(GeoPara, StartX, StopX, xSteps, StartY, StopY, ySteps, GeoPara->mPos, Topview, c0);
1063 h1HitPad->Fill(HiteRate);
1064 h1HitPadModule->Fill(HiteRate); // single module, goes into root file
1065 }
1066
1067 // fill HitRate into h2Layer histo
1068
1069 Int_t mStepX = Int_t((planeStartX - GeoPara->mPos[0]));
1070 Int_t mStepY = Int_t((planeStartY - GeoPara->mPos[1]));
1071
1072 for (Int_t stepY = int(planeStartY / mm2bin); stepY < int(planeStopY / mm2bin); stepY++) {
1073 // mStepX = Int_t((planeStartX-GeoPara->mPos[0]));
1074
1075 for (Int_t stepX = int(planeStartX / mm2bin); stepX < int(planeStopX / mm2bin); stepX++) {
1076
1077 if (Fast) {
1078 //DE
1079 // if (GeoPara->mPos[0]/*GeoPara->mPos[0]*/ > -1 && GeoPara->mPos[1]/*GeoPara->mPos[1]*/ > -1)
1080 if (GeoPara->mPos[0] /*GeoPara->mPos[0]*/ < 1 && GeoPara->mPos[1] /*GeoPara->mPos[1]*/ < 1) {
1081
1082 h2Layer->SetBinContent(stepX + int(winsize / mm2bin), stepY + int(winsize / mm2bin),
1083 HiteRate); // bin coordiantes
1084 h2Layer->SetBinContent(-1 * stepX + int(winsize / mm2bin), stepY + int(winsize / mm2bin),
1085 HiteRate); // bin coordiantes
1086 h2Layer->SetBinContent(stepX + int(winsize / mm2bin), -1 * stepY + int(winsize / mm2bin),
1087 HiteRate); // bin coordiantes
1088 h2Layer->SetBinContent(-1 * stepX + int(winsize / mm2bin), -1 * stepY + int(winsize / mm2bin),
1089 HiteRate); // bin coordiantes
1090
1091 if (GeoPara->moduleId == 5) {
1092 // cout << "filling 5 " << mStepX << " " << mStepY << " " << HiteRate << endl;
1093 h2Module->Fill(mStepX, mStepY,
1094 HiteRate); // single module, goes into root file
1095 }
1096 }
1097 }
1098
1099 else // slow
1100 {
1101
1102 h2Layer->SetBinContent(stepX + int(winsize / mm2bin), stepY + int(winsize / mm2bin),
1103 HiteRate); // bin coordinates
1104
1105 // cout << "StepX " << stepX << " StepY " << stepY << endl;
1106
1107 if (GeoPara->moduleId == 5) {
1108 cout << "filling 5 " << mStepX << " " << mStepY << " " << HiteRate << endl;
1109 h2Module->Fill(mStepX, mStepY,
1110 HiteRate); // single module, goes into root file
1111 }
1112 // cout << "mStepX " << mStepX << " mStepY " << mStepY << endl;
1113 }
1114
1115 mStepX += mm2bin;
1116 }
1117
1118 mStepY += mm2bin;
1119 }
1120
1121 //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);
1122
1123 // if (iC == GeoPara->sCol[iSecX]-1) // why -1 ??
1124 if (iC == GeoPara->sCol[iSecX]) {
1125 // cout << "iC " << iC << " " << GeoPara->nCol << " " << iSecX << endl;
1126 iSecX++;
1127 }
1128 StartX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1129 StopX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1130 planeStartX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0];
1131 planeStopX += GeoPara->stepDirection[0] * GeoPara->pSize[iSecX][0];
1132 // StartX += GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1133 // StopX += GeoPara->pSize[iSecX][0] * GeoPara->cosX;
1134 // planeStartX += GeoPara->pSize[iSecX][0];
1135 // planeStopX += GeoPara->pSize[iSecX][0];
1136 }
1137
1138 iSecX = 0;
1139
1140 // if (iR == GeoPara->sRow[iSecY]-1) // why -1 ??
1141 if (iR == GeoPara->sRow[iSecY]) {
1142 // cout << "iR " << iR << " " << GeoPara->nRow << " " << iSecY << endl;
1143 iSecY++;
1144 }
1145
1146 StartY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1147 StopY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1148 planeStartY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1];
1149 planeStopY += GeoPara->stepDirection[1] * GeoPara->pSize[iSecY][1];
1150 // StartY += GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1151 // StopY += GeoPara->pSize[iSecY][1] * GeoPara->cosY;
1152 // planeStartY += GeoPara->pSize[iSecY][1];
1153 // planeStopY += GeoPara->pSize[iSecY][1];
1154 }
1155
1156 //----------------------------------------------------------------------------------------
1157
1158 for (Int_t i = 0; i < 3; i++)
1159 Topview[i]->Write("", TObject::kOverwrite);
1160
1161 h2Layer->Write("", TObject::kOverwrite);
1162 h1HitPad->Write("", TObject::kOverwrite);
1163 h2Module->Write("", TObject::kOverwrite);
1164 h1HitPadModule->Write("", TObject::kOverwrite);
1165
1166 delete h2Module;
1167 delete h1HitPadModule;
1168}
1169
1170
1171void CbmTrdHitRateQa::DrawBorders(HitRateGeoPara* GeoPara, TH2F* /*Layer*/, TCanvas* c1)
1172{
1173 //----------------------Sector--------------------------------------
1174 Double_t SecXStart = 0.0;
1175 Double_t SecYStart = 0.0;
1176 Double_t SecXStop = 0.0;
1177 Double_t SecYStop = 0.0;
1178
1179 for (Int_t iSec = 0; iSec < GeoPara->nSec - 1; iSec++) // would be enough to iterate up to nSec-1
1180 {
1181
1182 if (GeoPara->sSize[iSec][0] < 2 * GeoPara->mSize[0] && GeoPara->sSize[iSec][0] > 0) // in x direction
1183 {
1184 SecXStart += GeoPara->sSize[iSec][0];
1185 SecXStop = SecXStart;
1186 }
1187 else {
1188 SecXStart = 0.0;
1189 SecXStop = GeoPara->sSize[iSec][0];
1190 }
1191
1192 if (GeoPara->sSize[iSec][1] < 2 * GeoPara->mSize[1] && GeoPara->sSize[iSec][1] > 0) // in x direction
1193 {
1194 SecYStart += GeoPara->sSize[iSec][1];
1195 SecYStop = SecYStart;
1196 }
1197 else {
1198 SecYStart = 0.0;
1199 SecYStop = GeoPara->sSize[iSec][1];
1200 }
1201
1202
1203 TLine* sector =
1204 new TLine(GeoPara->mPos[0] - GeoPara->mSize[0] + SecXStart, GeoPara->mPos[1] - GeoPara->mSize[1] + SecYStart,
1205 GeoPara->mPos[0] - GeoPara->mSize[0] + SecXStop, GeoPara->mPos[1] - GeoPara->mSize[1] + SecYStop);
1206 sector->SetLineColor(15);
1207 sector->SetLineStyle(2);
1208
1209 TPolyMarker* corner = new TPolyMarker(1);
1210 corner->SetPoint(0, GeoPara->mPos[0] - GeoPara->stepDirection[0] * GeoPara->mSize[0],
1211 GeoPara->mPos[1] - GeoPara->stepDirection[1] * GeoPara->mSize[1]);
1212 corner->SetMarkerStyle(21);
1213 corner->SetMarkerSize(.8);
1214
1215 if (fDraw) {
1216 c1->cd(1);
1217 sector->Draw("same");
1218 corner->Draw("same");
1219 }
1220 }
1221
1222
1223 //----------------------Module--------------------------------------
1224 TBox* module = new TBox(GeoPara->mPos[0] - GeoPara->mSize[0], GeoPara->mPos[1] - GeoPara->mSize[1],
1225 GeoPara->mPos[0] + GeoPara->mSize[0], GeoPara->mPos[1] + GeoPara->mSize[1]);
1226 module->SetFillStyle(0);
1227 module->SetLineColor(1);
1228
1229 TBox* M_inner =
1230 new TBox(GeoPara->mPos[0] - 100, GeoPara->mPos[1] - 100, GeoPara->mPos[0] + 100, GeoPara->mPos[1] + 100);
1231 // M_inner->SetUniqueID(ModuleId);
1232 M_inner->SetFillColor(kWhite);
1233
1234 if (fDraw) {
1235 c1->cd(1);
1236 module->Draw("same"); // black module frame
1237 // M_inner->Draw("same"); // white box in module center // do not draw fot the time being
1238 c1->Write("", TObject::kOverwrite);
1239 }
1240}
1241// --------------------------------------------------------------------
1242
1243
1244void CbmTrdHitRateQa::DrawPads(HitRateGeoPara* GeoPara, TH2F* /*Layer*/, TCanvas* c1)
1245{
1246 //----------------------Pad--------------------------------------
1247 // const Int_t nR = GeoPara->nRow;
1248 // const Int_t nC = GeoPara->nCol;
1249
1250 Int_t iSecX = 0;
1251 Int_t iSecY = 0;
1252
1253 Double_t StartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1254 Double_t StartY = GeoPara->mPos[1] - GeoPara->mSize[1];
1255
1256 Int_t SecCol = GeoPara->sCol[iSecX];
1257 Int_t SecRow = GeoPara->sRow[iSecY];
1258
1259 for (Int_t iR = 0; iR < GeoPara->nRow; iR++) // rows
1260 {
1261 StartX = GeoPara->mPos[0] - GeoPara->mSize[0];
1262 for (Int_t iC = 0; iC < GeoPara->nCol; iC++) // cols
1263 {
1264 TLine* Pa = new TLine(StartX, // - low
1265 StartY, StartX + GeoPara->pSize[iSecX][0], StartY);
1266 TLine* Pb = new TLine(StartX, // - high
1267 StartY + GeoPara->pSize[iSecY][1], StartX + GeoPara->pSize[iSecX][0],
1268 StartY + GeoPara->pSize[iSecY][1]);
1269 TLine* Pc = new TLine(StartX, // | left
1270 StartY, StartX, StartY + GeoPara->pSize[iSecY][1]);
1271 TLine* Pd = new TLine(StartX + GeoPara->pSize[iSecX][0], // | right
1272 StartY, StartX + GeoPara->pSize[iSecX][0], StartY + GeoPara->pSize[iSecY][1]);
1273 c1->cd(1);
1274 Pa->SetLineColor(17);
1275 //Pa->SetLineStyle(3);
1276 Pa->Draw("same");
1277
1278 Pb->SetLineColor(17);
1279 //Pb->SetLineStyle(3);
1280 Pb->Draw("same");
1281
1282 Pc->SetLineColor(17);
1283 //Pc->SetLineStyle(3);
1284 Pc->Draw("same");
1285
1286 Pd->SetLineColor(17);
1287 //Pd->SetLineStyle(3);
1288 Pd->Draw("same");
1289
1290 //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);
1291
1292 if (iC == SecCol - 1) {
1293 iSecX++;
1294 SecCol += GeoPara->sCol[iSecX];
1295 }
1296 StartX += GeoPara->pSize[iSecX][0];
1297 }
1298 iSecX = 0;
1299 if (iR == SecRow - 1) {
1300 iSecY++;
1301 SecRow += GeoPara->sRow[iSecY];
1302 }
1303 StartY += GeoPara->pSize[iSecY][1];
1304 }
1305}
1306// --------------------------------------------------------------------
1307
1308
1309// ------DrawDigi------------------------------------------------------
1310void CbmTrdHitRateQa::DrawDigi() {}
1311// --------------------------------------------------------------------
1312
1313
1314// ---- Register ------------------------------------------------------
1315void CbmTrdHitRateQa::Register()
1316{
1317 //FairRootManager::Instance()->Register("TrdDigi","Trd Digi", fDigiCollection, IsOutputBranchPersistent("TrdDigi"));
1318 //FairRootManager::Instance()->Register("TrdDigiMatch","Trd Digi Match", fDigiMatchCollection, IsOutputBranchPersistent("TrdDigiMatch"));
1319}
1320// --------------------------------------------------------------------
1321
1322ClassImp(CbmTrdHitRateQa)
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
HitRateGeoPara
struct HitRateGeoPara HitRateGeoPara
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
exp
friend fvec exp(const fvec &a)
Definition KfSimdPseudo.h:121
CbmTrdAddress::GetLayerId
static uint32_t GetLayerId(uint32_t address)
Return layer ID from address.
Definition CbmTrdAddress.h:69
CbmTrdGeoHandler::Init
void Init(Bool_t isSimulation=kFALSE)
Definition CbmTrdGeoHandler.cxx:52
CbmTrdHitRateQa::Exec
virtual void Exec(Option_t *option)
Definition CbmTrdHitRateQa.cxx:211
CbmTrdHitRateQa::GetModuleInformationFromDigiPar
void GetModuleInformationFromDigiPar(HitRateGeoPara *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 CbmTrdHitRateQa.cxx:570
CbmTrdHitRateQa::DrawBorders
void DrawBorders(HitRateGeoPara *GeoPara, TH2F *Layer, TCanvas *c1)
Definition CbmTrdHitRateQa.cxx:1171
CbmTrdHitRateQa::fDigiCollection
TClonesArray * fDigiCollection
Trd MC points.
Definition CbmTrdHitRateQa.h:187
CbmTrdHitRateQa::ReInit
virtual InitStatus ReInit()
Definition CbmTrdHitRateQa.cxx:155
CbmTrdHitRateQa::fAsicPar
CbmTrdParSetAsic * fAsicPar
MC Track information.
Definition CbmTrdHitRateQa.h:191
CbmTrdHitRateQa::fMCStacks
TClonesArray * fMCStacks
Corresponding MCPoints to TRD digis.
Definition CbmTrdHitRateQa.h:189
CbmTrdHitRateQa::HistoInit
void HistoInit(TCanvas *&c1, TCanvas *&c2, TH2F *&Layer, TH1F *&HitPad, Double_t ZRangeL, Double_t ZRangeU, Double_t mm2bin)
Definition CbmTrdHitRateQa.cxx:460
CbmTrdHitRateQa::Histo
void Histo(HitRateGeoPara *GeoPara, Bool_t Fast, TH2F *Layer, TCanvas *c1, TH1F *HitPad, TCanvas *c2, TH2F *Topview[3], TCanvas *c0, Double_t mm2bin)
Definition CbmTrdHitRateQa.cxx:974
CbmTrdHitRateQa::SetParContainers
virtual void SetParContainers()
Definition CbmTrdHitRateQa.cxx:139
CbmTrdHitRateQa::DrawPads
void DrawPads(HitRateGeoPara *GeoPara, TH2F *Layer, TCanvas *c1)
Definition CbmTrdHitRateQa.cxx:1244
CbmTrdHitRateQa::FinishEvent
virtual void FinishEvent()
Definition CbmTrdHitRateQa.cxx:520
CbmTrdHitRateQa::fDigiMap
std::map< std::pair< Int_t, std::pair< Int_t, Int_t > >, CbmTrdDigi * > fDigiMap
Definition CbmTrdHitRateQa.h:199
CbmTrdHitRateQa::CalcHitRate
Double_t CalcHitRate(HitRateGeoPara *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 CbmTrdHitRateQa.cxx:923
CbmTrdHitRateQa::fTrdPoints
TClonesArray * fTrdPoints
Definition CbmTrdHitRateQa.h:186
CbmTrdHitRateQa::fGeoPar
CbmTrdParSetGeo * fGeoPar
Definition CbmTrdHitRateQa.h:193
CbmTrdHitRateQa::fGeoHandler
CbmTrdGeoHandler * fGeoHandler
Definition CbmTrdHitRateQa.h:195
CbmTrdHitRateQa::Init
virtual InitStatus Init()
Definition CbmTrdHitRateQa.cxx:173
CbmTrdHitRateQa::fDigiMatchCollection
TClonesArray * fDigiMatchCollection
TRD digis.
Definition CbmTrdHitRateQa.h:188
CbmTrdHitRateQa::fDigiPar
CbmTrdParSetDigi * fDigiPar
Definition CbmTrdHitRateQa.h:192
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
HitRateGeoPara::pSize
Double_t pSize[3][3]
Definition CbmTrdHitRateQa.h:36
HitRateGeoPara::mSize
Double_t mSize[3]
Definition CbmTrdHitRateQa.h:33
HitRateGeoPara::vOrigin
Double_t vOrigin[3]
Definition CbmTrdHitRateQa.h:46
HitRateGeoPara::sSize
Double_t sSize[3][3]
Definition CbmTrdHitRateQa.h:35
HitRateGeoPara::stepDirection
Int_t stepDirection[3]
Definition CbmTrdHitRateQa.h:54
HitRateGeoPara::mPos
Double_t mPos[3]
Definition CbmTrdHitRateQa.h:32