CbmRoot
CbmMuchFindHitsGem.cxx
Go to the documentation of this file.
1 /*
2  * CbmMuchFindHitsGem.cxx
3  *
4  * Modified on 08/08/2019 : Hit reconstruction in Event (in time slice) and Time slice mode
5  * Default is time slice (kCbmTimeSlice) and it will run in event mode (kCbmEvent) if find event branch in the tree
6  * @authors Vikas Singhal and Ajit Kumar
7  *
8  //@author : Ekata Nandy (ekata@vecc.gov.in) since 21-06-19- Drift time correction for GEM and RPC have been done separately.
9  */
10 #include "CbmMuchFindHitsGem.h"
11 #include "CbmDigiManager.h"
12 #include "CbmMuchCluster.h"
13 #include "CbmMuchGeoScheme.h"
14 #include "CbmMuchModule.h"
15 #include "CbmMuchModuleGem.h"
16 #include "CbmMuchPad.h"
17 #include "CbmMuchPixelHit.h"
18 #include "FairRootManager.h"
19 #include "TClonesArray.h"
20 #include "TFile.h"
21 #include "TMath.h"
22 #include "TStopwatch.h"
23 //#include "CbmTimeSlice.h"
24 #include "CbmMuchAddress.h"
25 #include "CbmMuchBeamTimeDigi.h"
26 #include "CbmMuchDigi.h"
27 #include <algorithm>
28 #include <iomanip>
29 #include <iostream>
30 
31 using std::cout;
32 using std::endl;
33 using std::fixed;
34 using std::left;
35 using std::multimap;
36 using std::right;
37 using std::setw;
38 using std::vector;
39 
40 // -------------------------------------------------------------------------
41 CbmMuchFindHitsGem::CbmMuchFindHitsGem(const char* digiFileName, Int_t flag)
42  : FairTask("MuchFindHitsGem", 1)
43  , fDigiFile(digiFileName)
44  , fFlag(flag)
45  , fAlgorithm(3)
46  , fClusterSeparationTime(100.)
47  , fThresholdRatio(0.1)
48  , fEvent(0)
49  , fNofTimeslices(0)
50  ,
51  //fDigis(NULL),
52  fEvents(NULL)
53  , fClusterCharges()
54  , fLocalMax()
55  , fClusterPads()
56  , fNeighbours()
57  , fClusters(new TClonesArray("CbmMuchCluster", 1000))
58  , fHits(new TClonesArray("CbmMuchPixelHit", 1000))
59  , fGeoScheme(CbmMuchGeoScheme::Instance())
60  , fDigiIndices()
61  , fFiredPads()
62  ,
63  // fDaq(),
64  // fTimeSlice(NULL),
65  // fDigiData(),
66  fuClusters(0) {}
67 
68 // ----- Private method Init -------------------------------------------
69 InitStatus CbmMuchFindHitsGem::Init() {
70  FairRootManager* ioman = FairRootManager::Instance();
71  //if (fDaq) fTimeSlice = (CbmTimeSlice*) ioman->GetObject("TimeSlice.");
72  //else fDigis = (TClonesArray*) ioman->GetObject("MuchDigi");
73 
74  // --- Digi Manager for reading digis which were stored in vector
75  fDigiManager = CbmDigiManager::Instance();
76  if (bBeamTimeDigi) fDigiManager->UseMuchBeamTimeDigi();
77  fDigiManager->Init();
78 
79  // fDigis will not be used now. Just for checking. Need to remove
80  /*fDigis = (TClonesArray*) ioman->GetObject("MuchDigi");
81  if (! fDigis)
82  fDigis = (TClonesArray*) ioman->GetObject("MuchBeamTimeDigi");
83  if (! fDigis)
84  fDigis = (TClonesArray*) ioman->GetObject("CbmMuchBeamTimeDigi");
85  if (! fDigis)
86  fDigis = (TClonesArray*) ioman->GetObject("CbmMuchDigi");
87  if (! fDigis)
88  LOG(info) << "MuchFindHitsGem: No MuchDigi or MuchBeamTimeDigi or CbmMuchDigi or CbmMuchBeamTimeDigi exist";
89  */
90 
91  // Implementation of Event by event execution after TimeSlice Mode and Event Building should have one Event branch.
92 
93  fEvents = dynamic_cast<TClonesArray*>(
94  FairRootManager::Instance()->GetObject("Event"));
95  if (!fEvents)
96  fEvents = dynamic_cast<TClonesArray*>(
97  FairRootManager::Instance()->GetObject("CbmEvent"));
98 
99  if (!fEvents) {
100  LOG(info) << GetName() << ": No event branch present.";
101  // return kfatal;
102  } else {
103  fEventMode = kTRUE;
104  //fMode=kCbmEvent;
105  LOG(info)
106  << GetName()
107  << "TimeSlice: Event-by-event mode after Event Building selected. ";
108  }
109 
110  ioman->Register("MuchCluster",
111  "Cluster in MUCH",
112  fClusters,
113  IsOutputBranchPersistent("MuchCluster"));
114  ioman->Register("MuchPixelHit",
115  "Hit in MUCH",
116  fHits,
117  IsOutputBranchPersistent("MuchPixelHit"));
118 
119  // Initialize GeoScheme
120  TFile* oldfile = gFile;
121  TFile* file = new TFile(fDigiFile);
122  TObjArray* stations = (TObjArray*) file->Get("stations");
123  file->Close();
124  file->Delete();
125  gFile = oldfile;
126  fGeoScheme->Init(stations, fFlag);
127  return kSUCCESS;
128 }
129 // -------------------------------------------------------------------------
130 
131 // ----- Task execution ------------------------------------------------
132 void CbmMuchFindHitsGem::Exec(Option_t*) {
133  TStopwatch timer;
134  timer.Start();
135  // fDigiData.clear();
136  // Removing SetDaq functionality as Cluster and Hit Finder algorithm is same for both the Time Based and Event Based mode.
137  //if (fDaq) ;
138  //fDigiData = fTimeSlice->GetMuchData();
139  // else {
140  //LOG(debug)<<"Start Reading digi from a module ";
141 
142  /*for (Int_t iDigi = 0; iDigi < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); iDigi++) {
143  //Reading digi from CbmDigiManager which stors digis in vector
144  const auto * digi;
145  if(!bBeamTimeDigi) digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(iDigi);
146  else digi = (CbmMuchBeamTimeDigi*) fDigiManager->Get<CbmMuchBeamTimeDigi>(iDigi);
147  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(iDigi);
148  CbmMuchModule* module = fGeoScheme->GetModuleByDetId(digi->GetAddress()); //AZ
149  //std::cout << module->GetDetectorType() << std::endl; //AZ
150  if (module->GetDetectorType() == 2) continue; //AZ - skip 2-D straws
151  fDigiData.push_back(*digi);
152  }*/
153  //}
154 
155  // Clear output array
156  if (fHits) fHits->Clear();
157  if (fClusters)
158  fClusters->Delete(); //Clear(); // Delete because of memory escape
159 
160  fuClusters = 0;
161  // --- Time-slice mode: process entire array
162 
163  if (!fEventMode) {
164  //if ( fMode == kCbmTimeslice ){
165  ProcessData(nullptr);
166  LOG(info) << setw(20) << left << GetName() << ": processing time is "
167  << timer.RealTime() << " Time Slice Number is " << fNofTimeslices
168  << " digis "
169  << fDigiManager->GetNofDigis(ECbmModuleId::kMuch)
170  //<< "s digis " << fDigis->GetEntriesFast()
171  << " clusters " << fClusters->GetEntriesFast() << " total hits "
172  << fHits->GetEntriesFast();
173  }
174  // --- Event mode: loop over events
175  else {
176  assert(fEvents);
177  Int_t nEvents = fEvents->GetEntriesFast();
178  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
179  CbmEvent* event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
180  assert(event);
181  Int_t nDigis = (event ? event->GetNofData(ECbmDataType::kMuchDigi)
182  : fDigiManager->GetNofDigis(ECbmModuleId::kMuch));
183  //: fDigis->GetEntriesFast() );
184  //if (event) LOG(debug)<<" Timeslice "<< fNofTimeslices <<" event : " << event->GetNumber() <<" nDigi : " << nDigis;
185  ProcessData(event);
186  LOG(debug) << setw(20) << left
187  << GetName()
188  //<< ": Processing Time for an event is " << timer.RealTime()
189  << ": Time slice " << fNofTimeslices << " with " << nEvents
190  << (nEvents == 1 ? " event" : " events")
191  << " and processing event nubmer " << iEvent << " digis "
192  << nDigis
193  //<< "s digis " << fDigis->GetEntriesFast()
194  << " and created cluster "
195  << event->GetNofData(ECbmDataType::kMuchCluster)
196  << " and created hit "
197  << event->GetNofData(ECbmDataType::kMuchPixelHit);
198  } //# events
199  LOG(info) << setw(20) << left << GetName() << ": Processing Time is "
200  << timer.RealTime() << ": Time slice " << fNofTimeslices
201  << " with " << nEvents << (nEvents == 1 ? " event" : " events")
202  << "s digis "
203  << fDigiManager->GetNofDigis(ECbmModuleId::kMuch)
204  //<< "s digis " << fDigis->GetEntriesFast()
205  << " and event wise total "
206  << " clusters " << fClusters->GetEntriesFast() << " total hits "
207  << fHits->GetEntriesFast();
208 
209  } //? event mode
210  fNofTimeslices++;
211 }
212 // -------------------------------------------------------------------------
213 
214 // ----- Public method Exec --------------------------------------------
215 void CbmMuchFindHitsGem::ProcessData(CbmEvent* event) {
216  TStopwatch EventTimer;
217  EventTimer.Start();
218 
219  fEvent++;
220  //LOG(debug3)<<" Start creating cluster ";
221  // Find clusters
222  FindClusters(event);
223  Int_t NuOfClusterInEvent =
224  (event ? event->GetNofData(ECbmDataType::kMuchCluster)
225  : fClusters->GetEntriesFast());
226 
227  for (Int_t clusterIndex = 0; clusterIndex < NuOfClusterInEvent;
228  ++clusterIndex) {
229  UInt_t iCluster =
230  (event ? event->GetIndex(ECbmDataType::kMuchCluster, clusterIndex)
231  : clusterIndex);
232  CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(iCluster);
233  switch (fAlgorithm) {
234  // Hit
235  case 0:
236  // One hit per pad
237  case 1: {
238  // One hit per cluster
239  CreateHits(cluster, iCluster, event);
240  break;
241  }
242  case 2: {
243  // Simple cluster deconvolution
244  ExecClusteringSimple(cluster, iCluster, event);
245  break;
246  }
247  case 3: {
248  ExecClusteringPeaks(cluster, iCluster, event);
249  break;
250  }
251  default: {
252  Fatal("CbmMuchFindHitsGem::Exec:",
253  "The algorithm index does not exist.");
254  break;
255  }
256  }
257  }
258  fDigiIndices.clear();
259  fFiredPads.clear();
260 }
261 // -------------------------------------------------------------------------
262 
263 
264 // ----- Private method FindClusters ------------------------------------
265 void CbmMuchFindHitsGem::FindClusters(CbmEvent* event) {
266 
267  Int_t nDigis = (event ? event->GetNofData(ECbmDataType::kMuchDigi)
268  : fDigiManager->GetNofDigis(ECbmModuleId::kMuch));
269  //: fDigis->GetEntriesFast() );
270  if (event)
271  LOG(debug2) << " Timeslice " << fNofTimeslices
272  << " event : " << event->GetNumber() << " nDigi : " << nDigis;
273  if (nDigis < 0) return;
274  if (fAlgorithm == 0) {
275  for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
276  UInt_t digiIndex =
277  (event ? event->GetIndex(ECbmDataType::kMuchDigi, iDigi) : iDigi);
278  fDigiIndices.clear();
279  fDigiIndices.push_back(digiIndex);
280  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(digiIndex);
281  const CbmMuchDigi* digi;
282  if (!bBeamTimeDigi)
283  digi = static_cast<const CbmMuchDigi*>(
284  fDigiManager->Get<CbmMuchDigi>(digiIndex));
285  else
286  digi = static_cast<const CbmMuchDigi*>(
287  fDigiManager->Get<CbmMuchBeamTimeDigi>(digiIndex));
288  //const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(digiIndex));
289  CbmMuchCluster* cluster =
290  new ((*fClusters)[fuClusters++]) CbmMuchCluster();
291  Int_t address = CbmMuchAddress::GetAddress(
292  CbmMuchAddress::GetStationIndex(digi->GetAddress()),
293  CbmMuchAddress::GetLayerIndex(digi->GetAddress()),
294  CbmMuchAddress::GetLayerSideIndex(digi->GetAddress()),
295  CbmMuchAddress::GetModuleIndex(digi->GetAddress()));
296  cluster->SetAddress(address);
297  cluster->AddDigis(fDigiIndices);
298  // --- In event-by-event mode after event building: register clusters to event using ECbmDataType::kMuchCluster
299  //Uncomment below code
300  if (event) {
301  event->AddData(ECbmDataType::kMuchCluster, fuClusters - 1);
302  } //? Event object
303  }
304  return;
305  }
306 
307  vector<CbmMuchModuleGem*> modules = fGeoScheme->GetGemModules();
308 
309  // Clear array of digis in the modules
310  for (UInt_t m = 0; m < modules.size(); m++)
311  modules[m]->ClearDigis();
312 
313  // Fill array of digis in the modules. Digis are automatically sorted in time
314 
315  for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
316  UInt_t digiIndex =
317  (event ? event->GetIndex(ECbmDataType::kMuchDigi, iDigi) : iDigi);
318  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(digiIndex);
319  //const auto * digi;
320  const CbmMuchDigi* digi;
321  if (!bBeamTimeDigi)
322  digi = static_cast<const CbmMuchDigi*>(
323  fDigiManager->Get<CbmMuchDigi>(digiIndex));
324  else
325  digi = static_cast<const CbmMuchDigi*>(
326  fDigiManager->Get<CbmMuchBeamTimeDigi>(digiIndex));
327  //const CbmMuchDigi* digi =static_cast<const CbmMuchDigi*>(fDigis->At(digiIndex));
328  Double_t time = digi->GetTime();
329  // Double_t chanid = digi->GetChannelId();
330  UInt_t address = digi->GetAddress();
331  // UInt_t adc = digi->GetAdc();
332  fGeoScheme->GetModuleByDetId(address)->AddDigi(time, digiIndex);
333  }
334 
335  // Find clusters module-by-module
336  for (UInt_t m = 0; m < modules.size(); m++) {
337  CbmMuchModuleGem* module = modules[m];
338  multimap<Double_t, Int_t> digis = modules[m]->GetDigis();
339  multimap<Double_t, Int_t>::iterator it, itmin, itmax;
340 
341  // Split module digis into time slices according to fClusterSeparationTime
342  vector<multimap<Double_t, Int_t>::iterator> slices;
343  Double_t tlast = -10000;
344  // slices.push_back(digis.begin());
345  for (it = digis.begin(); it != digis.end(); ++it) {
346  Double_t t = it->first;
347  if (t > tlast + fClusterSeparationTime) slices.push_back(it);
348  tlast = t;
349  }
350  slices.push_back(it);
351  for (UInt_t s = 1; s < slices.size(); s++) {
352  fFiredPads.clear();
353  for (it = slices[s - 1]; it != slices[s]; it++) {
354  Int_t iDigi = it->second;
355  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(iDigi);
356  const CbmMuchDigi* digi;
357  //const auto * digi;
358  if (!bBeamTimeDigi)
359  digi = static_cast<const CbmMuchDigi*>(
360  fDigiManager->Get<CbmMuchDigi>(iDigi));
361  else
362  digi = static_cast<const CbmMuchDigi*>(
363  fDigiManager->Get<CbmMuchBeamTimeDigi>(iDigi));
364  //const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(iDigi));
365  CbmMuchPad* pad = module->GetPad(digi->GetAddress());
366  pad->SetDigiIndex(iDigi);
367  fFiredPads.push_back(pad);
368  }
369  // Loop over fired pads in a time slice of 100 ns
370  for (UInt_t p = 0; p < fFiredPads.size(); p++) {
371  fDigiIndices.clear();
372  CreateCluster(fFiredPads[p]);
373  if (fDigiIndices.size() == 0) continue;
374  //const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(fDigiIndices.front()));
375  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(fDigiIndices.front());
376  //const auto * digi;
377  const CbmMuchDigi* digi;
378  if (!bBeamTimeDigi)
379  digi = static_cast<const CbmMuchDigi*>(
380  fDigiManager->Get<CbmMuchDigi>(fDigiIndices.front()));
381  else
382  digi = static_cast<const CbmMuchDigi*>(
383  fDigiManager->Get<CbmMuchBeamTimeDigi>(fDigiIndices.front()));
384 
385  CbmMuchCluster* cluster =
386  new ((*fClusters)[fuClusters++]) CbmMuchCluster();
387  Int_t address = CbmMuchAddress::GetAddress(
388  CbmMuchAddress::GetStationIndex(digi->GetAddress()),
389  CbmMuchAddress::GetLayerIndex(digi->GetAddress()),
390  CbmMuchAddress::GetLayerSideIndex(digi->GetAddress()),
391  CbmMuchAddress::GetModuleIndex(digi->GetAddress()));
392 
393  cluster->SetAddress(address);
394  cluster->AddDigis(fDigiIndices);
395  // --- In event-by-event mode after event building: register clusters to event using ECbmDataType::kMuchCluster
396  if (event) {
397  event->AddData(ECbmDataType::kMuchCluster, fuClusters - 1);
398  } //? Event object
399  }
400  }
401  }
402 }
403 // -------------------------------------------------------------------------
404 
405 
406 // ----- Private method CreateCluster -----------------------------------
407 void CbmMuchFindHitsGem::CreateCluster(CbmMuchPad* pad) {
408  Int_t digiIndex = pad->GetDigiIndex();
409  if (digiIndex < 0) return;
410  fDigiIndices.push_back(digiIndex);
411  pad->SetDigiIndex(-1);
412  vector<CbmMuchPad*> neighbours = pad->GetNeighbours();
413  for (UInt_t i = 0; i < neighbours.size(); i++)
414  CreateCluster(neighbours[i]);
415 }
416 // -------------------------------------------------------------------------
417 
418 
419 // ----- Private method ExecClusteringSimple ----------------------------
420 void CbmMuchFindHitsGem::ExecClusteringSimple(CbmMuchCluster* cluster,
421  Int_t iCluster,
422  CbmEvent* event) {
423  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(cluster->GetDigi(0));
424  //const auto * digi;
425  const CbmMuchDigi* digi;
426  if (!bBeamTimeDigi)
427  digi = static_cast<const CbmMuchDigi*>(
428  fDigiManager->Get<CbmMuchDigi>(cluster->GetDigi(0)));
429  else
430  digi = static_cast<const CbmMuchDigi*>(
431  fDigiManager->Get<CbmMuchBeamTimeDigi>(cluster->GetDigi(0)));
432  //CbmMuchDigi* digi = static_cast<CbmMuchDigi*>(fDigis->At(cluster->GetDigi(0)));
433  CbmMuchModule* m = fGeoScheme->GetModuleByDetId(digi->GetAddress());
434  CbmMuchModuleGem* module = (CbmMuchModuleGem*) m;
435  // Int_t iStation = CbmMuchAddress::GetStationIndex(digi->GetAddress());
436 
437  Int_t maxCharge = 0;
438  for (Int_t iDigi = 0; iDigi < cluster->GetNofDigis(); iDigi++) {
439  Int_t digiIndex = cluster->GetDigi(iDigi);
440  //digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(digiIndex);
441  if (!bBeamTimeDigi)
442  digi = static_cast<const CbmMuchDigi*>(
443  fDigiManager->Get<CbmMuchDigi>(digiIndex));
444  else
445  digi = static_cast<const CbmMuchDigi*>(
446  fDigiManager->Get<CbmMuchBeamTimeDigi>(digiIndex));
447  //digi = static_cast<CbmMuchDigi*> (fDigis->At(digiIndex));
448  Int_t charge = digi->GetAdc();
449  if (charge > maxCharge) maxCharge = charge;
450  }
451 
452  UInt_t threshold = UInt_t(fThresholdRatio * maxCharge);
453 
454  // Fire pads which are higher than threshold in a cluster
455  fFiredPads.clear();
456  for (Int_t iDigi = 0; iDigi < cluster->GetNofDigis(); iDigi++) {
457  Int_t digiIndex = cluster->GetDigi(iDigi);
458  //digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(digiIndex);
459  if (!bBeamTimeDigi)
460  digi = static_cast<const CbmMuchDigi*>(
461  fDigiManager->Get<CbmMuchDigi>(digiIndex));
462  else
463  digi = static_cast<const CbmMuchDigi*>(
464  fDigiManager->Get<CbmMuchBeamTimeDigi>(digiIndex));
465  //digi = static_cast<CbmMuchDigi*> (fDigis->At(digiIndex));
466  if (digi->GetAdc() <= threshold) continue;
467  CbmMuchPad* pad = module->GetPad(digi->GetAddress());
468  pad->SetDigiIndex(digiIndex);
469  fFiredPads.push_back(pad);
470  }
471  for (UInt_t p = 0; p < fFiredPads.size(); p++) {
472  fDigiIndices.clear();
473  CreateCluster(fFiredPads[p]);
474  if (fDigiIndices.size() == 0) continue;
475  CbmMuchCluster cl;
476  cl.AddDigis(fDigiIndices);
477  CreateHits(&cl, iCluster, event);
478  }
479 }
480 // -------------------------------------------------------------------------
481 
482 
483 // -------------------------------------------------------------------------
484 void CbmMuchFindHitsGem::ExecClusteringPeaks(CbmMuchCluster* cluster,
485  Int_t iCluster,
486  CbmEvent* event) {
487  Int_t nDigis = cluster->GetNofDigis();
488  if (nDigis <= 2) {
489  CreateHits(cluster, iCluster, event);
490  return;
491  }
492  fClusterCharges.clear();
493  fClusterPads.clear();
494  fLocalMax.clear();
495  // for (Int_t i=0;i<fNeighbours.size();i++) fNeighbours[i].clear();
496  fNeighbours.clear();
497 
498  // Fill cluster map
499  for (Int_t i = 0; i < nDigis; i++) {
500  Int_t iDigi = cluster->GetDigi(i);
501  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(iDigi);
502  //const auto * digi;
503  const CbmMuchDigi* digi;
504  if (!bBeamTimeDigi)
505  digi =
506  static_cast<const CbmMuchDigi*>(fDigiManager->Get<CbmMuchDigi>(iDigi));
507  else
508  digi = static_cast<const CbmMuchDigi*>(
509  fDigiManager->Get<CbmMuchBeamTimeDigi>(iDigi));
510  //const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(iDigi));
511  UInt_t address = digi->GetAddress();
512  CbmMuchModuleGem* module =
513  (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(address);
514  CbmMuchPad* pad = module->GetPad(address);
515  Int_t adc = digi->GetAdc();
516  fClusterPads.push_back(pad);
517  fClusterCharges.push_back(adc);
518  fLocalMax.push_back(1);
519  }
520 
521  // Fill neighbours
522  for (Int_t i = 0; i < nDigis; i++) {
523  CbmMuchPad* pad = fClusterPads[i];
524  vector<CbmMuchPad*> neighbours = pad->GetNeighbours();
525  vector<Int_t> selected_neighbours;
526  for (UInt_t ip = 0; ip < neighbours.size(); ip++) {
527  CbmMuchPad* p = neighbours[ip];
528  vector<CbmMuchPad*>::iterator it =
529  find(fClusterPads.begin(), fClusterPads.end(), p);
530  if (it == fClusterPads.end()) continue;
531  selected_neighbours.push_back(it - fClusterPads.begin());
532  }
533  fNeighbours.push_back(selected_neighbours);
534  }
535 
536  // Flag local maxima
537  for (Int_t i = 0; i < nDigis; i++) {
538  Int_t c = fClusterCharges[i];
539  for (UInt_t n = 0; n < fNeighbours[i].size(); n++) {
540  Int_t in = fNeighbours[i][n];
541  Int_t cn = fClusterCharges[in];
542  if (cn < c) fLocalMax[in] = 0;
543  }
544  }
545 
546  // Fire pads corresponding to local maxima
547  fFiredPads.clear();
548  for (Int_t i = 0; i < nDigis; i++) {
549  if (fLocalMax[i] == 0) continue;
550  CbmMuchPad* pad = fClusterPads[i];
551  pad->SetDigiIndex(cluster->GetDigi(i));
552  fFiredPads.push_back(pad);
553  }
554 
555  // Create clusters
556  for (UInt_t p = 0; p < fFiredPads.size(); p++) {
557  fDigiIndices.clear();
558  CreateCluster(fFiredPads[p]);
559  if (fDigiIndices.size() == 0) continue;
560  CbmMuchCluster cl;
561  cl.AddDigis(fDigiIndices);
562  CreateHits(&cl, iCluster, event);
563  }
564 }
565 // -------------------------------------------------------------------------
566 
567 
568 // ----- Private method CreateHits --------------------------------------
569 void CbmMuchFindHitsGem::CreateHits(CbmMuchCluster* cluster,
570  Int_t iCluster,
571  CbmEvent* event) {
572  Int_t nDigis = cluster->GetNofDigis();
573  Double_t sumq = 0, sumx = 0, sumy = 0, sumt = 0, sumdx2 = 0, sumdy2 = 0,
574  sumdxy2 = 0, sumdt2 = 0;
575  Double_t q = 0, x = 0, y = 0, t = 0, z = 0, dx = 0, dy = 0, dxy = 0, dt = 0;
576  Double_t nX = 0, nY = 0, nZ = 0;
577  Int_t address = 0;
578  Int_t planeId = 0;
579  CbmMuchModuleGem* module = NULL;
580 
581  Double_t tmin = -1;
582 
583  for (Int_t i = 0; i < nDigis; i++) {
584  Int_t iDigi = cluster->GetDigi(i);
585  //const CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(iDigi);
586  //const auto * digi;
587  const CbmMuchDigi* digi;
588  if (!bBeamTimeDigi)
589  digi =
590  static_cast<const CbmMuchDigi*>(fDigiManager->Get<CbmMuchDigi>(iDigi));
591  else
592  digi = static_cast<const CbmMuchDigi*>(
593  fDigiManager->Get<CbmMuchBeamTimeDigi>(iDigi));
594  //const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(iDigi));
595  if (i == 0) {
596  address =
597  CbmMuchAddress::GetElementAddress(digi->GetAddress(), kMuchModule);
598  planeId = fGeoScheme->GetLayerSideNr(address);
599  module = (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(address);
600  z = module->GetPosition()[2];
601  }
602 
603  CbmMuchModule* moduleDet =
604  fGeoScheme->GetModuleByDetId(digi->GetAddress());
605 
606  CbmMuchPad* pad = module->GetPad(digi->GetAddress());
607  x = pad->GetX();
608  y = pad->GetY();
609 
611  Double_t gemDriftTimeCorrc =
612  15.0; // Drift time mean for GEM is 15 ns. Drift vel =100um/ns and drift width 3mm.
613  Double_t rpcDriftTimeCorrc =
614  8.33; // Drift time mean for RPC is 8.33 ns. Drift vel =120um/ns and drift width 2mm.
615  Double_t gemHitTimeError =
616  4.0; // Hit time error for GEM = 4.0 as residual dist width is 4, to make pull width 1.
617  Double_t rpcHitTimeError =
618  2.3; // Hit time error for RPC = 2.3 ns, as residual dist width is 2.3. This is made to make pull dist width ~1
619  Double_t timeShift =
620  0.5; // this is added because residual time dist is shifted by -0.5 from 0.
621 
622  if (moduleDet->GetDetectorType() == 3)
623  {
624  if (fFlag == 0)
625  t = digi->GetTime() - gemDriftTimeCorrc + timeShift;
626  else
627  t = digi->GetTime(); // Not correcting Drift Time for mCBM data
628  dt = gemHitTimeError;
629  }
630  if (moduleDet->GetDetectorType() == 4)
631  {
632  t = digi->GetTime() - rpcDriftTimeCorrc + timeShift;
633  dt = rpcHitTimeError;
634  }
635 
636  if (tmin < 0) tmin = t;
637  if (tmin < t) tmin = t;
638  q = digi->GetAdc();
639  dx = pad->GetDx();
640  dy = pad->GetDy();
641  dxy = pad->GetDxy();
642  //dt = 4.; // digi->GetDTime(); //TODO introduce time uncertainty determination
643  sumq += q;
644  sumx += q * x;
645  sumy += q * y;
646  sumt += q * t;
647  sumdx2 += q * q * dx * dx;
648  sumdy2 += q * q * dy * dy;
649  sumdxy2 += q * q * dxy * dxy;
650  sumdt2 += q * q * dt * dt;
651  //std::cout<<" i "<<i<<" q "<<q<<" sumq "<<sumq<<std::endl;
652  }
653 
654  x = sumx / sumq;
655  y = sumy / sumq;
656  // t = sumt/sumq;
657  t = tmin;
658  dx = sqrt(sumdx2 / 12) / sumq;
659  dy = sqrt(sumdy2 / 12) / sumq;
660  dxy = sqrt(sumdxy2 / 12) / sumq;
661  dt = sqrt(sumdt2) / sumq;
662  Int_t iHit = fHits->GetEntriesFast();
663 
664  //------------------------------Added by O. Singh 11.12.2017 for mCbm ---------------------------
665  Double_t tX = 18.5, tY = 80.5;
666  nX = x + tX; // Ajit + OS + Apar -> For miniMUCH setup in March 2019
667  nY = y + tY; // Ajit + OS + Apar -> For miniMUCH setup in March 2019
668  nZ = z;
669 
670  if (fFlag == 1) {
671  new ((*fHits)[iHit]) CbmMuchPixelHit(
672  address, nX, nY, nZ, dx, dy, 0, dxy, iCluster, planeId, t, dt); //mCbm
673  } else {
674  new ((*fHits)[iHit]) CbmMuchPixelHit(
675  address, x, y, z, dx, dy, 0, dxy, iCluster, planeId, t, dt); //Cbm
676  }
677  //Adding CbmMuchPixelHit entries in the CbmEvent
678  if (event) {
679  event->AddData(ECbmDataType::kMuchPixelHit, iHit);
680  } //? Event object
681 }
682 // ---------------------------------------------------------------------------------
683 
684 ClassImp(CbmMuchFindHitsGem)
CbmMuchModule.h
CbmMuchPad::GetNeighbours
std::vector< CbmMuchPad * > GetNeighbours() const
Definition: CbmMuchPad.h:38
CbmMuchDigi.h
CbmMuchFindHitsGem::fGeoScheme
CbmMuchGeoScheme * fGeoScheme
Definition: CbmMuchFindHitsGem.h:93
CbmMuchGeoScheme
Definition: CbmMuchGeoScheme.h:43
CbmMuchDigi::GetTime
virtual Double_t GetTime() const
Definition: CbmMuchDigi.h:78
CbmMuchFindHitsGem::fClusters
TClonesArray * fClusters
Definition: CbmMuchFindHitsGem.h:91
CbmMuchFindHitsGem::Exec
virtual void Exec(Option_t *opt)
Definition: CbmMuchFindHitsGem.cxx:132
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition: CbmDigiManager.cxx:71
CbmMuchCluster
Data container for MUCH clusters.
Definition: CbmMuchCluster.h:20
CbmMuchFindHitsGem::fAlgorithm
Int_t fAlgorithm
Definition: CbmMuchFindHitsGem.h:76
ECbmDataType::kMuchPixelHit
@ kMuchPixelHit
CbmMuchPad::GetY
Double_t GetY() const
Definition: CbmMuchPad.h:29
CbmMuchBeamTimeDigi
Definition: CbmMuchBeamTimeDigi.h:29
CbmMuchDigi::GetAdc
UShort_t GetAdc() const
Definition: CbmMuchDigi.h:74
CbmMuchFindHitsGem::ExecClusteringSimple
void ExecClusteringSimple(CbmMuchCluster *cluster, Int_t iCluster, CbmEvent *event)
Definition: CbmMuchFindHitsGem.cxx:420
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmMuchFindHitsGem::fClusterSeparationTime
Double_t fClusterSeparationTime
Definition: CbmMuchFindHitsGem.h:78
CbmMuchFindHitsGem::fEventMode
Bool_t fEventMode
Definition: CbmMuchFindHitsGem.h:97
CbmMuchPad::GetDy
Double_t GetDy() const
Definition: CbmMuchPad.h:31
CbmDigiManager::GetNofDigis
static Int_t GetNofDigis(ECbmModuleId systemId)
Definition: CbmDigiManager.cxx:62
CbmDigiManager::UseMuchBeamTimeDigi
void UseMuchBeamTimeDigi(Bool_t choice=kTRUE)
Use CbmMuchBeamTimeDigi instead of CbmMuchDigi for MUCH.
Definition: CbmDigiManager.h:130
CbmMuchGeoScheme::Init
void Init(TObjArray *stations, Int_t flag)
Definition: CbmMuchGeoScheme.cxx:121
CbmMuchFindHitsGem::fNeighbours
std::vector< std::vector< Int_t > > fNeighbours
Definition: CbmMuchFindHitsGem.h:89
CbmMuchGeoScheme::GetGemModules
std::vector< CbmMuchModuleGem * > GetGemModules() const
Definition: CbmMuchGeoScheme.cxx:872
CbmMuchPad::GetDx
Double_t GetDx() const
Definition: CbmMuchPad.h:30
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
CbmMuchModule::GetPosition
TVector3 GetPosition() const
Definition: CbmMuchModule.h:51
CbmMuchCluster.h
Data container for MUCH clusters.
CbmMuchAddress::GetModuleIndex
static Int_t GetModuleIndex(Int_t address)
Definition: CbmMuchAddress.h:112
CbmMuchPad::GetDxy
Double_t GetDxy() const
Definition: CbmMuchPad.h:32
CbmMuchPad::GetDigiIndex
Int_t GetDigiIndex() const
Definition: CbmMuchPad.h:34
CbmMuchFindHitsGem::fuClusters
UInt_t fuClusters
Definition: CbmMuchFindHitsGem.h:107
CbmMuchFindHitsGem::fEvents
TClonesArray * fEvents
Interface to digi branch.
Definition: CbmMuchFindHitsGem.h:85
CbmMuchPad::SetDigiIndex
void SetDigiIndex(Int_t iDigi)
Definition: CbmMuchPad.h:43
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmMuchFindHitsGem::fClusterCharges
std::vector< Int_t > fClusterCharges
Definition: CbmMuchFindHitsGem.h:86
CbmMuchFindHitsGem::fDigiIndices
std::vector< Int_t > fDigiIndices
Definition: CbmMuchFindHitsGem.h:95
CbmMuchModuleGem
Definition: CbmMuchModuleGem.h:24
CbmMuchBeamTimeDigi.h
CbmMuchFindHitsGem.h
CbmMuchDigi::GetAddress
virtual Int_t GetAddress() const
Definition: CbmMuchDigi.h:77
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmMuchFindHitsGem::bBeamTimeDigi
Bool_t bBeamTimeDigi
Definition: CbmMuchFindHitsGem.h:108
CbmMuchFindHitsGem::ExecClusteringPeaks
void ExecClusteringPeaks(CbmMuchCluster *cluster, Int_t iCluster, CbmEvent *event)
Definition: CbmMuchFindHitsGem.cxx:484
CbmMuchFindHitsGem::fHits
TClonesArray * fHits
Definition: CbmMuchFindHitsGem.h:92
CbmMuchFindHitsGem::fDigiFile
TString fDigiFile
Definition: CbmMuchFindHitsGem.h:74
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
CbmMuchPad.h
CbmMuchDigi
Definition: CbmMuchDigi.h:31
CbmMuchModule::GetDetectorType
Int_t GetDetectorType() const
Definition: CbmMuchModule.h:52
CbmMuchPad
Definition: CbmMuchPad.h:21
CbmMuchFindHitsGem::fLocalMax
std::vector< Bool_t > fLocalMax
Definition: CbmMuchFindHitsGem.h:87
CbmCluster::SetAddress
void SetAddress(Int_t address)
Definition: CbmCluster.h:94
CbmMuchAddress::GetElementAddress
static Int_t GetElementAddress(Int_t address, Int_t level)
Definition: CbmMuchAddress.h:122
CbmMuchFindHitsGem::fDigiManager
CbmDigiManager * fDigiManager
Definition: CbmMuchFindHitsGem.h:84
ECbmDataType::kMuchCluster
@ kMuchCluster
CbmMuchAddress.h
CbmDigiManager.h
CbmMuchFindHitsGem::fEvent
Int_t fEvent
Definition: CbmMuchFindHitsGem.h:81
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmMuchFindHitsGem
Definition: CbmMuchFindHitsGem.h:46
CbmMuchFindHitsGem::fThresholdRatio
Double_t fThresholdRatio
Definition: CbmMuchFindHitsGem.h:80
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
CbmMuchPixelHit.h
Class for pixel hits in MUCH detector.
CbmMuchFindHitsGem::fFiredPads
std::vector< CbmMuchPad * > fFiredPads
Definition: CbmMuchFindHitsGem.h:96
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition: CbmEvent.h:30
CbmMuchPad::GetX
Double_t GetX() const
Definition: CbmMuchPad.h:28
CbmMuchModule::GetDigis
std::multimap< Double_t, Int_t > GetDigis()
Definition: CbmMuchModule.h:70
ECbmModuleId::kMuch
@ kMuch
Muon detection system.
kMuchModule
@ kMuchModule
Module.
Definition: CbmMuchAddress.h:20
CbmMuchFindHitsGem::fFlag
Int_t fFlag
Definition: CbmMuchFindHitsGem.h:75
CbmMuchPixelHit
Definition: CbmMuchPixelHit.h:17
CbmMuchGeoScheme.h
CbmCluster::AddDigis
void AddDigis(const std::vector< Int_t > &indices)
Add array of digi to cluster.
Definition: CbmCluster.h:53
CbmMuchFindHitsGem::fClusterPads
std::vector< CbmMuchPad * > fClusterPads
Definition: CbmMuchFindHitsGem.h:88
CbmMuchFindHitsGem::fNofTimeslices
Int_t fNofTimeslices
Definition: CbmMuchFindHitsGem.h:82
CbmMuchFindHitsGem::Init
virtual InitStatus Init()
Definition: CbmMuchFindHitsGem.cxx:69
CbmMuchAddress::GetAddress
static UInt_t GetAddress(Int_t station=0, Int_t layer=0, Int_t side=0, Int_t module=0, Int_t sector=0, Int_t channel=0)
Definition: CbmMuchAddress.cxx:43
CbmMuchFindHitsGem::ProcessData
void ProcessData(CbmEvent *)
Definition: CbmMuchFindHitsGem.cxx:215
CbmMuchFindHitsGem::CbmMuchFindHitsGem
CbmMuchFindHitsGem(const char *digiFileName, Int_t flag)
Definition: CbmMuchFindHitsGem.cxx:41
CbmMuchGeoScheme::GetModuleByDetId
CbmMuchModule * GetModuleByDetId(Int_t detId) const
Definition: CbmMuchGeoScheme.cxx:278
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
CbmMuchModule::AddDigi
void AddDigi(Double_t time, Int_t id)
Definition: CbmMuchModule.h:64
CbmMuchFindHitsGem::CreateHits
void CreateHits(CbmMuchCluster *cluster, Int_t iCluster, CbmEvent *event)
Definition: CbmMuchFindHitsGem.cxx:569
ECbmDataType::kMuchDigi
@ kMuchDigi
CbmMuchModuleGem.h
CbmMuchAddress::GetLayerSideIndex
static Int_t GetLayerSideIndex(Int_t address)
Definition: CbmMuchAddress.h:109
CbmMuchGeoScheme::GetLayerSideNr
Int_t GetLayerSideNr(Int_t detId) const
Definition: CbmMuchGeoScheme.cxx:403
CbmMuchFindHitsGem::CreateCluster
void CreateCluster(CbmMuchPad *pad)
Definition: CbmMuchFindHitsGem.cxx:407
CbmMuchFindHitsGem::FindClusters
void FindClusters(CbmEvent *)
Definition: CbmMuchFindHitsGem.cxx:265