CbmRoot
CbmTrdRecoQa.cxx
Go to the documentation of this file.
1 // -----------------------------------------------------------------------
2 // ----- CbmTrdRecoQa -----
3 // ----- Created 24.02.07 by F. Uhlig -----
4 // -----------------------------------------------------------------------
5 
6 #include "CbmTrdRecoQa.h"
7 
8 #include "CbmMCTrack.h"
9 #include "CbmTrdCluster.h"
10 #include "CbmTrdDigi.h"
11 #include "CbmTrdHit.h"
12 #include "CbmTrdPoint.h"
13 
14 #include "CbmMCTrack.h"
15 #include "CbmTrdParModDigi.h"
16 #include "CbmTrdParSetDigi.h"
17 
18 #include "FairLogger.h"
19 #include "FairRootManager.h"
20 #include "FairRunAna.h"
21 #include "FairRuntimeDb.h"
22 
23 #include "CbmDigiManager.h"
24 #include "CbmTrdGeoHandler.h"
25 #include "CbmTrdUtils.h"
26 #include "TCanvas.h"
27 #include "TClonesArray.h"
28 #include "TColor.h"
29 #include "TFrame.h"
30 #include "TGeoManager.h"
31 #include "TGraph.h"
32 #include "TGraphErrors.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TLine.h"
36 #include "TMath.h"
37 #include "TMultiGraph.h"
38 #include "TPad.h"
39 #include "TPaveText.h"
40 #include "TPolyLine.h"
41 #include "TStopwatch.h"
42 #include "TVector3.h"
43 
44 
45 #include <iostream>
46 
47 using std::cout;
48 using std::endl;
49 using std::map;
50 
51 // ---- Default constructor -------------------------------------------------
52 
53 CbmTrdRecoQa::CbmTrdRecoQa() : CbmTrdRecoQa("TrdRecoQa", "") {}
54 
55 // ---- Standard constructor ------------------------------------------------
56 CbmTrdRecoQa::CbmTrdRecoQa(const char* name, const char*)
57  : FairTask(name)
58  , fTrianglePads(false)
59  , fTriggerTH(1.0e-6)
60  , fClusters(NULL)
61  , fHits(NULL)
62  , fMCPoints(NULL)
63  , fDigiPar(NULL)
64  , fModuleInfo(NULL)
65  , fGeoHandler(new CbmTrdGeoHandler())
66  , fModuleMap()
67  , fModuleMapPoint()
68  , fModuleMapDigi()
69  , fModuleMapCluster()
70  , fModuleMapHit()
71  , fModuleMapTrack() {}
72 // --------------------------------------------------------------------------
73 
74 
75 // ---- Destructor ----------------------------------------------------------
76 CbmTrdRecoQa::~CbmTrdRecoQa() {}
77 // --------------------------------------------------------------------------
78 void CbmTrdRecoQa::SetTriggerThreshold(Double_t minCharge) {
79  fTriggerTH = minCharge; // To be used for test beam data processing
80 }
81 
82 // ---- Initialisation ------------------------------------------------------
83 InitStatus CbmTrdRecoQa::Init() {
84  // Get pointer to the ROOT I/O manager
85  FairRootManager* rootMgr = FairRootManager::Instance();
86  if (NULL == rootMgr) {
87  cout << "-E- CbmTrdRecoQa::Init : "
88  << "ROOT manager is not instantiated !" << endl;
89  return kFATAL;
90  }
91 
92 
93  // Get pointer to TRD point array
94  fMCPoints = (TClonesArray*) rootMgr->GetObject("TrdPoint");
95  if (NULL == fMCPoints) {
96  cout << "-W- CbmTrdRecoQa::Init : "
97  << "no MC point array !" << endl;
98  return kERROR;
99  }
100 
101  CbmDigiManager::Instance()->Init();
102  if (!CbmDigiManager::Instance()->IsPresent(ECbmModuleId::kTrd)) LOG(fatal);
103 
104 
105  // Get pointer to TRD point array
106  fClusters = (TClonesArray*) rootMgr->GetObject("TrdCluster");
107  if (NULL == fClusters) {
108  cout << "-W- CbmTrdRecoQa::Init : "
109  << "no cluster array !" << endl;
110  return kERROR;
111  }
112 
113  // Get pointer to TRD point array
114  fHits = (TClonesArray*) rootMgr->GetObject("TrdHit");
115  if (NULL == fHits) {
116  cout << "-W- CbmTrdRecoQa::Init : "
117  << "no hit array !" << endl;
118  return kERROR;
119  }
120 
121  fGeoHandler->Init();
122 
123  return kSUCCESS;
124 }
125 
126 // --------------------------------------------------------------------------
127 
128 // ---- Initialisation ----------------------------------------------
129 void CbmTrdRecoQa::SetParContainers() {
130  cout << " * CbmTrdRecoQa * :: SetParContainers() " << endl;
131 
132 
133  // Get Base Container
134  FairRunAna* ana = FairRunAna::Instance();
135  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
136 
137  fDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
138 }
139 // --------------------------------------------------------------------
140 
141 // ---- ReInit -------------------------------------------------------
142 InitStatus CbmTrdRecoQa::ReInit() {
143 
144  cout << " * CbmTrdRecoQa * :: ReInit() " << endl;
145 
146 
147  FairRunAna* ana = FairRunAna::Instance();
148  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
149 
150  fDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
151 
152  return kSUCCESS;
153 }
154 void CbmTrdRecoQa::SetTriangularPads(Bool_t triangles) {
155  fTrianglePads = triangles;
156 }
157 /*
158 TPolyLine *CbmTrdRecoQa::CreateTriangularPad(Int_t column, Int_t row, Double_t content){
159  const Int_t nCoordinates = 4;
160  Double_t x[nCoordinates] = {column-0.5,column+0.5,column+0.5,column-0.5};
161  Double_t y[nCoordinates] = {row-0.5, row-0.5, row+0.5, row-0.5 };
162  if (row%2 != 0){
163  y[1] = row+0.5;
164  y[2] = row+0.5;
165  }
166  TPolyLine *pad = new TPolyLine(nCoordinates,x,y);
167  pad->SetFillColor(content);
168  return pad;
169 }
170 */
171 // ---- Task execution ------------------------------------------------------
172 void CbmTrdRecoQa::Exec(Option_t*) {
173  CbmTrdUtils* utils = new CbmTrdUtils();
174  TStopwatch timer;
175  timer.Start();
176  cout << "================CbmTrdRecoQa==============" << endl;
177  // Declare variables
178  CbmTrdPoint* point = NULL;
179  const CbmTrdDigi* digi = NULL;
180  CbmTrdCluster* cluster = NULL;
181  CbmTrdHit* hit = NULL;
182 
183  Int_t nPoints(0), nDigis(0), nClusters(0), nHits(0);
184  nDigis = CbmDigiManager::Instance()->GetNofDigis(ECbmModuleId::kTrd);
185  if (fClusters) nClusters = fClusters->GetEntries();
186  if (fMCPoints) nPoints = fMCPoints->GetEntries();
187  if (fHits) nHits = fHits->GetEntries();
188 
189  TH1D* digiMaxSpectrum =
190  new TH1D("digiMaxSpectrum", "digiMaxSpectrum", 10000, 0, 1e-4);
191  TH1D* digiSpectrum = new TH1D("digiSpectrum", "digiSpectrum", 10000, 0, 1e-4);
192 
193  LOG(info) << "CbmTrdRecoQa::Exec : MC-Points:" << nPoints;
194  LOG(info) << "CbmTrdRecoQa::Exec : Digis: " << nDigis;
195  LOG(info) << "CbmTrdRecoQa::Exec : Cluster: " << nClusters;
196  LOG(info) << "CbmTrdRecoQa::Exec : Hits: " << nHits;
197  Int_t moduleAddress(-1), moduleId(-1);
198 
199  // Double_t dummy_x(-1), dummy_y(-1);
200  TString name;
201  //cout << "Points" << endl;
202  Int_t pointCounter = 0;
203  for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
204  //cout << iPoint << endl;
205  point = (CbmTrdPoint*) fMCPoints->At(iPoint);
206  Double_t in[3] = {point->GetXIn(), point->GetYIn(), point->GetZIn()};
207  Double_t out[3] = {point->GetXOut(), point->GetYOut(), point->GetZOut()};
208  gGeoManager->FindNode(
209  (out[0] + in[0]) / 2, (out[1] + in[1]) / 2, (out[2] + in[2]) / 2);
210  if (!TString(gGeoManager->GetPath()).Contains("gas")) {
211  LOG(error) << "CbmTrdRecoQa::Exec: MC-track not in TRD! Node:"
212  << TString(gGeoManager->GetPath()).Data()
213  << " gGeoManager->MasterToLocal() failed!";
214  continue;
215  }
216  //std::map<Int_t, std::vector<TH2D*> >::iterator it = fModuleMap.find(moduleAddress);
217  moduleAddress = CbmTrdAddress::GetModuleAddress(point->GetDetectorID()); //
218  moduleId = CbmTrdAddress::GetModuleId(point->GetDetectorID());
219  //printf("Address:%i ID:%i\n",moduleAddress,moduleId);
220  fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(
221  moduleAddress);
222  //printf("Address:%i ID:%i\n",moduleAddress,moduleId);
223  if (fModuleInfo) {
224  std::map<Int_t, TGraphErrors*>::iterator it =
225  fModuleMapPoint.find(moduleAddress);
226  if (it == fModuleMapPoint.end()) {
227  name.Form("ModuleAddress%05i", moduleAddress);
228  fModuleMap[moduleAddress] = new TCanvas(name, name, 1000, 1000);
229  fModuleMap[moduleAddress]->Divide(2, 2);
230  //name.Form("ModuleAddress%iPoints",moduleAddress);
231  fModuleMapPoint[moduleAddress] = new TGraphErrors();
232  fModuleMapPoint[moduleAddress]->SetMarkerStyle(20);
233  fModuleMapPoint[moduleAddress]->SetMarkerSize(0.5);
234  fModuleMapPoint[moduleAddress]->SetMarkerColor(15);
235  //fModuleMapPoint[moduleAddress] = new TH2I(name,name,fModuleInfo->GetNofColumns(),-0.5,fModuleInfo->GetNofColumns()-0.5,fModuleInfo->GetNofRows(),-0.5,fModuleInfo->GetNofRows()-0.5);
236  TH2I* dummy = new TH2I(name,
237  name,
238  fModuleInfo->GetNofColumns(),
239  -0.5,
240  fModuleInfo->GetNofColumns() - 0.5,
241  fModuleInfo->GetNofRows(),
242  -0.5,
243  fModuleInfo->GetNofRows() - 0.5);
244  dummy->SetStats(kFALSE);
245  name.Form("ModuleAddress%05iDigis", moduleAddress);
246  fModuleMapDigi[moduleAddress] =
247  new TH2D(name,
248  name,
249  fModuleInfo->GetNofColumns(),
250  -0.5,
251  fModuleInfo->GetNofColumns() - 0.5,
252  fModuleInfo->GetNofRows(),
253  -0.5,
254  fModuleInfo->GetNofRows() - 0.5);
255  fModuleMapDigi[moduleAddress]->SetContour(99);
256  fModuleMapDigi[moduleAddress]->SetStats(kFALSE);
257  name.Form("ModuleAddress%05iClusters", moduleAddress);
258  fModuleMapCluster[moduleAddress] =
259  new TH2I(name,
260  name,
261  fModuleInfo->GetNofColumns(),
262  -0.5,
263  fModuleInfo->GetNofColumns() - 0.5,
264  fModuleInfo->GetNofRows(),
265  -0.5,
266  fModuleInfo->GetNofRows() - 0.5);
267  fModuleMapCluster[moduleAddress]->SetContour(99);
268  fModuleMapCluster[moduleAddress]->SetStats(kFALSE);
269  name.Form("ModuleAddress%05iHits", moduleAddress);
270  fModuleMapHit[moduleAddress] =
271  new TGraphErrors(); //name,name,fModuleInfo->GetNofColumns(),-0.5,fModuleInfo->GetNofColumns()-0.5,fModuleInfo->GetNofRows(),-0.5,fModuleInfo->GetNofRows()-0.5);
272  fModuleMapHit[moduleAddress]->SetMarkerStyle(24);
273  fModuleMapTrack[moduleAddress] = new std::vector<TLine*>;
274  fModuleMap[moduleAddress]->cd(1);
275  dummy->Draw(); /*
276  fModuleMapPoint[moduleAddress]->Draw("AP");
277  fModuleMapPoint[it->first]->SetMaximum(fModuleMapDigi[it->first]->GetYaxis()->GetXmax());
278  fModuleMapPoint[it->first]->SetMinimum(fModuleMapDigi[it->first]->GetYaxis()->GetXmin());
279  fModuleMapPoint[it->first]->GetXaxis()->SetLimits(fModuleMapDigi[it->first]->GetXaxis()->GetXmin(),
280  fModuleMapDigi[it->first]->GetXaxis()->GetXmax());
281  fModuleMapPoint[moduleAddress]->Draw("AP");
282  */
283  }
284  //cout << iPoint << endl;
285 
286 
287  Double_t local_in[3];
288  Double_t local_out[3];
289  gGeoManager->MasterToLocal(in, local_in);
290  gGeoManager->MasterToLocal(out, local_out);
291 
292  for (Int_t i = 0; i < 3; i++)
293  local_out[i] =
294  local_in[i]
295  + 0.975
296  * (local_out[i]
297  - local_in
298  [i]); // cut last 2.5% of tracklet, to move exit point within gas volume
299 
300  Int_t row_in(0), row_out(0), col_in(0), col_out(0), sec_in(0), sec_out(0);
301  Double_t x_in(0), y_in(0), x_out(0), y_out(0);
302  if (!fModuleInfo->GetPadInfo(local_in, sec_in, col_in, row_in)) continue;
303  //printf("a: local_in (%f,%f,%f) sec:%i, col:%i row:%i\n",local_in[0],local_in[1],local_in[2], sec_in, col_in, row_in);
304  if ((sec_in < 0) || (sec_in > 2)) {
305  cout << "sec_in:" << sec_in << endl;
306  continue;
307  }
308  fModuleInfo->TransformToLocalPad(local_in, x_in, y_in);
309  //printf("b: local_in (%f,%f,%f) sec:%i, col:%i row:%i (%f,%f)\n",local_in[0],local_in[1],local_in[2], sec_in, col_in, row_in, x_in, y_in);
310  row_in = fModuleInfo->GetModuleRow(sec_in, row_in);
311  //printf("c: local_in (%f,%f,%f) sec:%i, col:%i row:%i/%i (%f,%f)\n",local_in[0],local_in[1],local_in[2], sec_in, col_in, row_in, fModuleInfo->GetNofRows(), x_in, y_in);
312  if ((sec_out < 0) || (sec_out > 2)) {
313  cout << "sec_out:" << sec_out << endl;
314  continue;
315  }
316  //printf("d1: local_out (%f,%f,%f) sec:%i, col:%i row:%i\n",local_out[0],local_out[1],local_out[2], sec_out, col_out, row_out);
317  if (!fModuleInfo->GetPadInfo(local_out, sec_out, col_out, row_out))
318  continue;
319  //printf("d: local_out (%f,%f,%f) sec:%i, col:%i row:%i\n",local_out[0],local_out[1],local_out[2], sec_out, col_out, row_out);
320  fModuleInfo->TransformToLocalPad(local_out, x_out, y_out);
321  //printf("e: local_out (%f,%f,%f) sec:%i, col:%i row:%i (%f,%f)\n",local_out[0],local_out[1],local_out[2], sec_out, col_out, row_out, x_out, y_out);
322  row_out = fModuleInfo->GetModuleRow(sec_out, row_out);
323  //printf("f: local_out (%f,%f,%f) sec:%i, col:%i row:%i/%i (%f,%f)\n",local_out[0],local_out[1],local_out[2], sec_out, col_out,row_out, fModuleInfo->GetNofRows(), x_out, y_out);
324  Double_t W_in(fModuleInfo->GetPadSizeX(sec_in)),
325  W_out(fModuleInfo->GetPadSizeX(sec_out));
326  Double_t H_in(fModuleInfo->GetPadSizeY(sec_in)),
327  H_out(fModuleInfo->GetPadSizeY(sec_out));
328 
329 
330  pointCounter = fModuleMapPoint[moduleAddress]->GetN();
331  fModuleMapPoint[moduleAddress]->SetPoint(
332  pointCounter,
333  ((col_in + x_in / W_in) + (col_out + x_out / W_out)) / 2.,
334  ((row_in + y_in / H_in) + (row_out + y_out / H_out)) / 2.);
335  //pointCounter++;
336  fModuleMap[moduleAddress]->cd();
337  TLine* l = new TLine(col_in + x_in / W_in,
338  row_in + y_in / H_in,
339  col_out + x_out / W_out,
340  row_out + y_out / H_out);
341  l->SetLineWidth(2);
342  l->SetLineColor(15);
343  fModuleMapTrack[moduleAddress]->push_back(l);
344  fModuleMap[moduleAddress]->cd(1);
345  l->Draw("same");
346  /*
347  fModuleMap[moduleAddress]->cd(2);
348  l->Draw("same");
349  fModuleMap[moduleAddress]->cd(3);
350  l->Draw("same");
351  fModuleMap[moduleAddress]->cd(4);
352  l->Draw("same");
353  */
354  } else {
355  printf("Address:%i ID:%i\n", moduleAddress, moduleId);
356  break;
357  }
358  }
359  //fModuleMap[moduleAddress]->cd(1);
360  //name.Form("%i Points",pointCounter);
361  //TPaveText *textpoints = new TPaveText(10,2.5,50,3);
362  //textpoints->SetTextSize(0.035);
363  //textpoints->AddText(name);
364  //textpoints->Draw("same");
365  //cout << "Digis" << endl;
366  for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
367  digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(iDigi);
368  digiSpectrum->Fill(digi->GetCharge());
369  Int_t digiAddress = digi->GetAddress();
370  moduleAddress = CbmTrdAddress::GetModuleAddress(digiAddress);
371  fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(
372  moduleAddress);
373  Int_t sec(CbmTrdAddress::GetSectorId(digiAddress)),
374  row(CbmTrdAddress::GetRowId(digiAddress));
375  row = fModuleInfo->GetModuleRow(sec, row);
376  fModuleMapDigi[moduleAddress]->Fill(
377  CbmTrdAddress::GetColumnId(digiAddress), row, digi->GetCharge());
378  }
379  //cout << "Clusters" << endl;
380  Int_t lastModule(0), iCounter(0);
381  for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
382  Double_t charge(0), chargeMax(0);
383  cluster = (CbmTrdCluster*) fClusters->At(iCluster);
384  Int_t nDigisInCluster = cluster->GetNofDigis();
385  iCounter++;
386  for (Int_t iDigi = 0; iDigi < nDigisInCluster; iDigi++) {
387  digi =
388  CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(iDigi));
389  // digi = (CbmTrdDigi*)fDigis->At(cluster->GetDigi(iDigi));
390  Int_t digiAddress = digi->GetAddress();
391  moduleAddress = CbmTrdAddress::GetModuleAddress(digiAddress);
392  if (iDigi == 0 && lastModule != moduleAddress) {
393  iCounter = 0;
394  lastModule = moduleAddress;
395  }
396  //cout << moduleAddress << endl;
397  fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(
398  moduleAddress);
399  Int_t sec(CbmTrdAddress::GetSectorId(digiAddress)),
400  row(CbmTrdAddress::GetRowId(digiAddress));
401  row = fModuleInfo->GetModuleRow(sec, row);
402  fModuleMapCluster[moduleAddress]->Fill(
403  CbmTrdAddress::GetColumnId(digiAddress), row, iCounter + 1);
404  //fModuleMapCluster[moduleAddress]->SetBinContent(CbmTrdAddress::GetColumnId(digiAddress)+1, row+1, iCounter+1);
405  charge = digi->GetCharge();
406  if (charge > chargeMax) chargeMax = charge;
407  }
408  digiMaxSpectrum->Fill(chargeMax);
409  }
410 
411  //cout << "Hits" << endl;
412  Int_t modHit = 0;
413  for (Int_t iHit = 0; iHit < nHits; iHit++) {
414  hit = (CbmTrdHit*) fHits->At(iHit);
415 
416  moduleAddress = hit->GetAddress(); //GetDetectorID();//?????
417  fModuleInfo = (CbmTrdParModDigi*) fDigiPar->GetModulePar(
418  moduleAddress);
419  if (fModuleInfo == NULL) {
420  //printf("MA:%6i not found!\n",moduleAddress);
421  continue;
422  }
423 
424  Double_t pos[3] = {hit->GetX(), hit->GetY(), hit->GetZ()};
425  gGeoManager->FindNode(pos[0], pos[1], pos[2] - 0.3);
426  Double_t local_pos[3];
427  gGeoManager->MasterToLocal(pos, local_pos);
428  Int_t sec(-1), col(-1), row(-1);
429  fModuleInfo->GetPadInfo(local_pos, sec, col, row);
430  row = fModuleInfo->GetModuleRow(sec, row);
431  Double_t x(-1), y(-1);
432  fModuleInfo->TransformToLocalPad(local_pos, x, y);
433 
434  Double_t W(fModuleInfo->GetPadSizeX(sec)), H(fModuleInfo->GetPadSizeY(sec));
435  modHit = fModuleMapHit[moduleAddress]->GetN();
436  fModuleMapHit[moduleAddress]->SetPoint(modHit, col + x / W, row + y / H);
437 
438  if (hit->GetDx() <= W)
439  fModuleMapHit[moduleAddress]->SetPointError(
440  modHit, hit->GetDx() / W, hit->GetDy() / H);
441  else
442  fModuleMapHit[moduleAddress]->SetPointError(
443  modHit, hit->GetDy() / W, hit->GetDx() / H);
444  }
445 
446  gDirectory->mkdir("TrdRecoQA");
447  gDirectory->cd("TrdRecoQA");
448  digiMaxSpectrum->Write("", TObject::kOverwrite);
449  digiSpectrum->Write("", TObject::kOverwrite);
450  map<Int_t, TCanvas*>::iterator it;
451  for (it = fModuleMap.begin(); it != fModuleMap.end(); it++) {
452  it->second->cd(1);
453  name.Form("%i Points", fModuleMapPoint[it->first]->GetN());
454  fModuleMapPoint[it->first]->SetNameTitle(name, name);
455  fModuleMapPoint[it->first]->SetFillStyle(0);
456  TPaveText* ptext = new TPaveText(10, 3.5, 11, 4.0);
457  ptext->SetTextSize(0.035);
458  ptext->SetTextColor(15);
459  ptext->SetFillStyle(0);
460  ptext->SetLineColor(0);
461  ptext->AddText(name);
462  //fModuleMapPoint[it->first]->Draw("P");
463 
464 
465  //fModuleMapPoint[it->first]->Draw("P,same");
466  ptext->Draw("same");
467  it->second->cd(1)->Update();
468  it->second->cd(2);
469  fModuleMapDigi[it->first]->GetZaxis()->SetRangeUser(0, 0.0001);
470  fModuleMapDigi[it->first]->DrawCopy("colz");
471 
472  {
473  TPolyLine* pad = NULL;
474  const Int_t nRow = fModuleMapDigi[it->first]->GetNbinsY();
475  const Int_t nCol = fModuleMapDigi[it->first]->GetNbinsX();
476  const Double_t max_Range = fModuleMapDigi[it->first]->GetBinContent(
477  fModuleMapDigi[it->first]->GetMaximumBin());
478  for (Int_t iRow = 1; iRow <= nRow; iRow++) {
479  for (Int_t iCol = 1; iCol <= nCol; iCol++) {
480  Double_t charge =
481  fModuleMapDigi[it->first]->GetBinContent(iCol, iRow);
482  if (charge > 0.0) {
483  if (fTrianglePads)
484  pad = utils->CreateTriangularPad(
485  iCol - 1, iRow - 1, charge, 0.0, max_Range, false);
486  else
487  pad = utils->CreateRectangularPad(
488  iCol - 1, iRow - 1, charge, 0.0, max_Range, false);
489  pad->Draw("f,same");
490  }
491  }
492  }
493  }
494 
495  for (UInt_t t = 0; t < fModuleMapTrack[it->first]->size(); t++)
496  fModuleMapTrack[it->first]->at(t)->Draw("same");
497  /*
498  it->second->cd(3)->SetLogz(1);
499  fModuleMapDigi[it->first]->GetZaxis()->SetRangeUser(fTriggerTH,fModuleMapDigi[it->first]->GetBinContent(fModuleMapDigi[it->first]->GetMaximumBin()));
500  fModuleMapDigi[it->first]->DrawCopy("colz");
501  {
502  TPolyLine *pad = NULL;
503  const Int_t nRow = fModuleMapDigi[it->first]->GetNbinsY();
504  const Int_t nCol = fModuleMapDigi[it->first]->GetNbinsX();
505  const Double_t max_Range = fModuleMapDigi[it->first]->GetBinContent(fModuleMapDigi[it->first]->GetMaximumBin());
506  for (Int_t iRow = 1; iRow <= nRow; iRow++){
507  for (Int_t iCol = 1; iCol <= nCol; iCol++){
508  Double_t charge = fModuleMapDigi[it->first]->GetBinContent(iCol, iRow);
509  if (charge >= fTriggerTH){
510  if (fTrianglePads)
511  pad = utils->CreateTriangularPad(iCol-1, iRow-1, charge, 0, max_Range, true);
512  else
513  pad = utils->CreateRectangularPad(iCol-1, iRow-1, charge, 0.0, max_Range, false);
514  pad->Draw("f,same");
515  }
516  }
517  }
518  }
519  */
520 
521 
522  for (UInt_t t = 0; t < fModuleMapTrack[it->first]->size(); t++)
523  fModuleMapTrack[it->first]->at(t)->Draw("same");
524 
525  it->second->cd(3);
526  fModuleMapCluster[it->first]->DrawCopy("colz");
527  {
528  TPolyLine* pad = NULL;
529  const Int_t nRow = fModuleMapCluster[it->first]->GetNbinsY();
530  const Int_t nCol = fModuleMapCluster[it->first]->GetNbinsX();
531  const Double_t max_Range = fModuleMapCluster[it->first]->GetBinContent(
532  fModuleMapCluster[it->first]->GetMaximumBin());
533  for (Int_t iRow = 1; iRow <= nRow; iRow++) {
534  for (Int_t iCol = 1; iCol <= nCol; iCol++) {
535  Double_t clusterId =
536  fModuleMapCluster[it->first]->GetBinContent(iCol, iRow);
537  if (clusterId > 0) {
538  if (fTrianglePads)
539  pad = utils->CreateTriangularPad(
540  iCol - 1, iRow - 1, clusterId, 0, max_Range, false);
541  else
542  pad = utils->CreateRectangularPad(
543  iCol - 1, iRow - 1, clusterId, 0.0, max_Range, false);
544  pad->Draw("f,same");
545  }
546  }
547  }
548  }
549  for (UInt_t t = 0; t < fModuleMapTrack[it->first]->size(); t++)
550  fModuleMapTrack[it->first]->at(t)->Draw("same");
551  it->second->cd(3)->Update();
552  it->second->cd(4);
553  //fModuleMapPoint[it->first]->Draw();
554  //fModuleMapHit[it->first]->Draw("P,same");
555  if (fHits) {
556  /*
557  fModuleMapHit[it->first]->SetMaximum(fModuleMapDigi[it->first]->GetYaxis()->GetXmax());
558  fModuleMapHit[it->first]->SetMinimum(fModuleMapDigi[it->first]->GetYaxis()->GetXmin());
559  fModuleMapHit[it->first]->GetXaxis()->SetLimits(fModuleMapDigi[it->first]->GetXaxis()->GetXmin(),
560  fModuleMapDigi[it->first]->GetXaxis()->GetXmax());
561  */
562  name.Form("%i Hits", fModuleMapHit[it->first]->GetN());
563  fModuleMapHit[it->first]->SetNameTitle(name, name);
564  TPaveText* text = new TPaveText(10, 2.0, 11, 2.5);
565  text->SetFillStyle(0);
566  text->SetLineColor(0);
567  text->SetTextSize(0.035);
568  text->AddText(name);
569  TMultiGraph* mg = new TMultiGraph();
570 
571  //fModuleMapHit[it->first]->Draw("AP");
572  //fModuleMapPoint[it->first]->Draw("P,same");
573  mg->Draw("AP");
574  mg->Add(fModuleMapHit[it->first]);
575  mg->Add(fModuleMapPoint[it->first]);
576  mg->SetMaximum(fModuleMapDigi[it->first]->GetYaxis()->GetXmax());
577  mg->SetMinimum(fModuleMapDigi[it->first]->GetYaxis()->GetXmin());
578  mg->GetXaxis()->SetLimits(
579  fModuleMapDigi[it->first]->GetXaxis()->GetXmin(),
580  fModuleMapDigi[it->first]->GetXaxis()->GetXmax());
581  text->Draw("same");
582  ptext->Draw("same");
583  }
584  it->second->Update();
585  it->second->SaveAs("pics/" + (TString)(it->second->GetName()) + ".png");
586  it->second->Write("", TObject::kOverwrite);
587  it->second->Close();
588  }
589  gDirectory->cd("..");
590  timer.Stop();
591  Double_t rtime = timer.RealTime();
592  Double_t ctime = timer.CpuTime();
593 
594  printf("\n\n******************** Reading Test **********************\n");
595  printf(" RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime);
596  printf("*********************************************************\n\n");
597 }
598 
599 // --------------------------------------------------------------------------
600 //Int_t CbmTrdRecoQa::SecRowToGlobalRow(Int_t secRow)
601 //{
602 
603 //}
604 
605 // ---- Finish --------------------------------------------------------------
606 void CbmTrdRecoQa::Finish() { WriteHistograms(); }
607 // --------------------------------------------------------------------------
608 
609 // ---- Write test histograms ------------------------------------------------
610 
611 void CbmTrdRecoQa::WriteHistograms() {
612  /*
613  // Write histos to output
614  gDirectory->mkdir("TrdTracksPidQA");
615  gDirectory->cd("TrdTracksPidQA");
616 
617  if(WknPI) WknPI->Write();
618  if(WknEL) WknEL->Write();
619  if(WknALL) WknALL->Write();
620  if(WknLowPI) WknLowPI->Write();
621  if(WknLowEL) WknLowEL->Write();
622  if(WknLowALL) WknLowALL->Write();
623  if(WknHighPI) WknHighPI->Write();
624  if(WknHighEL) WknHighEL->Write();
625  if(WknHighALL) WknHighALL->Write();
626 
627  if(AnnPI) AnnPI->Write();
628  if(AnnEL) AnnEL->Write();
629  if(AnnALL) AnnALL->Write();
630  if(AnnLowPI) AnnLowPI->Write();
631  if(AnnLowEL) AnnLowEL->Write();
632  if(AnnLowALL) AnnLowALL->Write();
633  if(AnnHighPI) AnnHighPI->Write();
634  if(AnnHighEL) AnnHighEL->Write();
635  if(AnnHighALL) AnnHighALL->Write();
636 
637  if(LikePI) LikePI->Write();
638  if(LikeEL) LikeEL->Write();
639  if(LikeALL) LikeALL->Write();
640  if(LikeHighPI) LikeHighPI->Write();
641  if(LikeHighEL) LikeHighEL->Write();
642  if(LikeHighALL) LikeHighALL->Write();
643  if(LikeLowPI) LikeLowPI->Write();
644  if(LikeLowEL) LikeLowEL->Write();
645  if(LikeLowALL) LikeLowALL->Write();
646 
647  if(PartID) PartID->Write();
648  if(NrTRDHits) NrTRDHits->Write();
649  if(ELossPI) ELossPI->Write();
650  if(ELossEL) ELossEL->Write();
651  if(ELossALL) ELossALL->Write();
652  if(MomPI) MomPI->Write();
653  if(MomEL) MomEL->Write();
654  if(MomALL) MomALL->Write();
655  if(MOMvsELossPI) MOMvsELossPI->Write();
656  if(MOMvsELossEL) MOMvsELossEL->Write();
657  if(MOMvsELossALL) MOMvsELossALL->Write();
658 
659  gDirectory->cd("..");
660  */
661 }
662 
663 // --------------------------------------------------------------------------
664 
665 ClassImp(CbmTrdRecoQa);
CbmTrdRecoQa::Init
InitStatus Init()
Definition: CbmTrdRecoQa.cxx:83
CbmTrdParModDigi::GetPadSizeY
Double_t GetPadSizeY(Int_t i) const
Definition: CbmTrdParModDigi.h:37
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmTrdAddress::GetModuleAddress
static UInt_t GetModuleAddress(UInt_t address)
Return unique module ID from address.
Definition: CbmTrdAddress.h:117
CbmTrdPoint::GetYOut
Double_t GetYOut() const
Definition: CbmTrdPoint.h:67
CbmTrdDigi::GetAddress
Int_t GetAddress() const
Address getter for module in the format defined by CbmTrdDigi (format of CbmTrdAddress can be accesse...
Definition: CbmTrdDigi.h:92
CbmTrdRecoQa
Definition: CbmTrdRecoQa.h:42
CbmTrdPoint::GetZIn
Double_t GetZIn() const
Definition: CbmTrdPoint.h:65
CbmTrdUtils.h
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmTrdRecoQa::WriteHistograms
void WriteHistograms()
Definition: CbmTrdRecoQa.cxx:611
CbmTrdUtils
Definition: CbmTrdUtils.h:19
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmTrdRecoQa.h
CbmTrdParModDigi::GetNofRows
Int_t GetNofRows() const
Definition: CbmTrdParModDigi.cxx:340
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition: CbmDigiManager.cxx:71
rootMgr
static FairRootManager * rootMgr
Definition: CbmDeviceHitBuilderTof.cxx:72
CbmTrdRecoQa::CbmTrdRecoQa
CbmTrdRecoQa()
Definition: CbmTrdRecoQa.cxx:53
CbmTrdAddress::GetSectorId
static UInt_t GetSectorId(UInt_t address)
Return sector ID from address.
Definition: CbmTrdAddress.h:88
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmTrdHit
data class for a reconstructed Energy-4D measurement in the TRD
Definition: CbmTrdHit.h:35
CbmTrdPoint::GetXOut
Double_t GetXOut() const
Definition: CbmTrdPoint.h:66
CbmTrdAddress::GetModuleId
static UInt_t GetModuleId(UInt_t address)
Return module ID from address.
Definition: CbmTrdAddress.h:78
CbmTrdRecoQa::SetParContainers
void SetParContainers()
Definition: CbmTrdRecoQa.cxx:129
CbmDigiManager::GetNofDigis
static Int_t GetNofDigis(ECbmModuleId systemId)
Definition: CbmDigiManager.cxx:62
CbmTrdRecoQa::fModuleMapPoint
std::map< Int_t, TGraphErrors * > fModuleMapPoint
Definition: CbmTrdRecoQa.h:89
CbmTrdRecoQa::fModuleMapDigi
std::map< Int_t, TH2D * > fModuleMapDigi
Definition: CbmTrdRecoQa.h:90
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmTrdRecoQa::ReInit
InitStatus ReInit()
Definition: CbmTrdRecoQa.cxx:142
CbmTrdParModDigi::GetPadInfo
Bool_t GetPadInfo(const Double_t *local_point, Int_t &sectorId, Int_t &columnId, Int_t &rowId) const
Definition: CbmTrdParModDigi.cxx:436
CbmTrdGeoHandler
Definition: CbmTrdGeoHandler.h:29
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition: CbmDigiManager.h:93
CbmCluster::GetNofDigis
Int_t GetNofDigis() const
Number of digis in cluster.
Definition: CbmCluster.h:69
CbmTrdRecoQa::fModuleMap
std::map< Int_t, TCanvas * > fModuleMap
Definition: CbmTrdRecoQa.h:88
CbmTrdUtils::CreateRectangularPad
TPolyLine * CreateRectangularPad(Int_t column, Int_t row, Double_t value, Double_t min_range, Double_t max_range, Bool_t logScale)
Definition: CbmTrdUtils.cxx:101
CbmTrdCluster
Data Container for TRD clusters.
Definition: CbmTrdCluster.h:23
CbmTrdParModDigi::GetPadSizeX
Double_t GetPadSizeX(Int_t i) const
Definition: CbmTrdParModDigi.h:36
CbmTrdDigi.h
CbmTrdRecoQa::fDigiPar
CbmTrdParSetDigi * fDigiPar
Definition: CbmTrdRecoQa.h:75
CbmTrdRecoQa::SetTriggerThreshold
void SetTriggerThreshold(Double_t triggerthreshold)
Definition: CbmTrdRecoQa.cxx:78
CbmTrdRecoQa::fMCPoints
TClonesArray * fMCPoints
Definition: CbmTrdRecoQa.h:73
CbmTrdGeoHandler.h
Helper class to extract information from the GeoManager.
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmTrdRecoQa::fClusters
TClonesArray * fClusters
Definition: CbmTrdRecoQa.h:71
ClassImp
ClassImp(CbmTrdRecoQa)
CbmTrdParModDigi
Definition of chamber gain conversion for one TRD module.
Definition: CbmTrdParModDigi.h:14
CbmTrdGeoHandler::Init
void Init(Bool_t isSimulation=kFALSE)
Definition: CbmTrdGeoHandler.cxx:45
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
CbmTrdParModDigi.h
CbmTrdRecoQa::fHits
TClonesArray * fHits
Definition: CbmTrdRecoQa.h:72
CbmTrdRecoQa::~CbmTrdRecoQa
virtual ~CbmTrdRecoQa()
Definition: CbmTrdRecoQa.cxx:76
CbmTrdParSet::GetModulePar
virtual const CbmTrdParMod * GetModulePar(Int_t detId) const
Definition: CbmTrdParSet.cxx:45
CbmTrdHit.h
Class for hits in TRD detector.
CbmTrdRecoQa::fTriggerTH
Double_t fTriggerTH
Definition: CbmTrdRecoQa.h:69
CbmTrdPoint::GetZOut
Double_t GetZOut() const
Definition: CbmTrdPoint.h:68
CbmTrdParSetDigi.h
CbmTrdPoint
Definition: CbmTrdPoint.h:23
CbmTrdRecoQa::fGeoHandler
CbmTrdGeoHandler * fGeoHandler
Definition: CbmTrdRecoQa.h:77
ECbmModuleId::kTrd
@ kTrd
Transition Radiation Detector.
CbmTrdRecoQa::fTrianglePads
Bool_t fTrianglePads
Definition: CbmTrdRecoQa.h:68
CbmTrdParModDigi::TransformToLocalPad
void TransformToLocalPad(const Double_t *local_point, Double_t &posX, Double_t &posY) const
Definition: CbmTrdParModDigi.cxx:634
CbmTrdRecoQa::SetTriangularPads
void SetTriangularPads(Bool_t triangles)
Definition: CbmTrdRecoQa.cxx:154
CbmMCTrack.h
CbmTrdPoint::GetXIn
Double_t GetXIn() const
Definition: CbmTrdPoint.h:63
CbmDigiManager.h
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTrdRecoQa::fModuleMapHit
std::map< Int_t, TGraphErrors * > fModuleMapHit
Definition: CbmTrdRecoQa.h:92
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTrdPoint.h
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmTrdDigi
Definition: CbmTrdDigi.h:14
CbmTrdRecoQa::Finish
virtual void Finish()
Definition: CbmTrdRecoQa.cxx:606
CbmTrdParSetDigi
Definition: CbmTrdParSetDigi.h:15
CbmTrdRecoQa::fModuleMapCluster
std::map< Int_t, TH2I * > fModuleMapCluster
Definition: CbmTrdRecoQa.h:91
CbmTrdCluster.h
Data Container for TRD clusters.
CbmTrdParModDigi::GetModuleRow
Int_t GetModuleRow(Int_t &sectorId, Int_t &rowId) const
Definition: CbmTrdParModDigi.cxx:373
CbmTrdDigi::GetCharge
Double_t GetCharge() const
Charge getter for SPADIC.
Definition: CbmTrdDigi.cxx:133
CbmTrdPoint::GetYIn
Double_t GetYIn() const
Definition: CbmTrdPoint.h:64
CbmCluster::GetDigi
Int_t GetDigi(Int_t index) const
Get digi at position index.
Definition: CbmCluster.h:76
CbmTrdRecoQa::Exec
virtual void Exec(Option_t *option)
Definition: CbmTrdRecoQa.cxx:172
CbmTrdRecoQa::fModuleMapTrack
std::map< Int_t, std::vector< TLine * > * > fModuleMapTrack
Definition: CbmTrdRecoQa.h:93
CbmTrdParModDigi::GetNofColumns
Int_t GetNofColumns() const
Definition: CbmTrdParModDigi.cxx:321
CbmTrdAddress::GetRowId
static UInt_t GetRowId(UInt_t address)
Return row ID from address.
Definition: CbmTrdAddress.h:98
CbmTrdRecoQa::fModuleInfo
CbmTrdParModDigi * fModuleInfo
Definition: CbmTrdRecoQa.h:76
CbmTrdAddress::GetColumnId
static UInt_t GetColumnId(UInt_t address)
Return column ID from address.
Definition: CbmTrdAddress.h:107
CbmTrdUtils::CreateTriangularPad
TPolyLine * CreateTriangularPad(Int_t column, Int_t row, Double_t value, Double_t min_range, Double_t max_range, Bool_t logScale)
Definition: CbmTrdUtils.cxx:66