CbmRoot
CbmMuchHitFinderQa.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- CbmMuchHitProducerQa source file -----
3 // ----- Modified since 21/06/2019 by Ekata Nandy (ekata@vecc.gov.in) -Inclusion of RPC detector type
4 // ----- Modified 02/18 by Vikas Singhal -----
5 // ----- Created 16/11/07 by E. Kryshen -----
6 // -------------------------------------------------------------------------
7 // AUG 18 RELEASE VERSION
8 #include "CbmMuchHitFinderQa.h"
9 #include "CbmMuchCluster.h"
10 #include "CbmMuchDigi.h"
11 #include "CbmMuchPixelHit.h"
12 #include "CbmMuchPoint.h"
13 #include "FairRootManager.h"
14 
15 #include "CbmDigiManager.h"
16 #include "CbmMuchAddress.h"
17 #include "CbmMuchDigi.h"
18 #include "CbmMuchModuleGem.h"
19 #include "CbmMuchPad.h"
20 #include "CbmMuchSector.h"
21 #include "CbmMuchStation.h"
22 
23 #include "CbmMCTrack.h"
24 #include "CbmMatch.h"
25 #include "FairLogger.h"
26 #include "TDatabasePDG.h"
27 #include "TParticlePDG.h"
28 
29 #include "CbmMuchPointInfo.h"
30 #include "CbmMuchRecoDefs.h"
31 #include "TCanvas.h"
32 #include "TF1.h"
33 #include "TFile.h"
34 
35 #include "CbmMuchGeoScheme.h"
36 #include "TArrayI.h"
37 #include "TObjArray.h"
38 #include "TStyle.h"
39 
40 
41 #include "CbmGeoMuchPar.h"
42 #include "FairRuntimeDb.h"
43 #include "TGraph.h"
44 
45 #include <algorithm>
46 #include <cassert>
47 #include <map>
48 #include <vector>
49 
50 using std::cout;
51 using std::endl;
52 using std::map;
53 using std::vector;
54 // -------------------------------------------------------------------------
55 CbmMuchHitFinderQa::CbmMuchHitFinderQa(const char* name, Int_t verbose)
56  : FairTask(name, verbose)
57  , fGeoScheme(CbmMuchGeoScheme::Instance())
58  , fGeoFileName()
59  , fFileName()
60  , fSignalPoints(0)
61  , fSignalHits(0)
62  , fVerbose(verbose)
63  , fEvent(0)
64  , fFlag(0)
65  , fPoints(NULL)
66  ,
67  //fDigis(NULL),
68  //fDigiMatches(NULL),
69  fClusters(NULL)
70  , fHits(NULL)
71  , fMCTracks(NULL)
72  , fPointInfos(new TClonesArray("CbmMuchPointInfo", 10))
73  , fNstations(0)
74  , fChargeHistos(NULL)
75  , fhChargeEnergyLog(NULL)
76  , fhChargeEnergyLogPi(NULL)
77  , fhChargeEnergyLogPr(NULL)
78  , fhChargeEnergyLogEl(NULL)
79  , fhChargeTrackLength(NULL)
80  , fhChargeTrackLengthPi(NULL)
81  , fhChargeTrackLengthPr(NULL)
82  , fhChargeTrackLengthEl(NULL)
83  , fhChargeLog(NULL)
84  , fhChargePr_1GeV_3mm(NULL)
85  , fhNpadsVsS(NULL)
86  , fhCharge(NULL)
87  , fhOccupancyR(NULL)
88  , fhPadsTotalR(NULL)
89  , fhPadsFiredR(NULL)
90  , fhPullXpads1(NULL)
91  , fhPullYpads1(NULL)
92  , fhPullXpads2(NULL)
93  , fhPullYpads2(NULL)
94  , fhPullXpads3(NULL)
95  , fhPullYpads3(NULL)
96  , fnPadSizesX(0)
97  , fnPadSizesY(0)
98  , fNTimingPulls(8)
99  , fhPullX(NULL)
100  , fhPullY(NULL)
101  , fhPullT(NULL)
102  , fhResidualX(NULL)
103  , fhResidualY(NULL)
104  , fhResidualT(NULL)
105  , fhPointsInCluster(NULL)
106  , fhDigisInCluster(NULL)
107  , fhHitsInCluster(NULL)
108  , fNall(NULL)
109  , fNpr(NULL)
110  , fNpi(NULL)
111  , fNel(NULL)
112  , fNmu(NULL)
113  , fNka(NULL)
114  , fNprimary(NULL)
115  , fNsecondary(NULL)
116  , fPointsTotal(0)
117  , fPointsUnderCounted(0)
118  , fPointsOverCounted(0)
119  , fOccupancyQaOn(kTRUE)
120  , fPullsQaOn(kTRUE)
121  , fDigitizerQaOn(kTRUE)
122  , fStatisticsQaOn(kTRUE)
123  , fClusterDeconvQaOn(kTRUE)
124  , fPrintToFileOn(kTRUE)
125  , fPadMinLx(0.)
126  , fPadMinLy(0.)
127  , fPadMaxLx(0.)
128  , fPadMaxLy(0.)
129  , pointsFile(NULL)
130  , padsFile(NULL)
131 
132 {}
133 // -------------------------------------------------------------------------
134 
135 
136 // -------------------------------------------------------------------------
138 // -------------------------------------------------------------------------
139 
140 
141 // -------------------------------------------------------------------------
143  FairRootManager* fManager = FairRootManager::Instance();
144  fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack");
145  fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
146  fHits = (TClonesArray*) fManager->GetObject("MuchPixelHit");
147  //fDigis = (TClonesArray*) fManager->GetObject("MuchDigi");
148  //fDigiMatches = (TClonesArray*) fManager->GetObject("MuchDigiMatch"); ///added ekata
149  fClusters = (TClonesArray*) fManager->GetObject("MuchCluster");
150  // Reading Much Digis from CbmMuchDigiManager which are stored as vector
152  fDigiManager->Init();
153 
154  // printf(" %i",fMCTracks);
155  // printf(" %i",fPoints);
156  // printf(" %i",fHits);
157  // printf(" %i",fDigis);
158  // printf(" %i",fDigiMatches);
159  // printf(" %i",fClusters);
160  // printf("\n");
161 
162  TFile* f = new TFile(fGeoFileName, "R");
163  TObjArray* stations = (TObjArray*) f->Get("stations");
164  fGeoScheme->Init(stations, fFlag);
166  printf("Init: fNstations = %i\n", fNstations);
167 
168  fNall = new Int_t[fNstations];
169  fNpr = new Int_t[fNstations];
170  fNpi = new Int_t[fNstations];
171  fNel = new Int_t[fNstations];
172  fNmu = new Int_t[fNstations];
173  fNka = new Int_t[fNstations];
174  fNprimary = new Int_t[fNstations];
175  fNsecondary = new Int_t[fNstations];
176 
177  // Reset counters
178  for (Int_t i = 0; i < fNstations; i++) {
179  fNall[i] = 0;
180  fNpr[i] = 0;
181  fNpi[i] = 0;
182  fNel[i] = 0;
183  fNmu[i] = 0;
184  fNka[i] = 0;
185  fNprimary[i] = 0;
186  fNsecondary[i] = 0;
187  }
188 
189  fPointsTotal = 0;
191  fPointsOverCounted = 0;
192 
193 #define BINNING_CHARGE 100, 0, 3.0
194 #define BINNING_LENGTH 100, 0, 2.5
195 #define BINNING_CHARGE_LOG 100, 4, 8
196 #define BINNING_ENERGY_LOG 100, -0.5, 4.5
197 #define BINNING_ENERGY_LOG_EL 100, -0.5, 4.5
198 
199  // fhCharge = new TH1D("hCharge",
200  // "Charge distribution from tracks",
201  // BINNING_CHARGE);
202 
203  fhChargeLog = new TH1D(
204  "hChargeLog", "Charge (log.) distribution from tracks", BINNING_CHARGE_LOG);
205 
206  fhChargeEnergyLog = new TH2D("hChargeEnergyLog",
207  "Charge vs energy (log.) for all tracks",
210 
211  fhChargeEnergyLogPi = new TH2D("hChargeEnergyLogPi",
212  "Charge vs energy (log.) for pions",
215 
216  fhChargeEnergyLogPr = new TH2D("hChargeEnergyLogPr",
217  "Charge vs energy (log.) for protons",
220 
221  fhChargeEnergyLogEl = new TH2D("hChargeEnergyLogEl",
222  "Charge vs energy (log.) for electrons",
225 
226  fhChargeTrackLength = new TH2D("hChargeTrackLength",
227  "Charge vs length for all tracks",
230 
231  fhChargeTrackLengthPi = new TH2D("hChargeTrackLengthPi",
232  "Charge vs length for pions",
235  fhChargeTrackLengthPr = new TH2D("hChargeTrackLengthPr",
236  "Charge vs length for proton",
239 
240  fhChargeTrackLengthEl = new TH2D("hChargeTrackLengthEl",
241  "Charge vs length for electrons",
244 
246  new TH1D("hChargePr_1GeV_3mm", "Charge for 1 GeV protons", BINNING_CHARGE);
247 
248  // fhPointsInCluster = new TH1I("hPointsInCluster",
249  // "Number of tracks in cluster",
250  // 10,0,10);
251 
252  // fhDigisInCluster = new TH1I("hDigisInCluster",
253  // "Number of digis in cluster",
254  // 10,0,10);
255 
256  // fhHitsInCluster = new TH1I("hHitsInCluster",
257  // "Number of hits in cluster",
258  // 10,0,10);
259 
260 
261  // fhCharge ->GetXaxis()->SetTitle("Charge [10^{6} electrons]");
262  fhChargeLog->GetXaxis()->SetTitle("Lg (Charge) [Number of electrons]");
263  fhChargePr_1GeV_3mm->GetXaxis()->SetTitle("Charge [10^{6} electrons]");
264 
265 
266  fChargeHistos = new TObjArray();
275  // fChargeHistos->Add(fhCharge);
277 
278  for (Int_t i = 0; i < 8; i++) {
279  TH2D* histo = (TH2D*) fChargeHistos->At(i);
280  histo->SetStats(0);
281  histo->GetYaxis()->SetDecimals(kTRUE);
282  histo->GetYaxis()->SetTitleOffset(1.4);
283  histo->GetYaxis()->CenterTitle();
284  histo->GetYaxis()->SetTitle("Charge [10^{6} electrons]");
285  if (i < 4)
286  histo->GetXaxis()->SetTitle("Lg E_{kin} [MeV]");
287  else
288  histo->GetXaxis()->SetTitle("Track length [cm]");
289  }
290 
291  fhPointsInCluster = new TH1I*[fNstations];
292  fhDigisInCluster = new TH1I*[fNstations];
293  fhHitsInCluster = new TH1I*[fNstations];
294 
295  fhCharge = new TH1D*[fNstations];
296  fhOccupancyR = new TH1D*[fNstations];
297  fhPadsTotalR = new TH1D*[fNstations];
298  fhPadsFiredR = new TH1D*[fNstations];
299 
300  for (Int_t i = 0; i < fNstations; i++) {
301  CbmMuchStation* station = fGeoScheme->GetStation(i);
302  Double_t rMax = station->GetRmax();
303  Double_t rMin = station->GetRmin();
304  fhPadsTotalR[i] =
305  new TH1D(Form("hPadsTotal%i", i + 1),
306  Form("Number of pads vs radius: station %i;Radius [cm]", i + 1),
307  100,
308  0.6 * rMin,
309  1.2 * rMax);
310  fhPadsFiredR[i] = new TH1D(
311  Form("hPadsFired%i", i + 1),
312  Form("Number of fired pads vs radius: station %i;Radius [cm]", i + 1),
313  100,
314  0.6 * rMin,
315  1.2 * rMax);
316  fhOccupancyR[i] = new TH1D(
317  Form("hOccupancy%i", i + 1),
318  Form("Occupancy vs radius: station %i;Radius [cm];Occupancy, %%", i + 1),
319  100,
320  0.6 * rMin,
321  1.2 * rMax);
322 
323  fhCharge[i] =
324  new TH1D(Form("hCharge%i", i + 1),
325  Form("Charge : station %i; Charge(1e4); Count", i + 1),
326  200,
327  0,
328  500);
330  new TH1I(Form("hPointsInCluster%i", i + 1),
331  Form("Points in Cluster : Station %i ", i + 1),
332  10,
333  0,
334  10);
336  new TH1I(Form("hDigisInCluster%i", i + 1),
337  Form("Digis in Cluster : Station %i ", i + 1),
338  10,
339  0,
340  10);
341  fhHitsInCluster[i] = new TH1I(Form("hHitsInCluster%i", i + 1),
342  Form("Hits in Cluster : Station %i ", i + 1),
343  10,
344  0,
345  10);
346  }
347 
348  vector<CbmMuchModule*> modules = fGeoScheme->GetModules();
349  for (vector<CbmMuchModule*>::iterator im = modules.begin();
350  im != modules.end();
351  im++) {
352  CbmMuchModule* mod = (*im);
353  if (mod->GetDetectorType() == 4
354  || mod->GetDetectorType() == 3) { // modified for rpc
355  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
356  vector<CbmMuchPad*> pads = module->GetPads();
357  for (UInt_t ip = 0; ip < pads.size(); ip++) {
358  CbmMuchPad* pad = pads[ip];
359  Int_t stationId = CbmMuchAddress::GetStationIndex(pad->GetAddress());
360  Double_t x0 = pad->GetX();
361  Double_t y0 = pad->GetY();
362  Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
363  fhPadsTotalR[stationId]->Fill(r0);
364  }
365  }
366  }
367  /*
368  fPadMinLx = std::numeric_limits<Double_t>::max();
369  fPadMinLy = std::numeric_limits<Double_t>::max();
370  fPadMaxLx = std::numeric_limits<Double_t>::min();
371  fPadMaxLy = std::numeric_limits<Double_t>::min();
372 
373  Int_t nTotSectors = 0;
374  Int_t nTotChannels = 0;
375  printf("=========================================================================================\n");
376  printf(" Station Nr.\t| Sectors\t| Channels\t| Pad min size\t\t| Pad max length\t \n");
377  printf("-----------------------------------------------------------------------------------------\n");
378  */
379  Int_t nTotChannels = 0;
380  for (Int_t iStation = 0; iStation < fNstations; iStation++) {
381  /*
382  Double_t padMinLx = GetMinPadSize(iStation).X();
383  Double_t padMinLy = GetMinPadSize(iStation).Y();
384  Double_t padMaxLx = GetMaxPadSize(iStation).X();
385  Double_t padMaxLy = GetMaxPadSize(iStation).Y();
386  */
387  Int_t nChannels = GetNChannels(iStation);
388  /*
389  Int_t nSectors = GetNSectors(iStation);
390  if (fPadMinLx>padMinLx) fPadMinLx = padMinLx;
391  if (fPadMinLy>padMinLy) fPadMinLy = padMinLy;
392  if (fPadMaxLx<padMaxLx) fPadMaxLx = padMaxLx;
393  if (fPadMaxLy<padMaxLy) fPadMaxLy = padMaxLy;
394  printf("%i\t\t| %i\t\t| %i\t| %5.4fx%5.4f\t\t| %5.4fx%5.4f\n",iStation+1, nSectors, nChannels, padMinLx, padMinLy, padMaxLx, padMaxLy);
395  */
396  // nTotSectors += nSectors;
397  printf("%i\t\t| %i\t\t\n", iStation + 1, nChannels);
398  nTotChannels += nChannels;
399  }
400  printf("---------------------------------------------------------------------"
401  "--------------------\n");
402  printf(" Total:\t\t| %i\t\t\n", nTotChannels);
403  printf("====================================================================="
404  "====================\n");
405  /*
406  fnPadSizesX = TMath::CeilNint(TMath::Log2(fPadMaxLx/fPadMinLx)+1);
407  fnPadSizesY = TMath::CeilNint(TMath::Log2(fPadMaxLy/fPadMinLy)+1);
408  Info("Init","nPadSizesX=%i",fnPadSizesX);
409  Info("Init","nPadSizesY=%i",fnPadSizesY);
410  fhPullXpads1 = new TH1D*[fnPadSizesX];
411  fhPullYpads1 = new TH1D*[fnPadSizesY];
412  fhPullXpads2 = new TH1D*[fnPadSizesX];
413  fhPullYpads2 = new TH1D*[fnPadSizesY];
414  fhPullXpads3 = new TH1D*[fnPadSizesX];
415  fhPullYpads3 = new TH1D*[fnPadSizesY];
416  fhPullT = new TH1D*[fNTimingPulls];
417 
418  for (Int_t i=0;i<fnPadSizesX;i++){
419  fhPullXpads1[i] = new TH1D(Form("hPullXpads1%i",i),Form("Pull distribution X. Npads = 1 Size =%i; (x_{RC} - x_{MC}) / dx_{RC}",i),100,-5,5);
420  fhPullXpads2[i] = new TH1D(Form("hPullXpads2%i",i),Form("Pull distribution X. Npads = 2 Size =%i; (x_{RC} - x_{MC}) / dx_{RC}",i),100,-5,5);
421  fhPullXpads3[i] = new TH1D(Form("hPullXpads3%i",i),Form("Pull distribution X. Npads = 3 Size =%i; (x_{RC} - x_{MC}) / dx_{RC}",i),100,-5,5);
422  }
423 
424  for (Int_t i=0;i<fnPadSizesY;i++){
425  fhPullYpads1[i] = new TH1D(Form("hPullYpads1%i",i),Form("Pull distribution Y. Npads = 1 Size =%i; (y_{RC} - y_{MC}) / dy_{RC}",i),100,-5,5);
426  fhPullYpads2[i] = new TH1D(Form("hPullYpads2%i",i),Form("Pull distribution Y. Npads = 2 Size =%i; (y_{RC} - y_{MC}) / dy_{RC}",i),100,-5,5);
427  fhPullYpads3[i] = new TH1D(Form("hPullYpads3%i",i),Form("Pull distribution Y. Npads = 3 Size =%i; (y_{RC} - y_{MC}) / dy_{RC}",i),100,-5,5);
428  }
429 
430  fhPullT[0] = new TH1D("hPullT","Pull distribution T, all pads; (t_{RC} - t_{MC}) / dt_{RC}",100,-5,5);
431  for (Int_t i=1;i<fNTimingPulls;i++){
432  fhPullT[i] = new TH1D(Form("hPullT%i",i),Form("Pull distribution T. Npads = %i; (t_{RC} - t_{MC}) / dt_{RC}",i),100,-5,5);
433  }
434  */
435  //Pull Distribution
436  fhPullX = new TH1D(
437  "hPullX", "Pull distribution X;(x_{RC} - x_{MC}) / dx_{RC}", 500, -5, 5);
438  fhPullY = new TH1D(
439  "hPullY", "Pull distribution Y;(y_{RC} - y_{MC}) / dy_{RC}", 500, -5, 5);
440  fhPullT = new TH1D(
441  "hPullT", "Pull distribution T;(t_{RC} - t_{MC}) / dt_{RC}", 120, -10, 50);
442 
443  //Residual Distribution
444  fhResidualX = new TH1D(
445  "hResidualX", "Residual distribution X;(x_{RC} - x_{MC})(cm)", 500, -5, 5);
446  fhResidualY = new TH1D(
447  "hResidualY", "Residual distribution Y;(y_{RC} - y_{MC})(cm)", 500, -5, 5);
448  fhResidualT = new TH1D("hResidualT",
449  "Residual distribution T;(t_{RC} - t_{MC})(ns)",
450  140,
451  -20,
452  50);
453 
454  fhNpadsVsS = new TH2D("hNpadsVsS",
455  "Number of fired pads vs pad area:area:n pads",
456  10,
457  -5,
458  0,
459  10,
460  0.5,
461  10.5);
462  pointsFile = fopen("points.txt", "w");
463  padsFile = fopen("pads.txt", "w");
464  return kSUCCESS;
465 }
466 // -------------------------------------------------------------------------
467 
468 
469 // -------------------------------------------------------------------------
471  // Get run and runtime database
472  // FairRuntimeDb* db = FairRuntimeDb::instance();
473  // if ( ! db ) Fatal("SetParContainers", "No runtime database");
474  // Get MUCH geometry parameter container
475  // fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
476 }
477 // -------------------------------------------------------------------------
478 
479 
480 // -------------------------------------------------------------------------x
481 void CbmMuchHitFinderQa::Exec(Option_t*) {
482  fEvent++;
483  LOG(info) << "Event: " << fEvent;
484  fprintf(pointsFile, "Event %i\n", fEvent);
485  fprintf(padsFile, "Event %i\n", fEvent);
486 
487  if (fPullsQaOn) PullsQa();
492 
493  return;
494  for (int i = 0; i < fPoints->GetEntriesFast(); i++) {
495  CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(i);
496  Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
497  if (stId != 0) continue;
498  Int_t layerId = CbmMuchAddress::GetLayerIndex(point->GetDetectorID());
499  if (layerId != 0) continue;
500  printf("point %4i xin=%6.2f yin=%6.2f xout=%6.2f yout=%6.2f zin=%6.2f\n",
501  i,
502  point->GetXIn(),
503  point->GetYIn(),
504  point->GetXOut(),
505  point->GetYOut(),
506  point->GetZIn());
507  fprintf(pointsFile,
508  "%7.3f %7.3f %7.3f %7.3f\n",
509  point->GetXIn(),
510  point->GetYIn(),
511  point->GetXOut(),
512  point->GetYOut());
513  }
514 
515  // old code
516  /* for (Int_t i=0;i<fDigis->GetEntriesFast();i++){
517  CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(i); */
518  for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
520  UInt_t address = digi->GetAddress();
521  Int_t stId = CbmMuchAddress::GetStationIndex(address);
522  if (stId != 0) continue;
523  Int_t layerId = CbmMuchAddress::GetLayerIndex(address);
524  if (layerId != 0) continue;
525  CbmMuchModuleGem* module =
527  if (!module) continue;
528  CbmMuchPad* pad = module->GetPad(address);
529  Double_t x0 = pad->GetX();
530  Double_t y0 = pad->GetY();
531  UInt_t charge = digi->GetAdc();
532  printf("digi %4i x0=%5.1f y0=%5.1f charge=%3i\n", i, x0, y0, charge);
533  fprintf(padsFile, "%5.1f %5.1f %3i\n", x0, y0, charge);
534  }
535 }
536 // -------------------------------------------------------------------------
537 
538 
539 // -------------------------------------------------------------------------
541  printf("FinishTask\n");
542  fclose(pointsFile);
543  fclose(padsFile);
544 
545  gStyle->SetPaperSize(20, 20);
546  // Int_t NuOfColoums = 1;
547  // Int_t NuOfRows = 1;
548 
549  //Setting Canvas according to the Number of stations
550  /*
551  if(fNstations>4){ NuOfColoums = 3; NuOfRows = 2;}
552  else{
553  if(fNstations<=4 && fNstations >=3){ NuOfColoums = 2; NuOfRows = 2;}
554  else { NuOfColoums = fNstations; NuOfRows = 1;}
555  }
556 */
557 
558  Double_t errors[100];
559  for (Int_t i = 0; i < 100; i++)
560  errors[i] = 0;
561 
562  for (Int_t i = 0; i < fNstations; i++) {
563  fhPadsTotalR[i]->Sumw2();
564  fhPadsTotalR[i]->SetError(errors);
565  fhPadsFiredR[i]->Sumw2();
566  fhPadsFiredR[i]->Scale(1. / fEvent);
567  fhOccupancyR[i]->Divide(fhPadsFiredR[i], fhPadsTotalR[i]);
568  fhOccupancyR[i]->Scale(100);
569  fhCharge[i]->Scale(1. / fEvent);
570  fhPointsInCluster[i]->Scale(1. / fEvent);
571  fhDigisInCluster[i]->Scale(1. / fEvent);
572  fhHitsInCluster[i]->Scale(1. / fEvent);
573  }
574 
575  if (fPullsQaOn && fVerbose > 1) {
576  printf("===================================\n");
577  printf("PullsQa:\n");
578 
579  // TCanvas* cTiming = new TCanvas("cTiming","Timing pulls",250*4,250*2);
580  // cTiming->Divide(4,2);
581  // for (Int_t i=0;i<fNTimingPulls;i++){
582  // cTiming->cd(i+1);
583  // fhPullT[i]->Fit("gaus");
584  // fhPullT[i]->Draw("e");
585  // }
586  //
587  /* TCanvas* c4 = new TCanvas("c4","Pulls",800,400);
588  c4->Divide(3,1);
589  c4->cd(1);*/
590  fhPullX->Sumw2();
591  fhPullX->Fit("gaus", "Q");
592  fhPullX->GetFunction("gaus")->SetLineWidth(3);
593  fhPullX->GetFunction("gaus")->SetLineColor(kRed);
594  // fhPullX->Draw();
595  fhPullX->Write();
596  if (fPrintToFileOn) gPad->Print(".gif");
597  if (fPrintToFileOn) gPad->Print(".eps");
598  // c4->cd(2);
599  fhPullY->Sumw2();
600  fhPullY->Fit("gaus", "Q");
601  fhPullY->GetFunction("gaus")->SetLineWidth(3);
602  fhPullY->GetFunction("gaus")->SetLineColor(kRed);
603  // fhPullY->Draw();
604  fhPullY->Write();
605  if (fPrintToFileOn) gPad->Print(".gif");
606  if (fPrintToFileOn) gPad->Print(".eps");
607  // c4->cd(3);
608  fhPullT->Sumw2();
609  fhPullT->Fit("gaus", "Q");
610  fhPullT->GetFunction("gaus")->SetLineWidth(3);
611  fhPullT->GetFunction("gaus")->SetLineColor(kRed);
612  // fhPullT->Draw();
613  fhPullT->Write();
614  if (fPrintToFileOn) gPad->Print(".gif");
615  if (fPrintToFileOn) gPad->Print(".eps");
616  // c4->cd();
617 
618  // TCanvas* cR5 = new TCanvas("cR5","Residuals",800,400);
619  // cR5->Divide(3,1);
620  // cR5->cd(1);
621  fhResidualX->Sumw2();
622  fhResidualX->Fit("gaus", "Q");
623  fhResidualX->GetFunction("gaus")->SetLineWidth(3);
624  fhResidualX->GetFunction("gaus")->SetLineColor(kRed);
625  // fhResidualX->Draw();
626  fhResidualX->Write();
627  // if (fPrintToFileOn) gPad->Print(".gif");
628  // if (fPrintToFileOn) gPad->Print(".eps");
629  // cR5->cd(2);
630  fhResidualY->Sumw2();
631  fhResidualY->Fit("gaus", "Q");
632  fhResidualY->GetFunction("gaus")->SetLineWidth(3);
633  fhResidualY->GetFunction("gaus")->SetLineColor(kRed);
634  // fhResidualY->Draw();
635  fhResidualY->Write();
636  // if (fPrintToFileOn) gPad->Print(".gif");
637  // if (fPrintToFileOn) gPad->Print(".eps");
638  // cR5->cd(3);
639  fhResidualT->Sumw2();
640  fhResidualT->Fit("gaus", "Q");
641  fhResidualT->GetFunction("gaus")->SetLineWidth(3);
642  fhResidualT->GetFunction("gaus")->SetLineColor(kRed);
644  fhResidualT->Write();
645  // if (fPrintToFileOn) gPad->Print(".gif");
646  // if (fPrintToFileOn) gPad->Print(".eps");
647  // cR5->cd();
648 
649  // TCanvas* c_alone = new TCanvas("c_alone","Pulls",400,400);
650  // new TCanvas("c_alone","Pulls",400,400);
651  // fhPullX->Draw();
652  // gPad->Print(".gif");
653  //
654  // TCanvas* c4x = new TCanvas("c4x","X-pulls vs pad size and cluster size",fnPadSizesX*300,3*300);
655  // c4x->Divide(fnPadSizesX,3);
656  // for (Int_t i=0;i<fnPadSizesX;i++){
657  // c4x->cd(i+1);
658  // fhPullXpads1[i]->Sumw2();
659  // fhPullXpads1[i]->Fit("gaus","Q");
660  // fhPullXpads1[i]->GetFunction("gaus")->SetLineWidth(2);
661  // fhPullXpads1[i]->GetFunction("gaus")->SetLineColor(kBlue);
662  // fhPullXpads1[i]->Draw();
663  // if (fPrintToFileOn) gPad->Print(".gif");
664  // if (fPrintToFileOn) gPad->Print(".eps");
665  //
666  // c4x->cd(fnPadSizesX+i+1);
667  // fhPullXpads2[i]->Sumw2();
668  // fhPullXpads2[i]->Fit("gaus","Q");
669  // fhPullXpads2[i]->GetFunction("gaus")->SetLineWidth(2);
670  // fhPullXpads2[i]->GetFunction("gaus")->SetLineColor(kBlue);
671  // fhPullXpads2[i]->Draw();
672  // if (fPrintToFileOn) gPad->Print(".gif");
673  // if (fPrintToFileOn) gPad->Print(".eps");
674  //
675  // c4x->cd(2*fnPadSizesX+i+1);
676  // fhPullXpads3[i]->Sumw2();
677  // fhPullXpads3[i]->Fit("gaus","Q");
678  // fhPullXpads3[i]->GetFunction("gaus")->SetLineWidth(2);
679  // fhPullXpads3[i]->GetFunction("gaus")->SetLineColor(kBlue);
680  // fhPullXpads3[i]->Draw();
681  // if (fPrintToFileOn) gPad->Print(".gif");
682  // if (fPrintToFileOn) gPad->Print(".eps");
683  // }
684 
685  // TCanvas* c4y = new TCanvas("c4y","Y-pulls vs pad size and cluster size",fnPadSizesY*300,3*300);
686  // c4y->Divide(fnPadSizesY,3);
687  // for (Int_t i=0;i<fnPadSizesY;i++){
688  // c4y->cd(i+1);
689  // fhPullYpads1[i]->Sumw2();
690  // fhPullYpads1[i]->Fit("gaus","Q");
691  // fhPullYpads1[i]->GetFunction("gaus")->SetLineWidth(2);
692  // fhPullYpads1[i]->GetFunction("gaus")->SetLineColor(kBlue);
693  // fhPullYpads1[i]->Draw();
694  // if (fPrintToFileOn) gPad->Print(".gif");
695  // if (fPrintToFileOn) gPad->Print(".eps");
696  //
697  // c4y->cd(fnPadSizesY+i+1);
698  // fhPullYpads2[i]->Sumw2();
699  // fhPullYpads2[i]->Fit("gaus","Q");
700  // fhPullYpads2[i]->GetFunction("gaus")->SetLineWidth(2);
701  // fhPullYpads2[i]->GetFunction("gaus")->SetLineColor(kBlue);
702  // fhPullYpads2[i]->Draw();
703  // if (fPrintToFileOn) gPad->Print(".gif");
704  // if (fPrintToFileOn) gPad->Print(".eps");
705  //
706  // c4y->cd(2*fnPadSizesY+i+1);
707  // fhPullYpads3[i]->Sumw2();
708  // fhPullYpads3[i]->Fit("gaus","Q");
709  // fhPullYpads3[i]->GetFunction("gaus")->SetLineWidth(2);
710  // fhPullYpads3[i]->GetFunction("gaus")->SetLineColor(kBlue);
711  // fhPullYpads3[i]->Draw();
712  // if (fPrintToFileOn) gPad->Print(".gif");
713  // if (fPrintToFileOn) gPad->Print(".eps");
714  // }
715  }
716 
717 
718  if (fOccupancyQaOn && fVerbose > 1) {
719  printf("===================================\n");
720  printf("OccupancyQa:\n");
721  // TCanvas* c1 = new TCanvas("c1","All pad distributions",1200,800);
722  // c1->Divide(NuOfColoums,NuOfRows);
723  for (Int_t i = 0; i < fNstations; i++) {
724  // c1->cd(i+1);
725  // fhPadsTotalR[i]->SetStats(0);
726  // fhPadsTotalR[i]->Draw("hist");
727  // if (fPrintToFileOn) gPad->Print(".gif");
728  // if (fPrintToFileOn) gPad->Print(".eps");
729  }
730  // c1->cd();
731 
732  // TCanvas* c2 = new TCanvas("c2","Fired pad distributions",1200,800);
733  // c2->Divide(NuOfColoums,NuOfRows);
734  for (Int_t i = 0; i < fNstations; i++) {
735  // c2->cd(i+1);
736  // fhPadsFiredR[i]->SetStats(0);
737  // fhPadsFiredR[i]->Draw();
738  fhPadsFiredR[i]->Write();
739  // if (fPrintToFileOn) gPad->Print(".gif");
740  // if (fPrintToFileOn) gPad->Print(".eps");
741  }
742  // c2->cd();
743 
744  // TCanvas* c3 = new TCanvas("c3","Occupancy plots",2400,1600);
745  // c3->Divide(NuOfColoums,NuOfRows);
746  for (Int_t i = 0; i < fNstations; i++) {
747  // c3->cd(i+1);
748  // fhOccupancyR[i]->SetStats(0);
749  // fhOccupancyR[i]->Draw();
750  fhOccupancyR[i]->Write();
751  // if (fPrintToFileOn) gPad->Print(".gif");
752  // if (fPrintToFileOn) gPad->Print(".eps");
753  }
754  // c3->cd();
755 
756  // TCanvas *c4 = new TCanvas("c4","Charge plot",2400,1600); ////added for station charge
757  // c4->Divide(3,2);
758  for (Int_t i = 0; i < fNstations; i++) {
759  // c4->cd(i+1);
760  // fhCharge[i]->Draw();
761  fhCharge[i]->Write();
762  // fhPointsInCluster[i]->Draw();
763  fhPointsInCluster[i]->Write();
764  fhDigisInCluster[i]->Write();
765  fhHitsInCluster[i]->Write();
766  }
767 
768 
769  /*
770  // TCanvas* oc_1 = new TCanvas("oc_1","",1200,1000);
771  new TCanvas("oc_1","",1200,1000);
772  fhOccupancyR[0]->GetYaxis()->SetTitleOffset(1.6);
773  fhOccupancyR[0]->Draw();
774  if (fPrintToFileOn) gPad->Print(".gif");
775  // TCanvas* oc_2 = new TCanvas("oc_2","",1200,1000);
776  new TCanvas("oc_2","",1200,1000);
777  fhOccupancyR[1]->GetYaxis()->SetTitleOffset(1.6);
778  fhOccupancyR[1]->Draw();
779  if (fPrintToFileOn) gPad->Print(".gif");
780  // TCanvas* oc_3 = new TCanvas("oc_3","",1200,1000);
781  new TCanvas("oc_3","",1200,1000);
782  fhOccupancyR[2]->GetYaxis()->SetTitleOffset(1.6);
783  fhOccupancyR[2]->Draw();
784  if (fPrintToFileOn) gPad->Print(".gif");*/
785  }
786 
787  if (fDigitizerQaOn && fVerbose > 1) {
788  printf("===================================\n");
789  printf("DigitizerQa:\n");
790 
791  TF1* fit_el = new TF1("fit_el", LandauMPV, -0.5, 4.5, 1);
792  fit_el->SetParameter(0, 0.511);
793  fit_el->SetLineWidth(2);
794  fit_el->SetLineColor(kBlack);
795 
796  TF1* fit_pi = new TF1("fit_pi", LandauMPV, -0.5, 4.5, 1);
797  fit_pi->SetParameter(0, 139.57);
798  fit_pi->SetLineWidth(2);
799  fit_pi->SetLineColor(kBlack);
800 
801  TF1* fit_pr = new TF1("fit_pr", LandauMPV, -0.5, 4.5, 1);
802  fit_pr->SetParameter(0, 938.27);
803  fit_pr->SetLineWidth(2);
804  fit_pr->SetLineColor(kBlack);
805 
806  TH1D* hLength = fhChargeTrackLength->ProjectionX();
807  TH1D* hLengthEl = fhChargeTrackLengthEl->ProjectionX();
808  ;
809  TH1D* hLengthPi = fhChargeTrackLengthPi->ProjectionX();
810  ;
811  TH1D* hLengthPr = fhChargeTrackLengthPr->ProjectionX();
812  ;
813  hLength->SetTitle("Track length for all tracks");
814  hLengthEl->SetTitle("Track length for electrons");
815  hLengthPi->SetTitle("Track length for pions");
816  hLengthPr->SetTitle("Track length for protons");
817 
818  fChargeHistos->Add(hLength);
819  fChargeHistos->Add(hLengthEl);
820  fChargeHistos->Add(hLengthPi);
821  fChargeHistos->Add(hLengthPr);
822 
823  // TCanvas* c5 = new TCanvas("c5","Charge vs energy and length",2100,1300);
824  // c5->Divide(4,3);
825  for (Int_t i = 0; i < 8; i++) {
826  // c5->cd(i+1);
827  // gPad->Range(0,0,200,200);
828  // gPad->SetBottomMargin(0.11);
829  // gPad->SetRightMargin(0.14);
830  //gPad->SetLeftMargin(0.12);
831  //gPad->SetLogz();
832  // TH2D* histo = (TH2D*) fChargeHistos->At(i);
833  // histo->Draw("colz");
834  // if (i==1) fit_el->Draw("same");
835  // if (i==2) fit_pi->Draw("same");
836  // if (i==3) fit_pr->Draw("same");
837  // if (fPrintToFileOn) gPad->Print(".gif");
838  // if (fPrintToFileOn) gPad->Print(".eps");
839  }
840 
841  // TCanvas* ch_1 = new TCanvas("ch_1","",1200,1000);
842  // new TCanvas("ch_1","",1200,1000);
843  // gPad->SetBottomMargin(0.11);
844  // gPad->SetRightMargin(0.14);
845  // gPad->SetLeftMargin(0.12);
846  // gPad->SetLogz();
847  // TH2D* histo1 = (TH2D*) fChargeHistos->At(1);
848  // histo1->Draw("colz");
849  // fit_el->Draw("same");
850  // if (fPrintToFileOn) gPad->Print(".gif");
851 
852  // TCanvas* ch_2 = new TCanvas("ch_2","",1200,1000);
853  /* new TCanvas("ch_2","",1200,1000);
854  gPad->SetBottomMargin(0.11);
855  gPad->SetRightMargin(0.14);
856  gPad->SetLeftMargin(0.12);
857  gPad->SetLogz();*/
858  // TH2D* histo2 = (TH2D*) fChargeHistos->At(2);
859  /* histo2->Draw("colz");
860  fit_pi->Draw("same");
861  if (fPrintToFileOn) gPad->Print(".gif");*/
862 
863  // TCanvas* ch_3 = new TCanvas("ch_3","",1200,1000);
864  /* new TCanvas("ch_3","",1200,1000);
865  gPad->SetBottomMargin(0.11);
866  gPad->SetRightMargin(0.14);
867  gPad->SetLeftMargin(0.12);
868  gPad->SetLogz(); */
869  // TH2D* histo3 = (TH2D*) fChargeHistos->At(3);
870  /* histo3->Draw("colz");
871  fit_pr->Draw("same");
872  if (fPrintToFileOn) gPad->Print(".gif"); */
873 
874 
875  for (Int_t i = 10; i < 14; i++) {
876  // c5->cd(i-1);
877  // gPad->SetLogy();
878  // gStyle->SetOptStat(1110);
879  // TH1D* histo = (TH1D*) fChargeHistos->At(i);
880  // histo->Draw();
881  //if (fPrintToFileOn) gPad->Print(".gif");
882  //if (fPrintToFileOn) gPad->Print(".eps");
883  }
884 
885  // below for 31st CBM Collaboration Meeting
886  // TCanvas* c6 = new TCanvas("c6","Charge distribution",1200,400);
887  // c6->Divide(2,1);
888  // c6->cd(1);
889  // fhCharge->Draw();
890  // fhCharge->Write();
891  // if (fPrintToFileOn) gPad->Print(".gif");
892  // if (fPrintToFileOn) gPad->Print(".eps");
893 
894  //c6->cd(2);
895  // fhChargeLog->Draw();
896  fhChargeLog->Write();
897  // if (fPrintToFileOn) gPad->Print(".gif");
898  // if (fPrintToFileOn) gPad->Print(".eps");
899 
900 
901  /* Commented below for 31st CBM Collaboration Meeting
902  TCanvas* c6 = new TCanvas("c6","Charge distribution",1200,400);
903  c6->Divide(3,1);
904  c6->cd(1);
905  fhCharge->Draw();
906  if (fPrintToFileOn) gPad->Print(".gif");
907  if (fPrintToFileOn) gPad->Print(".eps");
908 
909  c6->cd(2);
910  fhChargeLog->Draw();
911  if (fPrintToFileOn) gPad->Print(".gif");
912  if (fPrintToFileOn) gPad->Print(".eps");
913 
914  c6->cd(3);
915  fhChargePr_1GeV_3mm->Draw();
916  if (fPrintToFileOn) gPad->Print(".gif");
917  if (fPrintToFileOn) gPad->Print(".eps");
918  */
919 
920  // TCanvas* c8 = new TCanvas("c8","Square vs nPads",800,400);
921  // c8->Divide(2,1);
922  //c8->cd(1);
923  // fhNpadsVsS->Draw("colz");
924  fhNpadsVsS->Write();
925  Double_t nMean[11];
926  // Double_t s[11];
927  for (Int_t iS = 1; iS <= 10; iS++) {
928  nMean[iS] = 0;
929  // s[iS]=-5.25+0.5*iS;
930  Double_t total = 0;
931  for (Int_t iN = 1; iN <= 10; iN++) {
932  nMean[iS] += iN * fhNpadsVsS->GetBinContent(iS, iN);
933  total += fhNpadsVsS->GetBinContent(iS, iN);
934  }
935  if (total > 0) nMean[iS] /= total;
936  //printf("%f %f\n",s[iS],nMean[iS]);
937  }
938  // c8->cd(2);
939  // TGraph* gNvsS = new TGraph(11,s,nMean);
940  //gNvsS->DrawClone();
941 
942 
943  printf("All tracks: ;");
944  for (Int_t i = 0; i < fNstations; i++)
945  printf("%8i;", fNall[i]);
946  printf("\n");
947  printf("------------;");
948  for (Int_t i = 0; i < fNstations; i++)
949  printf("---------");
950  printf("\n");
951  printf("Primary: ;");
952  for (Int_t i = 0; i < fNstations; i++)
953  printf("%8i;", fNprimary[i]);
954  printf("\n");
955  printf("Secondary: ;");
956  for (Int_t i = 0; i < fNstations; i++)
957  printf("%8i;", fNsecondary[i]);
958  printf("\n");
959  printf("-------------");
960  for (Int_t i = 0; i < fNstations; i++)
961  printf("---------");
962  printf("\n");
963  printf("Protons: ;");
964  for (Int_t i = 0; i < fNstations; i++)
965  printf("%8i;", fNpr[i]);
966  printf("\n");
967  printf("Pions: ;");
968  for (Int_t i = 0; i < fNstations; i++)
969  printf("%8i;", fNpi[i]);
970  printf("\n");
971  printf("Electrons: ;");
972  for (Int_t i = 0; i < fNstations; i++)
973  printf("%8i;", fNel[i]);
974  printf("\n");
975  printf("Muons: ;");
976  for (Int_t i = 0; i < fNstations; i++)
977  printf("%8i;", fNmu[i]);
978  printf("\n");
979  printf("Kaons: ;");
980  for (Int_t i = 0; i < fNstations; i++)
981  printf("%8i;", fNka[i]);
982  printf("\n");
983  }
984 
985 
986  if (fStatisticsQaOn) {
987  printf("===================================\n");
988  printf("StatisticsQa:\n");
989  // TCanvas* c7 = new TCanvas("c7","Cluster statistics",1200,400);
990  // c7->Divide(3,1);
991  // c7->cd(1);
992  // gStyle->SetOptStat(1110);
993  // gPad->SetLogy();
994  // fhPointsInCluster->Draw();
995  // fhPointsInCluster->Write();
996  // if (fPrintToFileOn) gPad->Print(".gif");
997  // if (fPrintToFileOn) gPad->Print(".eps");
998  // c7->cd(2);
999  // gStyle->SetOptStat(1110);
1000  // gPad->SetLogy();
1001  // fhDigisInCluster->Draw();
1002  // fhDigisInCluster->Write();
1003  // if (fPrintToFileOn) gPad->Print(".gif");
1004  // if (fPrintToFileOn) gPad->Print(".eps");
1005  // c7->cd(3);
1006  // gStyle->SetOptStat(1110);
1007  // gPad->SetLogy();
1008  // fhHitsInCluster->Draw();
1009  // fhHitsInCluster->Write();
1010  // if (fPrintToFileOn) gPad->Print(".gif");
1011  // if (fPrintToFileOn) gPad->Print(".eps");
1012 
1013  printf("Total number of points: %i\n", fPointsTotal);
1014  printf("Points overcounted: %i\n", fPointsOverCounted);
1015  printf("Points undercounted: %i\n", fPointsUnderCounted);
1016  }
1017 
1018  if (fClusterDeconvQaOn) {
1019  printf("Signal points: %i\n", fSignalPoints);
1020  printf("Signal hits: %i\n", fSignalHits);
1021  }
1022 
1023  // TFile* performanceFile = new TFile(fFileName, "recreate");
1024  //
1025  // for (Int_t i=0;i<fNstations;i++) {
1026  // fhPadsTotalR[i]->Write();
1027  // fhPadsFiredR[i]->Write();
1028  // fhOccupancyR[i]->Write();
1029  // }
1030  //
1031  // fhPullX->Write();
1032  // fhPullY->Write();
1033  //
1034  // fhPointsInCluster->Write();
1035  // fhDigisInCluster->Write();
1036  // fhHitsInCluster->Write();
1037  //
1038  // for (Int_t i=0;i<10;i++) fChargeHistos->At(i)->Write();
1039  //
1040  // fhNpadsVsS->Write();
1041  // performanceFile->Close();
1042 }
1043 // -------------------------------------------------------------------------
1044 
1045 
1046 // -------------------------------------------------------------------------
1047 Double_t LandauMPV(Double_t* lg_x, Double_t* par) {
1048  Double_t gaz_gain_mean = 1.7e+4;
1049  Double_t scale = 1.e+6;
1050  gaz_gain_mean /= scale;
1051  Double_t mass = par[0]; // mass in MeV
1052  Double_t x = TMath::Power(10, lg_x[0]);
1053  return gaz_gain_mean * MPV_n_e(x, mass);
1054 }
1055 // -------------------------------------------------------------------------
1056 
1057 
1058 // -------------------------------------------------------------------------
1060  Bool_t verbose = (fVerbose > 2);
1061  TVector3 vIn; // in position of the track
1062  TVector3 vOut; // out position of the track
1063  TVector3 p; // track momentum
1064 
1065  fPointInfos->Clear();
1066 
1067  for (Int_t i = 0; i < fPoints->GetEntriesFast(); i++) {
1068  CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(i);
1069  Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
1070 
1071  // Check if the point corresponds to a certain MC Track
1072  Int_t trackId = point->GetTrackID();
1073  if (trackId < 0) {
1074  new ((*fPointInfos)[i]) CbmMuchPointInfo(0, 0, 0, 0, 0);
1075  continue;
1076  }
1077  CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(trackId);
1078  if (!mcTrack) {
1079  new ((*fPointInfos)[i]) CbmMuchPointInfo(0, 0, 0, 0, 0);
1080  continue;
1081  }
1082 
1083  Int_t motherId = mcTrack->GetMotherId();
1084 
1085  // Get mass
1086  Int_t pdgCode = mcTrack->GetPdgCode();
1087  TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
1088  if (!particle || pdgCode == 22 || // photons
1089  pdgCode == 2112) // neutrons
1090  {
1091  new ((*fPointInfos)[i]) CbmMuchPointInfo(0, 0, 0, 0, 0);
1092  continue;
1093  }
1094 
1095  if (CbmMuchAddress::GetLayerIndex(point->GetDetectorID()) == 0) {
1096  fNall[stId]++;
1097 
1098  if (pdgCode == 2212)
1099  fNpr[stId]++;
1100  else if (pdgCode == -211 || pdgCode == 211)
1101  fNpi[stId]++;
1102  else if (pdgCode == -11 || pdgCode == 11)
1103  fNel[stId]++;
1104  else if (pdgCode == -13 || pdgCode == 13)
1105  fNmu[stId]++;
1106  else if (pdgCode == -321 || pdgCode == 321)
1107  fNka[stId]++;
1108 
1109  if (motherId == -1)
1110  fNprimary[stId]++;
1111  else
1112  fNsecondary[stId]++;
1113  }
1114 
1115  Double_t mass = particle->Mass();
1116 
1117  point->PositionIn(vIn);
1118  point->PositionOut(vOut);
1119  point->Momentum(p);
1120  Double_t length = (vOut - vIn).Mag(); // Track length
1121  Double_t kine = sqrt(p.Mag2() + mass * mass) - mass; // Kinetic energy
1122  // Get mother pdg code
1123  Int_t motherPdg = 0;
1124  CbmMCTrack* motherTrack = NULL;
1125  if (motherId != -1) motherTrack = (CbmMCTrack*) fMCTracks->At(motherId);
1126  if (motherTrack) motherPdg = motherTrack->GetPdgCode();
1127  new ((*fPointInfos)[i])
1128  CbmMuchPointInfo(pdgCode, motherPdg, kine, length, stId);
1129  }
1130 
1131 
1132  // Filling generated charge for each point
1133  // Changing below loop from DigiMatches to Digis
1134 
1135  //for (Int_t i=0;i<fDigiMatches->GetEntriesFast();i++){
1136  // old code
1137  /*for (Int_t i=0;i<fDigis->GetEntriesFast();i++){
1138  // Get pad area
1139  CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(i); */
1140  for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
1142  assert(digi);
1143  //CbmMatch* match = (CbmMatch*) fDigiMatches->At(i);
1144  CbmMatch* match =
1146  assert(match);
1148  assert(module);
1149  if (!module) continue;
1150  Double_t area = 0;
1151  if (module->GetDetectorType() != 4 && module->GetDetectorType() != 3)
1152  continue;
1153  LOG(debug) << GetName() << " Processing MuchDigi " << i << " at address "
1154  << digi->GetAddress() << " Module number "
1155  << module->GetDetectorType();
1156 
1157  CbmMuchModuleGem* module1 = (CbmMuchModuleGem*) module;
1158  assert(module1);
1159  CbmMuchPad* pad = module1->GetPad(digi->GetAddress());
1160  assert(pad);
1161  area = pad->GetDx() * pad->GetDy();
1162  Int_t nofLinks = match->GetNofLinks();
1163  for (Int_t pt = 0; pt < nofLinks; pt++) {
1164  Int_t pointId = match->GetLink(pt).GetIndex();
1165  Int_t charge = match->GetLink(pt).GetWeight();
1166  CbmMuchPointInfo* info = (CbmMuchPointInfo*) fPointInfos->At(pointId);
1167  if (info->GetPdgCode() == 0) continue;
1168  info->AddCharge(charge);
1169  info->AddArea(area);
1170  }
1171  }
1172 
1173  //Filling digitizer performance plots
1174  for (Int_t i = 0; i < fPointInfos->GetEntriesFast(); i++) {
1176  if (verbose) {
1177  printf("%i", i);
1178  info->Print();
1179  }
1180  Double_t length = info->GetLength();
1181  Double_t kine = 1000 * (info->GetKine());
1182  Double_t charge = info->GetCharge();
1183  Int_t pdg = info->GetPdgCode();
1184  if (pdg == 0) continue;
1185  if (charge <= 0) continue;
1186  Double_t log_kine = TMath::Log10(kine);
1187  Double_t log_charge = TMath::Log10(charge);
1188  charge = charge / 1.e+4; // measure charge in 10^4 electrons
1189 
1190  Int_t nPads = info->GetNPads();
1191  Double_t area = info->GetArea() / nPads;
1192  //printf("%f %i\n",TMath::Log2(area),nPads);
1193  fhNpadsVsS->Fill(TMath::Log2(area), nPads);
1194  fhCharge[info->GetStationId()]->Fill(charge);
1195  fhChargeLog->Fill(log_charge);
1196  fhChargeEnergyLog->Fill(log_kine, charge);
1197  fhChargeTrackLength->Fill(length, charge);
1198 
1199  if (pdg == 2212)
1200  fhChargeEnergyLogPr->Fill(log_kine, charge);
1201  else if (pdg == 211 || pdg == -211)
1202  fhChargeEnergyLogPi->Fill(log_kine, charge);
1203  else if (pdg == 11 || pdg == -11)
1204  fhChargeEnergyLogEl->Fill(log_kine, charge);
1205 
1206  if (pdg == 2212)
1207  fhChargeTrackLengthPr->Fill(length, charge);
1208  else if (pdg == 211 || pdg == -211)
1209  fhChargeTrackLengthPi->Fill(length, charge);
1210  else if (pdg == 11 || pdg == -11)
1211  fhChargeTrackLengthEl->Fill(length, charge);
1212 
1213  if (pdg == 2212 && length > 0.3 && length < 0.32 && kine > 1000
1214  && kine < 1020)
1215  fhChargePr_1GeV_3mm->Fill(charge);
1216  }
1217 }
1218 // -------------------------------------------------------------------------
1219 
1220 
1221 // -------------------------------------------------------------------------
1223  // Bool_t verbose = (fVerbose>2);
1224  // Filling occupancy plots
1225  // old code
1226  /*for (Int_t i=0;i<fDigis->GetEntriesFast();i++){
1227  CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(i); */
1228  for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
1230  assert(digi);
1231  UInt_t address = digi->GetAddress();
1232  CbmMuchModule* module = fGeoScheme->GetModuleByDetId(address);
1233  if (!module) continue;
1234  Double_t r0 = 0;
1235  if (module->GetDetectorType() != 4 && module->GetDetectorType() != 3)
1236  continue;
1237  CbmMuchModuleGem* module1 = (CbmMuchModuleGem*) module;
1238  CbmMuchPad* pad = module1->GetPad(
1239  address); //fGeoScheme->GetPadByDetId(detectorId, channelId);
1240  assert(pad);
1241  Double_t x0 = pad->GetX();
1242  Double_t y0 = pad->GetY();
1243  r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
1244  fhPadsFiredR[CbmMuchAddress::GetStationIndex(address)]->Fill(r0);
1245  }
1246 }
1247 // -------------------------------------------------------------------------
1248 
1249 
1250 // -------------------------------------------------------------------------
1252  // Bool_t verbose = (fVerbose>2);
1253  Int_t nClusters = fClusters->GetEntriesFast();
1254  TArrayI hitsInCluster;
1255  TArrayI pointsInCluster;
1256  hitsInCluster.Set(nClusters);
1257  pointsInCluster.Set(nClusters);
1258  cout << " start Stat QA " << endl;
1259  for (Int_t i = 0; i < nClusters; i++)
1260  hitsInCluster[i] = 0;
1261  for (Int_t i = 0; i < nClusters; i++)
1262  pointsInCluster[i] = 0;
1263 
1264  for (Int_t i = 0; i < fHits->GetEntriesFast(); i++) {
1265  CbmMuchPixelHit* hit = (CbmMuchPixelHit*) fHits->At(i);
1266  //cout<<" hit index "<<i<<" plane id "<<hit->GetPlaneId()<<" x "<<hit->GetX()<<" y "<<hit->GetY()<<" z "<<hit->GetZ()<<" cluster Id "<< hit->GetRefId()<<endl;
1267 
1268  // cout<<" hit index "<<i<<" plane id "<<hit->GetPlaneId()<<endl;
1269  Int_t clusterId = hit->GetRefId();
1270  hitsInCluster[clusterId]++;
1271  }
1272 
1273  for (Int_t i = 0; i < nClusters; i++) {
1274  CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(i);
1275 
1276  map<Int_t, Int_t> map_points;
1277  Int_t nDigis = cluster->GetNofDigis();
1278 
1279  auto address = cluster->GetAddress();
1280  auto StationIndex = CbmMuchAddress::GetStationIndex(address);
1281  //cout<<" station index "<<StationIndex<<endl;
1282 
1283  fhDigisInCluster[StationIndex]->Fill(nDigis);
1284 
1285  for (Int_t digiId = 0; digiId < nDigis; digiId++) {
1286  Int_t index = cluster->GetDigi(digiId);
1287  //Access Match from CbmDigi only
1288  // CbmMuchDigi* digi= (CbmMuchDigi*) fDigis->At(index);
1289  //CbmMatch* match = (CbmMatch*) fDigiMatches->At(index);
1290  CbmMatch* match =
1292  Int_t nPoints = match->GetNofLinks();
1293  for (Int_t iRefPoint = 0; iRefPoint < nPoints; iRefPoint++) {
1294  Int_t pointId = match->GetLink(iRefPoint).GetIndex();
1295  map_points[pointId] = 1;
1296  }
1297  }
1298  pointsInCluster[i] = map_points.size();
1299  map_points.clear();
1300  }
1301 
1302 
1303  for (Int_t i = 0; i < nClusters; i++) {
1304  // added
1305  CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(i);
1306  auto address = cluster->GetAddress();
1307  auto StationIndex = CbmMuchAddress::GetStationIndex(address);
1308  //cout<<" station index "<<StationIndex<<endl;
1310 
1311  Int_t nPts = pointsInCluster[i];
1312  Int_t nHits = hitsInCluster[i];
1313  fhPointsInCluster[StationIndex]->Fill(nPts);
1314  // fhPointsInCluster->Fill(nPts);
1315  fhHitsInCluster[StationIndex]->Fill(nHits);
1316  if (nPts > nHits) fPointsUnderCounted += (nPts - nHits);
1317  if (nHits > nPts) fPointsOverCounted += (nHits - nPts);
1318  fPointsTotal += nPts;
1319  }
1320 }
1321 // -------------------------------------------------------------------------
1322 
1323 
1324 // -------------------------------------------------------------------------
1326  Bool_t verbose = (fVerbose > 2);
1327  // Filling residuals and pools for hits at the first layer
1328  for (Int_t i = 0; i < fHits->GetEntriesFast(); i++) {
1329  CbmMuchPixelHit* hit = (CbmMuchPixelHit*) fHits->At(i);
1330  // Select hits from the first station only
1331  Int_t iStation = CbmMuchAddress::GetStationIndex(hit->GetAddress());
1332  Int_t iLayer = CbmMuchAddress::GetLayerIndex(hit->GetAddress());
1333  // if((iStation !=0 && iStation !=1))continue;
1334  // if((iStation !=0 && iStation !=1))continue;
1335  // cout<<" PULLS QA STATION INDEX "<<iStation<<endl;
1336  //Earlier finding for only one station
1337  //if(!(iStation == 0)) continue;
1338  // if(!(iStation == 3 && iLayer == 0)) continue;
1339  if (verbose)
1340  printf(" Hit %i, station %i, layer %i ", i, iStation, iLayer);
1341 
1342  Int_t clusterId = hit->GetRefId();
1343  Bool_t hit_unique = 1;
1344  for (Int_t j = i + 1; j < fHits->GetEntriesFast(); j++) {
1345  CbmMuchPixelHit* hit1 = (CbmMuchPixelHit*) fHits->At(j);
1346  //if (hit1->GetStationNr()>stationNr) break;
1347  if (hit1->GetRefId() == clusterId) {
1348  hit_unique = 0;
1349  break;
1350  }
1351  }
1352  if (verbose) printf("hit_unique=%i", hit_unique);
1353  if (!hit_unique) {
1354  if (verbose) printf("\n");
1355  continue;
1356  }
1357 
1358  // Select hits with clusters formed by only one MC Point
1359  CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(clusterId);
1360  Bool_t point_unique = 1;
1361  Int_t pointId = -1;
1362 
1363  /* Double_t xmax = 0;
1364  Double_t xmin = 0;
1365  Double_t ymax = 0;
1366  Double_t ymin = 0;
1367  Double_t dxmin = 0;
1368  Double_t dymin = 0;
1369  */
1370  // if (cluster->GetNDigis()>1) {if (verbose) printf("\n"); continue;}
1371  // cout<<" digi number per cluster "<<cluster->GetNofDigis()<<endl;
1372  for (Int_t digiId = 0; digiId < cluster->GetNofDigis(); digiId++) {
1373  Int_t index = cluster->GetDigi(digiId);
1374  // printf("%i\n",index);
1375  // CbmMuchDigi* digi= (CbmMuchDigi*) fDigis->At(index);
1376  CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(index);
1377  // cout<<" check 1"<<endl;
1378  if (!digi) continue;
1379  if (index < 0) continue;
1380  //CbmMatch* match = (CbmMatch*) fDigiMatches->At(index);
1381  CbmMatch* match =
1383  if (!match) continue;
1384  // cout<<" check 2 "<<endl;
1385  //CbmMuchDigiMatch* match = (CbmMuchDigiMatch*) fDigiMatches->At(index);
1386  // Not unique if the pad has several mcPoint references
1387  if (verbose) printf(" n=%i", match->GetNofLinks());
1388  if (match->GetNofLinks() == 0) {
1389  printf(" noise hit");
1390  point_unique = 0;
1391  break;
1392  }
1393  if (match->GetNofLinks() > 1) {
1394  point_unique = 0;
1395  break;
1396  }
1397  Int_t currentPointId = match->GetLink(0).GetIndex();
1398  CbmMuchModuleGem* module =
1400  if (!module) continue;
1401  // CbmMuchPad* pad = module->GetPad(digi->GetAddress());//fGeoScheme->GetPadByDetId(digi->GetDetectorId(), digi->GetChannelId());
1402  /*Double_t x = pad->GetX();
1403  Double_t y = pad->GetY();
1404  Double_t dx = pad->GetDx();
1405  Double_t dy = pad->GetDy();
1406  if (digiId==0 || dxmin>dx ) dxmin=dx;
1407  if (digiId==0 || dymin>dy ) dymin=dy;
1408  if (digiId==0 || xmin>x-dx/2) xmin=x-dx/2;
1409  if (digiId==0 || xmax<x+dx/2) xmax=x+dx/2;
1410  if (digiId==0 || ymin>y-dy/2) ymin=y-dy/2;
1411  if (digiId==0 || ymax<y+dy/2) ymax=y+dy/2;*/
1412  if (digiId == 0) {
1413  pointId = currentPointId;
1414  continue;
1415  }
1416  // Not unique if mcPoint references differ for different digis
1417  if (currentPointId != pointId) {
1418  point_unique = 0;
1419  break;
1420  }
1421  }
1422 
1423 
1424  if (verbose) printf(" point_unique=%i", point_unique);
1425  if (!point_unique) {
1426  if (verbose) printf("\n");
1427  continue;
1428  }
1429  //printf(" %f %f %f %f %f %f\n",xmin,xmax,ymin,ymax,dxmin,dymin);
1430  // Int_t nPadsX=Int_t((xmax-xmin)/dxmin);
1431  // Int_t nPadsY=Int_t((ymax-ymin)/dymin);
1432  //printf("nPadsX=%i nPadsY=%i\n",nPadsX,nPadsY);
1433 
1434  if (verbose) printf(" pointId=%i", pointId);
1435  CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(pointId);
1436 
1437  Double_t xMC = 0.5 * (point->GetXIn() + point->GetXOut());
1438  Double_t yMC = 0.5 * (point->GetYIn() + point->GetYOut());
1439  Double_t tMC = point->GetTime();
1440  // cout<<" MC point time "<<tMC<<" z "<<point->GetZ()<<endl;
1441  Double_t xRC = hit->GetX();
1442  Double_t yRC = hit->GetY();
1443  Double_t dx = hit->GetDx();
1444  Double_t dy = hit->GetDy();
1445 
1446  Double_t tRC = hit->GetTime();
1447  Double_t dt = hit->GetTimeError();
1448  // cout<<" Rec Hit time "<<tRC<<endl;
1449 
1450  if (dx < 1.e-10) {
1451  printf("Anomalously small dx\n");
1452  continue;
1453  }
1454  if (dy < 1.e-10) {
1455  printf("Anomalously small dy\n");
1456  continue;
1457  }
1458  fhPullX->Fill((xRC - xMC) / dx);
1459  fhPullY->Fill((yRC - yMC) / dy);
1460  fhPullT->Fill((tRC - tMC) / dt);
1461  fhResidualX->Fill((xRC - xMC));
1462  fhResidualY->Fill((yRC - yMC));
1463  fhResidualT->Fill((tRC - tMC));
1464  // Int_t nDigis = cluster->GetNDigis();
1465  // if (nDigis<=fNTimingPulls) fhPullT[nDigis]->Fill((tRC-tMC)/dt);
1466 
1467  if (verbose) printf("\n");
1468 
1469  // Int_t index = cluster->GetDigiIndex(0);
1470  //
1471  // // printf("index=%i\n",index);
1472  // CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(index);
1473  //
1474  //
1475  // Int_t padSizeX = Int_t(TMath::Log2(dxmin/fPadMinLx));
1476  // Int_t padSizeY = Int_t(TMath::Log2(dymin/fPadMinLy));
1477  // if (padSizeX>=fnPadSizesX || padSizeX<0) { printf("wrong x pad size\n"); continue; }
1478  // if (padSizeY>=fnPadSizesY || padSizeY<0) { printf("wrong y pad size\n"); continue; }
1479  // if (nPadsX==1 && nPadsY==1) fhPullXpads1[padSizeX]->Fill((xRC-xMC)/dx);
1480  // if (nPadsY==1 && nPadsX==1) fhPullYpads1[padSizeY]->Fill((yRC-yMC)/dy);
1481  // if (nPadsX==2 && nPadsY==1) fhPullXpads2[padSizeX]->Fill((xRC-xMC)/dx);
1482  // if (nPadsY==2 && nPadsX==1) fhPullYpads2[padSizeY]->Fill((yRC-yMC)/dy);
1483  // if (nPadsX==3 && nPadsY==1) fhPullXpads3[padSizeX]->Fill((xRC-xMC)/dx);
1484  // if (nPadsY==3 && nPadsX==1) fhPullYpads3[padSizeY]->Fill((yRC-yMC)/dy);
1485  }
1486 }
1487 // -------------------------------------------------------------------------
1488 
1489 
1490 // -------------------------------------------------------------------------
1492  Int_t nPoints = fPoints->GetEntriesFast();
1493  // Int_t nMatches = fDigiMatches->GetEntriesFast();
1494  Int_t nClusters = fClusters->GetEntriesFast();
1495  vector<Int_t> pIndices;
1496  vector<Int_t>::iterator it;
1497 
1498  for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
1499  if (IsSignalPoint(iPoint)) fSignalPoints++;
1500  }
1501 
1502  for (Int_t iCluster = 0; iCluster < nClusters; ++iCluster) {
1503  CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(iCluster);
1504  if (!cluster) continue;
1505  Int_t nDigis = cluster->GetNofDigis();
1506  for (Int_t id = 0; id < nDigis; ++id) {
1507  Int_t iDigi = cluster->GetDigi(id);
1508  // CbmMuchDigi* digi= (CbmMuchDigi*) fDigis->At(iDigi);
1509  //CbmMatch* match = (CbmMatch*) fDigiMatches->At(iDigi);
1510  CbmMatch* match =
1512 
1513  //CbmMuchDigiMatch* match = (CbmMuchDigiMatch*) fDigiMatches->At(iDigi);
1514  if (!match) continue;
1515  for (Int_t ip = 0; ip < match->GetNofLinks(); ++ip) {
1516  Int_t iPoint = match->GetLink(ip).GetIndex();
1517  it = find(pIndices.begin(), pIndices.end(), iPoint);
1518  if (it != pIndices.end()) continue;
1519  pIndices.push_back(iPoint);
1520  if (IsSignalPoint(iPoint)) fSignalHits++;
1521  }
1522  }
1523  }
1524 }
1525 // -------------------------------------------------------------------------
1526 
1528  Int_t nPoints = fPoints->GetEntriesFast();
1529  Int_t nTracks = fMCTracks->GetEntriesFast();
1530  CbmMuchPoint* point = (iPoint < 0 || iPoint > nPoints - 1)
1531  ? NULL
1532  : (CbmMuchPoint*) fPoints->At(iPoint);
1533  if (!point) return kFALSE;
1534  Int_t iTrack = point->GetTrackID();
1535  CbmMCTrack* track = (iTrack < 0 || iTrack > nTracks - 1)
1536  ? NULL
1537  : (CbmMCTrack*) fMCTracks->At(iTrack);
1538  if (!track) return kFALSE;
1539  if (iTrack != 0 && iTrack != 1)
1540  return kFALSE; // Signal tracks must be fist ones
1541  // Verify if it is a signal muon
1542  if (track->GetMotherId() < 0) { // No mother for signal
1543  Int_t pdgCode = track->GetPdgCode();
1544  if (TMath::Abs(pdgCode) == 13) { // If it is a muon
1545  return kTRUE;
1546  }
1547  }
1548  return kFALSE;
1549 }
1550 
1551 // -------------------------------------------------------------------------
1552 Int_t CbmMuchHitFinderQa::GetNChannels(Int_t iStation) {
1553  Int_t nChannels = 0;
1554  vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
1555  for (UInt_t im = 0; im < modules.size(); im++) {
1556  CbmMuchModule* mod = modules[im];
1557  if (mod->GetDetectorType() != 4 && mod->GetDetectorType() != 3) continue;
1558  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
1559  if (!module) continue;
1560  nChannels += module->GetNPads();
1561  }
1562  return nChannels;
1563 }
1564 
1565 // -------------------------------------------------------------------------
1566 Int_t CbmMuchHitFinderQa::GetNSectors(Int_t iStation) {
1567  Int_t nSectors = 0;
1568  vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
1569  for (UInt_t im = 0; im < modules.size(); im++) {
1570  CbmMuchModule* mod = modules[im];
1571  if (mod->GetDetectorType() != 4 && mod->GetDetectorType() != 3) continue;
1572  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
1573  if (!module) continue;
1574  nSectors += module->GetNSectors();
1575  }
1576  return nSectors;
1577 }
1578 
1579 // -------------------------------------------------------------------------
1580 TVector2 CbmMuchHitFinderQa::GetMinPadSize(Int_t iStation) {
1581  Double_t padMinLx = std::numeric_limits<Double_t>::max();
1582  Double_t padMinLy = std::numeric_limits<Double_t>::max();
1583  vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
1584  for (UInt_t im = 0; im < modules.size(); im++) {
1585  CbmMuchModule* mod = modules[im];
1586  if (mod->GetDetectorType() != 1) continue;
1587  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
1588  vector<CbmMuchPad*> pads = module->GetPads();
1589  for (UInt_t ip = 0; ip < pads.size(); ip++) {
1590  CbmMuchPad* pad = pads[ip];
1591  if (pad->GetDx() < padMinLx) padMinLx = pad->GetDx();
1592  if (pad->GetDy() < padMinLy) padMinLy = pad->GetDy();
1593  }
1594  }
1595  return TVector2(padMinLx, padMinLy);
1596 }
1597 // -------------------------------------------------------------------------
1598 
1599 // -------------------------------------------------------------------------
1600 TVector2 CbmMuchHitFinderQa::GetMaxPadSize(Int_t iStation) {
1601  Double_t padMaxLx = std::numeric_limits<Double_t>::min();
1602  Double_t padMaxLy = std::numeric_limits<Double_t>::min();
1603  vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
1604  for (UInt_t im = 0; im < modules.size(); im++) {
1605  CbmMuchModule* mod = modules[im];
1606  if (mod->GetDetectorType() != 1) continue;
1607  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
1608  vector<CbmMuchPad*> pads = module->GetPads();
1609  for (UInt_t ip = 0; ip < pads.size(); ip++) {
1610  CbmMuchPad* pad = pads[ip];
1611  if (pad->GetDx() > padMaxLx) padMaxLx = pad->GetDx();
1612  if (pad->GetDy() > padMaxLy) padMaxLy = pad->GetDy();
1613  }
1614  }
1615  return TVector2(padMaxLx, padMaxLy);
1616 }
1617 // -------------------------------------------------------------------------
1618 
1619 // -------------------------------------------------------------------------
1620 Double_t MPV_n_e(Double_t Tkin, Double_t mass) {
1621  Double_t logT;
1622  TF1 fPol6("fPol6", "pol6", -5, 10);
1623  if (mass < 0.1) {
1624  logT = log(Tkin * 0.511 / mass);
1625  if (logT > 9.21034) logT = 9.21034;
1626  if (logT < min_logT_e) logT = min_logT_e;
1627  return fPol6.EvalPar(&logT, mpv_e);
1628  } else if (mass >= 0.1 && mass < 0.2) {
1629  logT = log(Tkin * 105.658 / mass);
1630  if (logT > 9.21034) logT = 9.21034;
1631  if (logT < min_logT_mu) logT = min_logT_mu;
1632  return fPol6.EvalPar(&logT, mpv_mu);
1633  } else {
1634  logT = log(Tkin * 938.272 / mass);
1635  if (logT > 9.21034) logT = 9.21034;
1636  if (logT < min_logT_p) logT = min_logT_p;
1637  return fPol6.EvalPar(&logT, mpv_p);
1638  }
1639 }
1640 // -------------------------------------------------------------------------
1641 
CbmMuchPoint::GetXOut
Double_t GetXOut() const
Definition: CbmMuchPoint.h:73
CbmMuchHitFinderQa::fNpi
Int_t * fNpi
Definition: CbmMuchHitFinderQa.h:137
CbmMuchHitFinderQa::GetNChannels
Int_t GetNChannels(Int_t iStation)
Definition: CbmMuchHitFinderQa.cxx:1552
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMuchHitFinderQa::fNall
Int_t * fNall
Definition: CbmMuchHitFinderQa.h:135
CbmMatch
Definition: CbmMatch.h:22
CbmMuchHitFinderQa::fPointsOverCounted
Int_t fPointsOverCounted
Definition: CbmMuchHitFinderQa.h:146
CbmMuchHitFinderQa::fStatisticsQaOn
Bool_t fStatisticsQaOn
Definition: CbmMuchHitFinderQa.h:151
CbmMuchHitFinderQa::FinishTask
virtual void FinishTask()
Definition: CbmMuchHitFinderQa.cxx:540
CbmMuchDigi.h
CbmMuchGeoScheme
Definition: CbmMuchGeoScheme.h:43
CbmMuchHitFinderQa::fNel
Int_t * fNel
Definition: CbmMuchHitFinderQa.h:138
CbmMuchPoint
Definition: CbmMuchPoint.h:21
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
CbmMuchHitFinderQa::fVerbose
Int_t fVerbose
Definition: CbmMuchHitFinderQa.h:78
CbmMuchHitFinderQa::fNprimary
Int_t * fNprimary
Definition: CbmMuchHitFinderQa.h:141
CbmMuchHitFinderQa::SetParContainers
virtual void SetParContainers()
Definition: CbmMuchHitFinderQa.cxx:470
CbmMuchStation
Definition: CbmMuchStation.h:22
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
CbmMuchPointInfo::GetNPads
Int_t GetNPads()
Definition: CbmMuchPointInfo.h:41
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmMuchHitFinderQa::fDigiManager
CbmDigiManager * fDigiManager
Definition: CbmMuchHitFinderQa.h:83
BINNING_CHARGE_LOG
#define BINNING_CHARGE_LOG
CbmMuchHitFinderQa::fhPullY
TH1D * fhPullY
Definition: CbmMuchHitFinderQa.h:123
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmMuchPointInfo::AddArea
void AddArea(Double_t s)
Definition: CbmMuchPointInfo.h:28
CbmMuchHitFinderQa::padsFile
FILE * padsFile
Definition: CbmMuchHitFinderQa.h:162
CbmMuchHitFinderQa::fNpr
Int_t * fNpr
Definition: CbmMuchHitFinderQa.h:136
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition: CbmDigiManager.cxx:71
CbmMuchModuleGem::GetPads
std::vector< CbmMuchPad * > GetPads()
Definition: CbmMuchModuleGem.cxx:41
BINNING_CHARGE
#define BINNING_CHARGE
CbmMuchCluster
Data container for MUCH clusters.
Definition: CbmMuchCluster.h:20
CbmMuchHitFinderQa::fPointInfos
TClonesArray * fPointInfos
Definition: CbmMuchHitFinderQa.h:88
CbmMuchHitFinderQa::fhChargeEnergyLog
TH2D * fhChargeEnergyLog
Definition: CbmMuchHitFinderQa.h:93
CbmMuchHitFinderQa::fNka
Int_t * fNka
Definition: CbmMuchHitFinderQa.h:140
CbmMuchHitFinderQa::fhResidualY
TH1D * fhResidualY
Definition: CbmMuchHitFinderQa.h:128
CbmMuchHitFinderQa::CbmMuchHitFinderQa
CbmMuchHitFinderQa(const char *name="MuchHitFinderQa", Int_t verbose=1)
Definition: CbmMuchHitFinderQa.cxx:55
CbmMuchHitFinderQa::IsSignalPoint
Bool_t IsSignalPoint(Int_t iPoint)
Definition: CbmMuchHitFinderQa.cxx:1527
CbmMuchPad::GetY
Double_t GetY() const
Definition: CbmMuchPad.h:29
CbmMuchHitFinderQa::fhChargeTrackLengthEl
TH2D * fhChargeTrackLengthEl
Definition: CbmMuchHitFinderQa.h:100
CbmMuchDigi::GetAdc
UShort_t GetAdc() const
Definition: CbmMuchDigi.h:74
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
mpv_p
constexpr double mpv_p[]
Definition: CbmMuchRecoDefs.h:27
CbmMuchHitFinderQa::fPointsUnderCounted
Int_t fPointsUnderCounted
Definition: CbmMuchHitFinderQa.h:145
CbmMuchHitFinderQa::fhChargePr_1GeV_3mm
TH1D * fhChargePr_1GeV_3mm
Definition: CbmMuchHitFinderQa.h:103
CbmMuchPad::GetDy
Double_t GetDy() const
Definition: CbmMuchPad.h:31
CbmMuchPointInfo::GetStationId
Int_t GetStationId()
Definition: CbmMuchPointInfo.h:39
CbmMuchHitFinderQa::fGeoScheme
CbmMuchGeoScheme * fGeoScheme
Definition: CbmMuchHitFinderQa.h:73
CbmDigiManager::GetNofDigis
static Int_t GetNofDigis(ECbmModuleId systemId)
Definition: CbmDigiManager.cxx:62
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
CbmMuchGeoScheme::Init
void Init(TObjArray *stations, Int_t flag)
Definition: CbmMuchGeoScheme.cxx:121
BINNING_ENERGY_LOG
#define BINNING_ENERGY_LOG
CbmMuchHitFinderQa::fhNpadsVsS
TH2D * fhNpadsVsS
Definition: CbmMuchHitFinderQa.h:104
LandauMPV
Double_t LandauMPV(Double_t *lg_x, Double_t *par)
Definition: CbmMuchHitFinderQa.cxx:1047
CbmMatch.h
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmHit::GetTimeError
Double_t GetTimeError() const
Definition: CbmHit.h:76
CbmMuchPad::GetDx
Double_t GetDx() const
Definition: CbmMuchPad.h:30
CbmMuchHitFinderQa::fhPointsInCluster
TH1I ** fhPointsInCluster
Definition: CbmMuchHitFinderQa.h:131
CbmMuchHitFinderQa.h
CbmMuchHitFinderQa::fOccupancyQaOn
Bool_t fOccupancyQaOn
Definition: CbmMuchHitFinderQa.h:148
CbmMuchHitFinderQa::fhChargeTrackLength
TH2D * fhChargeTrackLength
Definition: CbmMuchHitFinderQa.h:97
CbmMuchPointInfo
Definition: CbmMuchPointInfo.h:11
CbmMuchHitFinderQa::fSignalHits
Int_t fSignalHits
Definition: CbmMuchHitFinderQa.h:77
CbmMuchHitFinderQa::OccupancyQa
void OccupancyQa()
Definition: CbmMuchHitFinderQa.cxx:1222
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
BINNING_LENGTH
#define BINNING_LENGTH
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:33
CbmMuchPoint::GetYOut
Double_t GetYOut() const
Definition: CbmMuchPoint.h:74
CbmMuchPointInfo::GetCharge
Int_t GetCharge()
Definition: CbmMuchPointInfo.h:38
CbmMuchCluster.h
Data container for MUCH clusters.
CbmMuchHitFinderQa::fhChargeEnergyLogEl
TH2D * fhChargeEnergyLogEl
Definition: CbmMuchHitFinderQa.h:96
CbmMuchHitFinderQa::fPullsQaOn
Bool_t fPullsQaOn
Definition: CbmMuchHitFinderQa.h:149
min_logT_p
constexpr double min_logT_p
Definition: CbmMuchRecoDefs.h:31
CbmMuchHitFinderQa::fhChargeLog
TH1D * fhChargeLog
Definition: CbmMuchHitFinderQa.h:102
CbmMuchRecoDefs.h
CbmMuchHitFinderQa::fhPullT
TH1D * fhPullT
Definition: CbmMuchHitFinderQa.h:124
CbmMuchHitFinderQa::fhChargeTrackLengthPi
TH2D * fhChargeTrackLengthPi
Definition: CbmMuchHitFinderQa.h:98
CbmMuchHitFinderQa::fhDigisInCluster
TH1I ** fhDigisInCluster
Definition: CbmMuchHitFinderQa.h:132
CbmMuchGeoScheme::GetNStations
Int_t GetNStations() const
Definition: CbmMuchGeoScheme.h:79
CbmMuchHitFinderQa::pointsFile
FILE * pointsFile
Definition: CbmMuchHitFinderQa.h:161
BINNING_ENERGY_LOG_EL
#define BINNING_ENERGY_LOG_EL
CbmMuchHitFinderQa::fNmu
Int_t * fNmu
Definition: CbmMuchHitFinderQa.h:139
CbmMuchPoint.h
min_logT_mu
constexpr double min_logT_mu
Definition: CbmMuchRecoDefs.h:30
CbmGeoMuchPar.h
CbmMuchPointInfo::GetLength
Double_t GetLength()
Definition: CbmMuchPointInfo.h:35
CbmMuchModuleGem::GetNPads
Int_t GetNPads()
Definition: CbmMuchModuleGem.cxx:58
CbmMuchHitFinderQa
Definition: CbmMuchHitFinderQa.h:26
CbmMuchPoint::PositionIn
void PositionIn(TVector3 &pos) const
Definition: CbmMuchPoint.h:79
CbmMuchHitFinderQa::GetMinPadSize
TVector2 GetMinPadSize(Int_t iStation)
Definition: CbmMuchHitFinderQa.cxx:1580
CbmMuchHitFinderQa::fhChargeEnergyLogPr
TH2D * fhChargeEnergyLogPr
Definition: CbmMuchHitFinderQa.h:95
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmMuchStation::GetRmax
Double_t GetRmax() const
Definition: CbmMuchStation.h:52
CbmMuchHitFinderQa::fGeoFileName
TString fGeoFileName
Definition: CbmMuchHitFinderQa.h:74
CbmMuchHitFinderQa::Exec
virtual void Exec(Option_t *option)
Definition: CbmMuchHitFinderQa.cxx:481
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmMuchPad::GetAddress
Int_t GetAddress() const
Definition: CbmMuchPad.h:27
min_logT_e
constexpr double min_logT_e
Definition: CbmMuchRecoDefs.h:29
MPV_n_e
Double_t MPV_n_e(Double_t Tkin, Double_t mass)
Definition: CbmMuchHitFinderQa.cxx:1620
CbmMuchHitFinderQa::fPrintToFileOn
Bool_t fPrintToFileOn
Definition: CbmMuchHitFinderQa.h:154
CbmMuchPointInfo::Print
void Print(Option_t *="") const
Definition: CbmMuchPointInfo.cxx:44
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
CbmMuchModuleGem
Definition: CbmMuchModuleGem.h:24
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
CbmMuchHitFinderQa::fhPullX
TH1D * fhPullX
Definition: CbmMuchHitFinderQa.h:122
CbmMuchHitFinderQa::fhResidualT
TH1D * fhResidualT
Definition: CbmMuchHitFinderQa.h:129
CbmDigiManager::GetMatch
const CbmMatch * GetMatch(ECbmModuleId systemId, UInt_t index) const
Get a match object.
Definition: CbmDigiManager.cxx:54
CbmMuchHitFinderQa::fPoints
TClonesArray * fPoints
Definition: CbmMuchHitFinderQa.h:81
CbmMuchSector.h
CbmMuchHitFinderQa::fChargeHistos
TObjArray * fChargeHistos
Definition: CbmMuchHitFinderQa.h:92
CbmMuchHitFinderQa::fSignalPoints
Int_t fSignalPoints
Definition: CbmMuchHitFinderQa.h:76
CbmCluster::GetAddress
Int_t GetAddress() const
Definition: CbmCluster.h:90
CbmMuchHitFinderQa::fFlag
Int_t fFlag
Definition: CbmMuchHitFinderQa.h:80
CbmMuchPoint::PositionOut
void PositionOut(TVector3 &pos) const
Definition: CbmMuchPoint.h:80
CbmMuchPoint::GetZIn
Double_t GetZIn() const
Definition: CbmMuchPoint.h:72
CbmMuchDigi::GetAddress
virtual Int_t GetAddress() const
Definition: CbmMuchDigi.h:77
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
mpv_e
constexpr double mpv_e[]
Definition: CbmMuchRecoDefs.h:23
CbmMuchHitFinderQa::fNstations
Int_t fNstations
Definition: CbmMuchHitFinderQa.h:90
CbmMuchModule
Definition: CbmMuchModule.h:24
CbmMuchAddress::GetLayerIndex
static Int_t GetLayerIndex(Int_t address)
Definition: CbmMuchAddress.h:106
CbmMuchAddress::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchAddress.h:103
CbmMuchHitFinderQa::fClusters
TClonesArray * fClusters
Definition: CbmMuchHitFinderQa.h:85
CbmMuchHitFinderQa::PullsQa
void PullsQa()
Definition: CbmMuchHitFinderQa.cxx:1325
CbmMuchPad.h
CbmMuchDigi
Definition: CbmMuchDigi.h:31
CbmMuchModule::GetDetectorType
Int_t GetDetectorType() const
Definition: CbmMuchModule.h:52
CbmMuchPad
Definition: CbmMuchPad.h:21
mpv_mu
constexpr double mpv_mu[]
Definition: CbmMuchRecoDefs.h:25
CbmMuchPointInfo::GetKine
Double_t GetKine()
Definition: CbmMuchPointInfo.h:34
CbmMCTrack.h
CbmMuchHitFinderQa::fhPadsTotalR
TH1D ** fhPadsTotalR
Definition: CbmMuchHitFinderQa.h:108
CbmMuchAddress.h
CbmMuchHitFinderQa::fhChargeEnergyLogPi
TH2D * fhChargeEnergyLogPi
Definition: CbmMuchHitFinderQa.h:94
CbmDigiManager.h
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmMuchHitFinderQa::GetMaxPadSize
TVector2 GetMaxPadSize(Int_t iStation)
Definition: CbmMuchHitFinderQa.cxx:1600
CbmMuchPoint::GetYIn
Double_t GetYIn() const
Definition: CbmMuchPoint.h:71
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmMuchHitFinderQa::fhCharge
TH1D ** fhCharge
Definition: CbmMuchHitFinderQa.h:106
CbmMuchHitFinderQa::fMCTracks
TClonesArray * fMCTracks
Definition: CbmMuchHitFinderQa.h:87
CbmMuchPixelHit.h
Class for pixel hits in MUCH detector.
CbmMuchHitFinderQa::fPointsTotal
Int_t fPointsTotal
Definition: CbmMuchHitFinderQa.h:144
CbmMuchHitFinderQa::fClusterDeconvQaOn
Bool_t fClusterDeconvQaOn
Definition: CbmMuchHitFinderQa.h:152
CbmMuchHitFinderQa::StatisticsQa
void StatisticsQa()
Definition: CbmMuchHitFinderQa.cxx:1251
CbmMuchHitFinderQa::fhResidualX
TH1D * fhResidualX
Definition: CbmMuchHitFinderQa.h:127
CbmMuchHitFinderQa::fhChargeTrackLengthPr
TH2D * fhChargeTrackLengthPr
Definition: CbmMuchHitFinderQa.h:99
CbmMuchPoint::GetXIn
Double_t GetXIn() const
Definition: CbmMuchPoint.h:70
CbmMuchHitFinderQa::Init
virtual InitStatus Init()
Definition: CbmMuchHitFinderQa.cxx:142
CbmMuchHitFinderQa::fhHitsInCluster
TH1I ** fhHitsInCluster
Definition: CbmMuchHitFinderQa.h:133
CbmMuchPad::GetX
Double_t GetX() const
Definition: CbmMuchPad.h:28
CbmMuchHitFinderQa::DigitizerQa
void DigitizerQa()
Definition: CbmMuchHitFinderQa.cxx:1059
CbmMuchGeoScheme::GetStation
CbmMuchStation * GetStation(Int_t iStation) const
Definition: CbmMuchGeoScheme.cxx:209
CbmMuchPointInfo::AddCharge
void AddCharge(Int_t charge)
Definition: CbmMuchPointInfo.h:27
ECbmModuleId::kMuch
@ kMuch
Muon detection system.
CbmMuchGeoScheme::GetModules
std::vector< CbmMuchModule * > GetModules() const
Definition: CbmMuchGeoScheme.cxx:852
CbmMuchPixelHit
Definition: CbmMuchPixelHit.h:17
CbmMuchHitFinderQa::fhPadsFiredR
TH1D ** fhPadsFiredR
Definition: CbmMuchHitFinderQa.h:109
CbmMuchGeoScheme.h
CbmMuchPointInfo.h
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
CbmMuchStation::GetRmin
Double_t GetRmin() const
Definition: CbmMuchStation.h:51
CbmMuchHitFinderQa::fNsecondary
Int_t * fNsecondary
Definition: CbmMuchHitFinderQa.h:142
CbmMuchGeoScheme::GetModuleByDetId
CbmMuchModule * GetModuleByDetId(Int_t detId) const
Definition: CbmMuchGeoScheme.cxx:278
CbmMuchStation.h
CbmCluster::GetDigi
Int_t GetDigi(Int_t index) const
Get digi at position index.
Definition: CbmCluster.h:76
CbmMuchModuleGem::GetPad
CbmMuchPad * GetPad(Int_t address)
Definition: CbmMuchModuleGem.cxx:79
CbmMuchHitFinderQa::ClusterDeconvQa
void ClusterDeconvQa()
Definition: CbmMuchHitFinderQa.cxx:1491
CbmMuchPointInfo::GetPdgCode
Int_t GetPdgCode()
Definition: CbmMuchPointInfo.h:36
CbmMuchModuleGem::GetNSectors
Int_t GetNSectors() const
Definition: CbmMuchModuleGem.h:57
CbmMuchHitFinderQa::fhOccupancyR
TH1D ** fhOccupancyR
Definition: CbmMuchHitFinderQa.h:107
CbmMuchHitFinderQa::fEvent
Int_t fEvent
Definition: CbmMuchHitFinderQa.h:79
CbmMuchHitFinderQa::fHits
TClonesArray * fHits
Definition: CbmMuchHitFinderQa.h:86
CbmMuchHitFinderQa::fDigitizerQaOn
Bool_t fDigitizerQaOn
Definition: CbmMuchHitFinderQa.h:150
CbmMuchHitFinderQa::GetNSectors
Int_t GetNSectors(Int_t iStation)
Definition: CbmMuchHitFinderQa.cxx:1566
CbmMuchModuleGem.h
CbmMuchPointInfo::GetArea
Double_t GetArea()
Definition: CbmMuchPointInfo.h:40
CbmMuchHitFinderQa::~CbmMuchHitFinderQa
virtual ~CbmMuchHitFinderQa()
Definition: CbmMuchHitFinderQa.cxx:137