FairRoot/PandaRoot
PndGemFindClustersTB.cxx
Go to the documentation of this file.
1 /* $Id: */
2 
3 // -------------------------------------------------------------------------
4 // ----- PndGemFindClustersTB source file -----
5 // ----- Created 15/06/2013 by R. Karabowicz -----
6 // -------------------------------------------------------------------------
7 
18 #include "PndGemFindClustersTB.h"
19 
20 #include "PndGemCluster.h"
21 #include "PndGemDigi.h"
22 #include "PndGemDigiPar.h"
23 #include "PndGemSensor.h"
24 #include "PndGemStation.h"
25 
26 #include "PndDetectorList.h"
27 
28 #include "FairRootManager.h"
29 #include "FairRunAna.h"
30 #include "FairRuntimeDb.h"
31 
32 #include "TClonesArray.h"
33 #include "TMath.h"
34 #include "TRandom2.h"
35 
36 #include <iomanip>
37 
38 using std::cout;
39 using std::cerr;
40 using std::endl;
41 using std::flush;
42 using std::fixed;
43 using std::right;
44 using std::left;
45 using std::setw;
46 using std::setprecision;
47 using std::set;
48 using std::map;
49 
50 // ----- Default constructor ------------------------------------------
51 PndGemFindClustersTB::PndGemFindClustersTB()
52  : FairTask("GEMFindClusters", 0),
53  fDigiPar (NULL),
54  fDigis (NULL),
55  fClusters (NULL),
56  fTNofEvents (0),
57  fTNofDigis (0),
58  fTNofClusters(0),
59  fFunctor (NULL),
60  fPrepTime (0.),
61  fCreateTime (0.),
62  fSortTime (0.),
63  fCheckTime (0.),
64  fAnaTime (0.),
65  fWriteTime (0.),
66  fAllTime (0.),
67  fInBranchName("GEMDigi")
68 {
69 }
70 // -------------------------------------------------------------------------
71 
72 
73 
74 // ----- Standard constructor ------------------------------------------
75 PndGemFindClustersTB::PndGemFindClustersTB(Int_t iVerbose)
76  : FairTask("GEMFindClusters", iVerbose),
77  fDigiPar (NULL),
78  fDigis (NULL),
79  fClusters (NULL),
80  fTNofEvents (0),
81  fTNofDigis (0),
82  fTNofClusters(0),
83  fFunctor (NULL),
84  fPrepTime (0.),
85  fCreateTime (0.),
86  fSortTime (0.),
87  fCheckTime (0.),
88  fAnaTime (0.),
89  fWriteTime (0.),
90  fAllTime (0.),
91  fInBranchName("GEMDigi")
92 {
93 }
94 // -------------------------------------------------------------------------
95 
96 
97 
98 // ----- Constructor with name -----------------------------------------
99 PndGemFindClustersTB::PndGemFindClustersTB(const char* name, Int_t iVerbose)
100  : FairTask(name, iVerbose),
101  fDigiPar (NULL),
102  fDigis (NULL),
103  fClusters (NULL),
104  fTNofEvents (0),
105  fTNofDigis (0),
106  fTNofClusters(0),
107  fFunctor (NULL),
108  fPrepTime (0.),
109  fCreateTime (0.),
110  fSortTime (0.),
111  fCheckTime (0.),
112  fAnaTime (0.),
113  fWriteTime (0.),
114  fAllTime (0.),
115  fInBranchName("GEMDigi") {
116 }
117 // -------------------------------------------------------------------------
118 
119 
120 
121 // ----- Destructor ----------------------------------------------------
122 PndGemFindClustersTB::~PndGemFindClustersTB() {
123  if ( 0!=fClusters ) {
124  fClusters->Delete();
125  delete fClusters;
126  }
127  if ( 0!=fFunctor ) delete fFunctor;
128 }
129 // -------------------------------------------------------------------------
130 
131 
132 
133 // ----- Public method Exec --------------------------------------------
134 void PndGemFindClustersTB::Exec(Option_t*) {
135  // cout << endl << "======== PndGemFindClustersTB::Exec(Event = " << fTNofEvents << " ) ====================" << endl;
136  fTimer.Start();
137 
138  if ( fInBranchName.Contains("Sorted") ) {
139  fDigis->Clear();
140  fDigis = FairRootManager::Instance()->GetData(fInBranchName, fFunctor, 20.); //FairRootManager::Instance()->GetEventTime() +
141  if(fVerbose > 1)
142  std::cout << "-I- PndGemFindClustersTB::Exec Digis: " << fDigis->GetEntries() << std::endl;
143  }
144  else
145  fDigis = (TClonesArray*)FairRootManager::Instance()->GetObject(fInBranchName);
146 
147  fClusters = FairRootManager::Instance()->GetTClonesArray("GEMCluster");
148  if ( ! fClusters ) Fatal("Exec", "No GEM Cluster Array");
149  fClusters->Delete();
150 
151  fPrepTime+=fTimer.RealTime();
152  fTimer.Continue();
153 
154  Int_t nDigis = fDigis->GetEntries();
155  Int_t nClusters = 0;
156 
157  CreateClusters();
158  fCreateTime+=fTimer.RealTime();
159  fTimer.Continue();
160 
161  Int_t nofC2W = 0;
162  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ )
163  if ( fDigiClusters[idc].cluPos > -0.5 )
164  nofC2W++;
165  // cout << ">>>> created " << fDigiClusters.size() << " clusters, " << nofC2W << " to write" << endl;
166 
167  // cout << "sorting" << endl;
168  SortClusters();
169  fSortTime+=fTimer.RealTime();
170  fTimer.Continue();
171  // cout << "done" << endl;
172 
173  // cout << "checking" << endl;
174  CheckClusters();
175  fCheckTime+=fTimer.RealTime();
176  fTimer.Continue();
177  // cout << "done" << endl;
178 
179  // cout << "analyzing" << endl;
180  AnalyzeClusters();
181  fAnaTime+=fTimer.RealTime();
182  fTimer.Continue();
183  // cout << "done" << endl;
184 
185  nofC2W = 0;
186  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ )
187  if ( fDigiClusters[idc].cluPos > -0.5 )
188  nofC2W++;
189  // cout << ">>>> after cleaning " << fDigiClusters.size() << " clusters, " << nofC2W << " to write" << endl;
190 
191  // PrintClusters();
192 
193  nClusters = WriteClusters();
194  fWriteTime+=fTimer.RealTime();
195  fTimer.Continue();
196 
197  /*
198  // cout << " PndGemFindClustersTB::Exec. Created " << nClusters << " clusters out of " << nDigis << " digis" << flush;
199  // in " << fTimer.RealTime() << "s" << flush;
200  if ( fDigis->GetEntries() > 0 ) {
201  cout << ", ( " << ((PndGemDigi*)fDigis->At(0))->GetTimeStamp()
202  << " -- " << ((PndGemDigi*)fDigis->At(nDigis-1))->GetTimeStamp()
203  << " )" << endl;
204  }
205  else {
206  cout << "." << endl;
207  }
208  */
209 
210  fTNofEvents += 1;
211  fTNofDigis += nDigis;
212  fTNofClusters += nClusters;
213 
214  fAllTime+=fTimer.RealTime();
215 
216  fTimer.Stop();
217 
218 }
219 // -------------------------------------------------------------------------
220 
221 // -------------------------------------------------------------------------
222 Int_t PndGemFindClustersTB::CreateClusters() {
223  fDigiClusters.clear();
224 
225  PndGemDigi* digi = NULL;
226  for ( Int_t idigi = 0 ; idigi < fDigis->GetEntries() ; idigi++ ) {
227  digi = (PndGemDigi*) fDigis->At(idigi);
228  // cout << "digi " << idigi << " in " << digi->GetDetectorId() << ", there are already " << fDigiClusters.size() << " clusters" << endl;
229  // cout << idigi << " ---> " << digi->GetDetectorId() << " @ " << digi->GetChannelNr() << " @ " << digi->GetTimeStamp() << endl;
230 
231  if ( CompareDigiToClusters(idigi) == kTRUE ) continue;
232 
233  DigiClusterTB tempDC;
234  tempDC.detId = digi->GetDetectorId();
235  tempDC.digiNr.push_back(idigi);
236  tempDC.chanNr.push_back((Int_t)digi->GetChannelNr());
237  tempDC.sigADC.push_back(digi->GetCharge());
238  tempDC.sigTDC.push_back(digi->GetTimeStamp());
239  tempDC.cluTDC = digi->GetTimeStamp();
240  tempDC.cluADC = digi->GetCharge();
241  tempDC.cluPos = digi->GetChannelNr();
242  tempDC.cluPMn = digi->GetChannelNr();
243  tempDC.cluPMx = digi->GetChannelNr();
244  tempDC.cluMPs = digi->GetChannelNr();
245  tempDC.cluMVl = digi->GetCharge();
246  fDigiClusters.push_back(tempDC);
247 
248  }
249 
250  // if ( fTNofEvents == 9 ) {
251  // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
252  // if ( fDigiClusters[idc].detId == 568328832 || fDigiClusters[idc].detId == 568329088 ) {
253  // cout << "CRE CLUSTER " << idc << " in " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].cluPMn << " to " << fDigiClusters[idc].cluPMx << " at " << fDigiClusters[idc].cluPos << endl;
254  // }
255  // }
256  // }
257 
258  // cout << endl;
259  return fDigiClusters.size();
260 }
261 // -------------------------------------------------------------------------
262 
263 // ----- Private method CompareDigiToClusters() ------------------------------------------
264 Bool_t PndGemFindClustersTB::CompareDigiToClusters(Int_t digiNumber) {
265  Bool_t printInfo = kFALSE;
266 
267  PndGemSensor* sensor;
268  PndGemDigi* digi = (PndGemDigi*) fDigis->At(digiNumber);
269 
270  // if ( fTNofEvents == 9 && digi->GetDetectorId() == 568328832 && TMath::Abs(digi->GetChannelNr()-3100)<100 ) printInfo = kTRUE;
271 
272  if ( printInfo )
273  cout << "digi [" << digiNumber << "] in detector " << digi->GetDetectorId() << " channel " << digi->GetChannelNr() << " @ " << digi->GetTimeStamp() << endl;
274 
275  Int_t digiUsed = -1;
276  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
277  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
278  if ( digi->GetDetectorId() != fDigiClusters[idc].detId ) continue;
279 
280  // if ( fDigiClusters[idc].detId == 568329024 ) {
281  //if ( fDigiClusters[idc].detId == -1 ) {
282  if ( printInfo ) {
283  cout << "Trying to match digi "
284  << digiNumber << " @ "
285  << digi->GetChannelNr() << " to cluster @ "
286  << fDigiClusters[idc].cluPos << " (from "
287  << fDigiClusters[idc].cluPMn << " to "
288  << fDigiClusters[idc].cluPMx << ") @ "
289  << fDigiClusters[idc].cluTDC << " with "
290  << fDigiClusters[idc].sigADC.size() << " digis" << endl;
291  }
292 
293  Int_t stationNr = digi->GetStationNr();
294  Int_t sensorNr = digi->GetSensorNr();
295  Int_t iSide = digi->GetSide();
296 
297  sensor = fDigiPar->GetSensor(stationNr, sensorNr);
298 
299  // This GetDistance(side,ca,cb,ct) returns:
300  // 0 if ct is between ca and cb
301  // -1 if ct is on different plane that ca/cb
302  // n the distance of ct to ca or cb (smaller one)
303  Int_t digiClustDist = sensor->GetDistance(iSide,fDigiClusters[idc].cluPMn,fDigiClusters[idc].cluPMx,digi->GetChannelNr());
304  Int_t clusterLength = sensor->GetDistance(iSide,fDigiClusters[idc].cluPMn,fDigiClusters[idc].cluPMx)+1;
305  if ( printInfo )
306  cout << "GOT DISTANCE: " << digiClustDist << " WHILE THE CLUSTER LENGTH = " << clusterLength << " CHECK AT " << 19/clusterLength+1 << endl;
307  if ( digiClustDist < 0 ) continue; // there's no connection
308  if ( digiClustDist > 20 ) continue; // they are too far
309  Double_t digiTimeDist = digi->GetTimeStamp() - fDigiClusters[idc].cluTDC;
310  // cout << "TIME DIST: " << digiTimeDist << endl;
311  if ( TMath::Abs(digiTimeDist) > 11. ) continue; // they are too far in time
312 
313  if ( fVerbose > 1 )
314  cout << "digi " << digiNumber << " ( " << digi->GetDetectorId() << " / " << digi->GetChannelNr() << " ) could attach to cluster "
315  << idc << " ( " << fDigiClusters[idc].cluPos << " ) with distance " << digiClustDist << " . " << fDigiClusters[idc].digiNr.size() << " digis in cluster" << flush;
316  // cout << "ADDING DIGI @ " << digi->GetTimeStamp() << " TO CLUSTER[" << idc << "] AT " << fDigiClusters[idc].cluTDC << flush;
317 
318  AddDigiToCluster(digiNumber,idc);
319  // cout << " MOVES IT TO " << fDigiClusters[idc].cluTDC << endl;
320 
321  if ( fVerbose > 1 )
322  cout << " >>> new pos = " << fDigiClusters[idc].cluPos << endl;
323 
324  if ( printInfo ) {
325  // if ( fDigiClusters[idc].detId == 568329024 ) {
326  cout << "DIGI " << digiNumber << " AT " << digi->GetChannelNr() << " MATCHES TO CLUSTER " << endl;
327  }
328 
329  if ( digiUsed >= 0 ) {
330  //cout << "DIGI " << digiNumber << " MATCHES TO CLUSTER " << digiUsed << " and " << idc << endl;
331  JoinTwoClusters(idc,digiUsed);
332  }
333  digiUsed = idc;
334  }
335  return (digiUsed==-1?kFALSE:kTRUE);
336 
337 }
338 // -------------------------------------------------------------------------
339 
340 // ----- Private method JoinTwoClusters ---------------------------------
341 void PndGemFindClustersTB::JoinTwoClusters(Int_t clus1, Int_t clus2) {
342  Bool_t printInfo = kFALSE;
343  // if ( fTNofEvents == 9 && fDigiClusters[clus1].detId == 568328832 ) printInfo = kTRUE;
344 
345  if ( printInfo ) {
346  cout << "JoinTwoClusters-------------------------------------" << endl;
347  SortCluster(clus1);
348  PrintCluster(clus1);
349  SortCluster(clus2);
350  PrintCluster(clus2);
351  }
352 
353  // go through digis in clus2
354  for ( size_t id2 = 0 ; id2 < fDigiClusters[clus2].digiNr.size() ; id2++ ) {
355  Int_t digi2 = fDigiClusters[clus2].digiNr[id2];
356 
357  Int_t addDigi = kTRUE;
358  for ( size_t id1 = 0 ; id1 < fDigiClusters[clus1].digiNr.size() ; id1++ )
359  if ( digi2 == fDigiClusters[clus1].digiNr[id1] )
360  addDigi = kFALSE;
361  //if ( addDigi ) cout << "SHOULD " << (addDigi?"":"NOT ") << "ADD DIGI " << digi2 << " ( " << fDigiClusters[clus2].chanNr[id2] << " )" << endl;
362  if ( addDigi ) AddDigiToCluster(digi2,clus1);
363  }
364  fDigiClusters[clus2].cluPos = -1.;
365 }
366 // -------------------------------------------------------------------------
367 
368 // ----- Private method AddDigiToCluster ---------------------------------
369 void PndGemFindClustersTB::AddDigiToCluster(Int_t digiNr, Int_t clusNr) {
370  PndGemDigi* digi = (PndGemDigi*) fDigis->At(digiNr);
371  Int_t stationNr = digi->GetStationNr();
372  Int_t sensorNr = digi->GetSensorNr();
373  Int_t iSide = digi->GetSide();
374 
375  PndGemSensor* sensor = fDigiPar->GetSensor(stationNr, sensorNr);
376 
377  if ( fDigiClusters[clusNr].sigTDC.size() == 0 ) {
378  // cout << "NEW CLUSTER " << clusNr << " CREATED WITH AddDigiToCluster(" << digiNr << "," << clusNr << ") IN DETECTOR " << digi->GetDetectorId() << endl;
379  fDigiClusters[clusNr].detId = digi->GetDetectorId();
380  fDigiClusters[clusNr].cluPos = digi->GetChannelNr();
381  fDigiClusters[clusNr].cluADC = digi->GetCharge();
382  fDigiClusters[clusNr].cluTDC = digi->GetTimeStamp();
383  fDigiClusters[clusNr].cluPMn = digi->GetChannelNr();
384  fDigiClusters[clusNr].cluPMx = digi->GetChannelNr();
385  fDigiClusters[clusNr].cluMVl = digi->GetCharge();
386  fDigiClusters[clusNr].cluMPs = digi->GetChannelNr();
387  }
388  else {
389  fDigiClusters[clusNr].cluTDC = (fDigiClusters[clusNr].cluTDC*((Double_t)(fDigiClusters[clusNr].sigADC.size()-1.))+digi->GetTimeStamp())/((Double_t)(fDigiClusters[clusNr].sigADC.size()));
390  fDigiClusters[clusNr].cluPos = sensor->GetMeanChannel(iSide,fDigiClusters[clusNr].cluPos,fDigiClusters[clusNr].cluADC,digi->GetChannelNr(),digi->GetCharge());
391  if ( fDigiClusters[clusNr].cluPMn > digi->GetChannelNr() ) fDigiClusters[clusNr].cluPMn = digi->GetChannelNr();
392  if ( fDigiClusters[clusNr].cluPMx < digi->GetChannelNr() ) fDigiClusters[clusNr].cluPMx = digi->GetChannelNr();
393  fDigiClusters[clusNr].cluADC = fDigiClusters[clusNr].cluADC+digi->GetCharge();
394  if ( digi->GetCharge() > fDigiClusters[clusNr].cluMVl ) {
395  fDigiClusters[clusNr].cluMVl = digi->GetCharge();
396  fDigiClusters[clusNr].cluMPs = digi->GetChannelNr();
397  }
398  }
399 
400  fDigiClusters[clusNr].digiNr.push_back(digiNr);
401  fDigiClusters[clusNr].chanNr.push_back((Int_t)digi->GetChannelNr());
402  fDigiClusters[clusNr].sigADC.push_back(digi->GetCharge());
403  fDigiClusters[clusNr].sigTDC.push_back(digi->GetTimeStamp());
404 }
405 // -------------------------------------------------------------------------
406 
407 // ----- Private method AnalyzeClusters --------------------------------------
408 void PndGemFindClustersTB::AnalyzeClusters() {
409  PndGemDigi* digi = NULL;
410  PndGemSensor* sensor = NULL;
411 
412  // cout << "CALLING AnalyzeClusters() in EVENT " << fTNofEvents << ", and have " << fDigiClusters.size() << endl;
413 
414  Bool_t printInfo = kFALSE;
415 
416  // if ( fTNofEvents == 496 ) printInfo = kTRUE;
417 
418  if ( printInfo ) cout << "AnalyzeClusters" << endl;
419 
420  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
421  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
422 
423  // if ( fTNofEvents == 9 && fDigiClusters[idc].detId == 568329088 ) printInfo = kTRUE;
424 
425  if ( printInfo ) {
426  cout << "CLUSTER " << idc << " IN " << fDigiClusters[idc].detId << " AT " << fDigiClusters[idc].cluPos
427  << " FROM " << fDigiClusters[idc].digiNr.size() << " DIGIS ("
428  << fDigiClusters[idc].cluPMn << " TO " << fDigiClusters[idc].cluPMx << "), MAX("
429  << fDigiClusters[idc].cluMPs << ") = " << fDigiClusters[idc].cluMVl << endl;
430  }
431 
432  digi = (PndGemDigi*) fDigis->At(fDigiClusters[idc].digiNr[0]);
433 
434  Int_t stationNr = digi->GetStationNr();
435  Int_t sensorNr = digi->GetSensorNr();
436  Int_t iSide = digi->GetSide();
437 
438  if ( printInfo )
439  cout << "DIGI [" << fDigiClusters[idc].digiNr[0] << "] in station " << stationNr << " sensor " << sensorNr << " side " << iSide << " channel " << digi->GetChannelNr() << endl;
440 
441  sensor = fDigiPar->GetSensor(stationNr, sensorNr);
442 
443  if ( printInfo )
444  cout << "CLUSTERS IN " << fDigiClusters[idc].detId << " COVERS "
445  << sensor->GetDistance(iSide,fDigiClusters[idc].cluPMn,fDigiClusters[idc].cluPMx)+1 << " STRIPS FROM "
446  << fDigiClusters[idc].cluPMn << " TO "
447  << fDigiClusters[idc].cluPMx << " @ "
448  << fDigiClusters[idc].cluTDC << "ns, "
449  << fDigiClusters[idc].digiNr.size() << " STRIPS FIRED." << endl;
450 
451  if ( sensor->GetDistance(iSide,fDigiClusters[idc].cluPMn,fDigiClusters[idc].cluPMx)+1 == fDigiClusters[idc].digiNr.size() ) {
452  if ( printInfo )
453  cout << "The cluster has as many digis as possible strips (" << fDigiClusters[idc].digiNr.size() << ")." << endl;
454  continue;
455  }
456 
457  // PrintCluster(idc);
458 
459  Int_t thisChan = -1;
460  Int_t lastChan = -1;
461  Int_t chanDist = -1;
462  for ( size_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
463  thisChan = fDigiClusters[idc].chanNr[idigi];
464  digi = (PndGemDigi*) fDigis->At(fDigiClusters[idc].digiNr[idigi]);
465 
466  if ( lastChan != -1 )
467  chanDist = sensor->GetDistance(iSide,thisChan,lastChan);
468 
469  if ( lastChan == -1 || chanDist > 1 ) {
470  DigiClusterTB tempDC;
471  tempDC.detId = digi->GetDetectorId();
472  tempDC.digiNr.push_back(fDigiClusters[idc].digiNr[idigi]);
473  tempDC.chanNr.push_back((Int_t)digi->GetChannelNr());
474  tempDC.sigADC.push_back(digi->GetCharge());
475  tempDC.sigTDC.push_back(digi->GetTimeStamp());
476  tempDC.cluTDC = digi->GetTimeStamp();
477  tempDC.cluADC = digi->GetCharge();
478  tempDC.cluPos = digi->GetChannelNr();
479  tempDC.cluPMn = digi->GetChannelNr();
480  tempDC.cluPMx = digi->GetChannelNr();
481  tempDC.cluMPs = digi->GetChannelNr();
482  tempDC.cluMVl = digi->GetCharge();
483  if ( printInfo )
484  cout<< "CREATING ClUSTER " << fDigiClusters.size() << " IN DETECTOR " << tempDC.detId << " / " << digi->GetSide() << " at " << digi->GetChannelNr() << endl;
485  fDigiClusters.push_back(tempDC);
486  }
487  else {
488  AddDigiToCluster(fDigiClusters[idc].digiNr[idigi],fDigiClusters.size()-1);
489  if ( printInfo ) {
490  cout << " ADDED CHANNEL " << fDigiClusters[idc].chanNr[idigi]
491  << " @ " << fDigiClusters[idc].sigTDC[idigi] << "ns"
492  << " TO CLUSTER " << fDigiClusters.size()-1
493  << " TO HAVE NOW " << fDigiClusters[fDigiClusters.size()-1].digiNr.size()
494  << " DIGIS MOVED IT TO: " << fDigiClusters[fDigiClusters.size()-1].cluPos << endl;
495  }
496  }
497  lastChan = thisChan;
498  }
499  fDigiClusters[idc].cluPos = -1.;
500  }
501 
502  // if ( fTNofEvents == 9 ) {
503  // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
504  // if ( fDigiClusters[idc].detId == 568328832 || fDigiClusters[idc].detId == 568329088 ) {
505  // cout << "ANA CLUSTER " << idc << " in " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].cluPMn << " to " << fDigiClusters[idc].cluPMx << " at " << fDigiClusters[idc].cluPos << endl;
506  // }
507  // }
508  // }
509 
510 }
511 // -------------------------------------------------------------------------
512 
513 // ----- Private method WriteClusters --------------------------------------
514 Int_t PndGemFindClustersTB::WriteClusters() {
515 
516  Int_t nClusters = 0;
517  std::vector<Int_t> clusterRefs;
518  PndGemDigi* digi;
519 
520  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
521 
522  // if ( fDigiClusters[idc].cluADC < 1. ) continue;
523  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
524 
525  clusterRefs.clear();
526 
527  for ( size_t id = 0 ; id < fDigiClusters[idc].digiNr.size() ; id++ ) {
528  digi = (PndGemDigi*)fDigis->At(fDigiClusters[idc].digiNr[id]);
529 
530  if ( !digi ) cout << "there is no digi number " << fDigiClusters[idc].digiNr[id] << ", cause there are only " << fDigis->GetEntries() << " of them" << endl;
531  for ( Int_t irf = digi->GetNIndices()-1 ; irf >= 0 ; irf-- ) {
532  Int_t refIndex = digi->GetIndex(irf);
533 
534  Bool_t alreadyKnown = kFALSE;
535  for ( Int_t icr = clusterRefs.size()-1 ; icr >= 0 ; icr-- ) {
536  if ( refIndex == clusterRefs[icr] ) {
537  alreadyKnown = kTRUE;
538  break;
539  }
540  }
541  if ( !alreadyKnown ) clusterRefs.push_back(refIndex);
542  }
543  }
544 
545  // cout << "CREATING CLUSTER " << idc << " IN DETECTOR " << fDigiClusters[idc].detId << " AT " << fDigiClusters[idc].cluPos << " WITH TIME " << fDigiClusters[idc].cluTDC << endl;
546  // PndGemCluster(Int_t iDetectorId, Double_t iChannel, Int_t bChannel, Int_t eChannel, Double_t signal, Double_t time, std::vector<Int_t> index);
547 
548  new ((*fClusters)[nClusters]) PndGemCluster(fDigiClusters[idc].detId,
549  fDigiClusters[idc].cluPos,
550  fDigiClusters[idc].cluPMn,
551  fDigiClusters[idc].cluPMx,
552  fDigiClusters[idc].cluADC,
553  fDigiClusters[idc].cluTDC,
554  clusterRefs);
555 
556  // for ( Int_t icr = clusterRefs.size()-1 ; icr > 0 ; icr-- )
557  // digi->AddIndex(clusterRefs[icr]);
558  if ( fVerbose > 1 )
559  cout << "creating cluster at " << fDigiClusters[idc].cluPos << " with height of " << fDigiClusters[idc].cluADC << " /// compare maximum at " << fDigiClusters[idc].cluMPs << endl;
560  nClusters++;
561  }
562  return nClusters;
563 }
564 // -------------------------------------------------------------------------
565 
566 // ----- Private method CheckClusters -------------------------------
567 void PndGemFindClustersTB::CheckClusters() {
568  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
569  CheckCluster(idc);
570  }
571 }
572 // -------------------------------------------------------------------------
573 
574 // ----- Private method CheckCluster -------------------------------
575 void PndGemFindClustersTB::CheckCluster(Int_t clus) {
576  // cout << "CHECKING CLUSTER " << clus << endl;
577 
578  PndGemDigi* digi = (PndGemDigi*) fDigis->At(fDigiClusters[clus].digiNr[0]);
579 
580  Int_t stationNr = digi->GetStationNr();
581  Int_t sensorNr = digi->GetSensorNr();
582  Int_t iSide = digi->GetSide();
583 
584  PndGemSensor* sensor = fDigiPar->GetSensor(stationNr, sensorNr);
585 
586  Int_t minPart = sensor->GetSensorPart(iSide,fDigiClusters[clus].cluPMn);
587  if ( minPart == -1 ) return;
588 
589  Int_t maxPart = sensor->GetSensorPart(iSide,fDigiClusters[clus].cluPMx);
590  if ( minPart == maxPart ) return;
591 
592  Int_t tempPart = -1;
593  for ( size_t ichan = 1 ; ichan < fDigiClusters[clus].digiNr.size()-1 ; ichan++ ) {
594  tempPart = sensor->GetSensorPart(iSide,fDigiClusters[clus].chanNr[ichan]);
595  if ( tempPart != minPart && tempPart != maxPart ) {
596  // cout << "PARTS ARE DIFFERENT: "
597  // << fDigiClusters[clus].cluPMn << " (" << minPart << ") / "
598  // << fDigiClusters[clus].chanNr[ichan] << " (" << tempPart << ") / "
599  // << fDigiClusters[clus].cluPMx << " (" << maxPart << ")"
600  // << endl;
601  tempPart = 666;
602  break;
603  }
604  }
605 
606  if ( tempPart < 600 ) return;
607 
608  // the temp was 666
609  // it means that the digis are from 3 different parts of the sensor (0,1 and 2),
610  // while clusters can be created of digis from at most 2 parts (0,1 or 0,2).
611  // Have to create TWO clusters consisting of digis from parts (0,1) and (0,2).
612 
613  DigiClusterTB tempDC;
614  fDigiClusters.push_back(tempDC);
615  fDigiClusters.push_back(tempDC);
616 
617  for ( size_t ichan = 0 ; ichan < fDigiClusters[clus].digiNr.size() ; ichan++ ) {
618  tempPart = sensor->GetSensorPart(iSide,fDigiClusters[clus].chanNr[ichan]);
619 
620  if ( tempPart == 0 || tempPart == 1 ) {
621  // cout << "ADDING DIGI " << fDigiClusters[clus].digiNr[ichan]
622  // << " TO CLUSTER " << fDigiClusters.size()-1
623  // << " WITH " << fDigiClusters[fDigiClusters.size()-1].digiNr.size() << endl;
624  AddDigiToCluster(fDigiClusters[clus].digiNr[ichan],fDigiClusters.size()-1);
625  }
626  if ( tempPart == 0 || tempPart == 2 ) {
627  // cout << "ADDING DIGI " << fDigiClusters[clus].digiNr[ichan]
628  // << " TO CLUSTER " << fDigiClusters.size()-2
629  // << " WITH " << fDigiClusters[fDigiClusters.size()-2].digiNr.size() << endl;
630  AddDigiToCluster(fDigiClusters[clus].digiNr[ichan],fDigiClusters.size()-2);
631  }
632  }
633 
634  fDigiClusters[clus].cluPos = -1.;
635 }
636 // -------------------------------------------------------------------------
637 
638 // ----- Private method SortClusters -------------------------------
639 void PndGemFindClustersTB::SortClusters() {
640  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
641  SortCluster(idc);
642  }
643 }
644 // -------------------------------------------------------------------------
645 
646 // ----- Private method SortClusters -------------------------------
647 void PndGemFindClustersTB::SortCluster(Int_t clus) {
648  Int_t channelPosition = 0;
649 
650  for ( Int_t ich = fDigiClusters[clus].cluPMn ; ich <= fDigiClusters[clus].cluPMx ; ich++ ) {
651  for ( size_t idigi = channelPosition ; idigi < fDigiClusters[clus].digiNr.size() ; idigi++ ) {
652  if ( fDigiClusters[clus].chanNr[idigi] == ich ) {
653  Int_t dN = fDigiClusters[clus].digiNr[idigi];
654  Int_t cN = fDigiClusters[clus].chanNr[idigi];
655  Double_t sA = fDigiClusters[clus].sigADC[idigi];
656  Double_t sT = fDigiClusters[clus].sigTDC[idigi];
657  fDigiClusters[clus].digiNr[idigi] = fDigiClusters[clus].digiNr[channelPosition];
658  fDigiClusters[clus].chanNr[idigi] = fDigiClusters[clus].chanNr[channelPosition];
659  fDigiClusters[clus].sigADC[idigi] = fDigiClusters[clus].sigADC[channelPosition];
660  fDigiClusters[clus].sigTDC[idigi] = fDigiClusters[clus].sigTDC[channelPosition];
661  fDigiClusters[clus].digiNr[channelPosition] = dN;
662  fDigiClusters[clus].chanNr[channelPosition] = cN;
663  fDigiClusters[clus].sigADC[channelPosition] = sA;
664  fDigiClusters[clus].sigTDC[channelPosition] = sT;
665  channelPosition++;
666  }
667  }
668  }
669 }
670 // -------------------------------------------------------------------------
671 
672 // ----- Private method PrintClusters -------------------------------
673 void PndGemFindClustersTB::PrintClusters() {
674 
675  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
676  // if ( fDigiClusters[idc].detId != 568329024 ) continue;
677  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
678  PrintCluster(idc);
679  }
680 }
681 // -------------------------------------------------------------------------
682 
683 // ----- Private method PrintCluster -------------------------------
684 void PndGemFindClustersTB::PrintCluster(Int_t clus) {
685 
686  cout << "CLUSTER " << clus << " IN " << fDigiClusters[clus].detId << " AT " << fDigiClusters[clus].cluPos
687  << " CREATED FROM " << fDigiClusters[clus].digiNr.size() << " DIGIS, FROM "
688  << fDigiClusters[clus].cluPMn << " TO " << fDigiClusters[clus].cluPMx << " WITH MAXIMUM OF "
689  << fDigiClusters[clus].cluMVl << " AT " << fDigiClusters[clus].cluMPs << endl;
690 
691  // for ( Int_t idigi = 0 ; idigi < fDigiClusters[clus].digiNr.size() ; idigi++ ) {
692  // cout << idigi << " > "
693  // << fDigiClusters[clus].digiNr[idigi] << " @ "
694  // << fDigiClusters[clus].chanNr[idigi] << " h "
695  // << fDigiClusters[clus].sigADC[idigi] << " t "
696  // << fDigiClusters[clus].sigTDC[idigi] << endl;
697  // }
698 
699  for ( Int_t icol = 0 ; icol < fDigiClusters[clus].cluPMx-fDigiClusters[clus].cluPMn+3 ; icol++ ) cout << "-" << flush; cout<<endl;
700  for ( Double_t fh = 0.99 ; fh >= 0. ; fh-=0.247475 ) {
701  cout << "|" << flush;
702  Int_t thisDigi = 0;
703  Double_t thisADC = 0.;
704  for ( Int_t ichan = fDigiClusters[clus].cluPMn ; ichan <= fDigiClusters[clus].cluPMx ; ichan++ ) {
705  thisADC = 0.;
706  if ( fDigiClusters[clus].chanNr[thisDigi] == ichan ) {
707  thisADC = fDigiClusters[clus].sigADC[thisDigi];
708  thisDigi++;
709  }
710  if ( thisADC > fh*fDigiClusters[clus].cluMVl ) cout << "*" << flush;
711  else cout << " " << flush;
712  }
713  cout << "|" << endl;
714  }
715  for ( Int_t icol = 0 ; icol < fDigiClusters[clus].cluPMx-fDigiClusters[clus].cluPMn+3 ; icol++ ) cout << "-" << flush; cout<<endl;
716 }
717 // -------------------------------------------------------------------------
718 
719 // ----- Private method SetParContainers -------------------------------
720 void PndGemFindClustersTB::SetParContainers() {
721 
722  // Get run and runtime database
723  FairRunAna* run = FairRunAna::Instance();
724  if ( ! run ) Fatal("SetParContainers", "No analysis run");
725 
726  FairRuntimeDb* db = run->GetRuntimeDb();
727  if ( ! db ) Fatal("SetParContainers", "No runtime database");
728 
729  // Get GEM digitisation parameter container
730  fDigiPar = (PndGemDigiPar*) db->getContainer("PndGemDetectors");
731 }
732 // -------------------------------------------------------------------------
733 
734 
735 
736 
737 // ----- Private method Init -------------------------------------------
738 InitStatus PndGemFindClustersTB::Init() {
739 
740  // Get input array
741  FairRootManager* ioman = FairRootManager::Instance();
742  if ( ! ioman ) Fatal("Init", "No FairRootManager");
743  fDigis = (TClonesArray*) ioman->GetObject(fInBranchName.Data());
744 
745  fFunctor = new TimeGap();
746 
747  // Register output array
748  // fClusters = new TClonesArray("PndGemDigi", 1000);
749  // ioman->Register("GEMCluster", "PndGEM", fClusters, kTRUE);
750  fClusters = ioman->Register("GEMCluster", "PndGemCluster", "PndGEM", kTRUE);
751 
752  cout << "-I- " << fName.Data() << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
753  cout << "-I- " << fName.Data() << "::Init(). Initialization succesfull." << endl;
754 
755  PndGemStation* station = (PndGemStation*)fDigiPar->GetStation(0);
756  PndGemSensor* sensor = (PndGemSensor*)station->GetSensor(0);
757  cout << "sensor out rad is " << sensor->GetOuterRadius() << endl;
758 
759 
760  return kSUCCESS;
761 }
762 // -------------------------------------------------------------------------
763 
764 
765 
766 
767 // ----- Private method ReInit -----------------------------------------
768 InitStatus PndGemFindClustersTB::ReInit() {
769 
770  // Create sectorwise digi sets
771  // MakeSets();
772 
773  return kSUCCESS;
774 }
775 // -------------------------------------------------------------------------
776 
777 
778 
779 // ----- Public method Finish ------------------------------------------
780 void PndGemFindClustersTB::Finish() {
781  if ( fClusters ) fClusters->Clear();
782 
783  cout << "-------------------- " << fName.Data() << " : Summary -----------------------" << endl;
784  cout << " Events: " << setw(10) << fTNofEvents << endl;
785  cout << " Digis: " << setw(10) << fTNofDigis << " ( " << (Double_t)fTNofDigis /((Double_t)fTNofEvents) << " per event )" << endl;
786  cout << " Clusters: " << setw(10) << fTNofClusters << " ( " << (Double_t)fTNofClusters/((Double_t)fTNofEvents) << " per event )" << endl;
787  cout << " --> ( " << (Double_t)fTNofDigis /((Double_t)fTNofClusters) << " digis per cluster )" << endl;
788  cout << "---------------------------------------------------------------------" << endl;
789  cout << " >>> CF >>> prep time = " << fPrepTime << "s (get data from input array)" << endl;
790  cout << " >>> CF >>> create time = " << fCreateTime-fPrepTime << "s (create clusters, " << fCreateTime << ")" << endl;
791  cout << " >>> CF >>> sort time = " << fSortTime-fCreateTime << "s (sort clusters, " << fSortTime << ")" << endl;
792  cout << " >>> CF >>> check time = " << fCheckTime-fSortTime << "s (check if geometrically fine, " << fCheckTime << ")" << endl;
793  cout << " >>> CF >>> ana time = " << fAnaTime-fCheckTime << "s (analyze if not splittable, " << fAnaTime << ")" << endl;
794  cout << " >>> CF >>> write time = " << fWriteTime-fAnaTime << "s (wrote clusters to output, " << fWriteTime << ")" << endl;
795  cout << " >>> CF >>> all time = " << fAllTime-fWriteTime << "s (all time spent in Exec, " << fAllTime << ")" << endl;
796  cout << "---------------------------------------------------------------------" << endl;
797 }
798 // -------------------------------------------------------------------------
799 
800 
801 ClassImp(PndGemFindClustersTB)
802 
fVerbose
int fVerbose
Definition: poormantracks.C:24
digi
TClonesArray * digi
Definition: checkdedx_helixhit.C:20
run
Int_t run
Definition: autocutx.C:47
PndGemDigi::GetStationNr
Int_t GetStationNr() const
Definition: PndGemDigi.h:84
map
PndTransMap * map
Definition: sim_emc_apd.C:99
DigiClusterTB::sigTDC
std::vector< Double_t > sigTDC
Definition: PndGemFindClustersTB.h:43
PndGemFindClustersTB::Init
virtual InitStatus Init()
Definition: PndGemFindClustersTB.cxx:738
sensor
TGeoVolume * sensor
Definition: createDummySensor.C:21
PndGemDigiPar
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
PndGemDigi::GetDetectorId
Int_t GetDetectorId() const
Definition: PndGemDigi.h:79
PndGemFindClustersTB::JoinTwoClusters
void JoinTwoClusters(Int_t clus1, Int_t clus2)
Definition: PndGemFindClustersTB.cxx:341
PndGemStation::GetSensor
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
PndGemFindClustersTB::Exec
virtual void Exec(Option_t *opt)
Definition: PndGemFindClustersTB.cxx:134
PndGemFindClustersTB::AddDigiToCluster
void AddDigiToCluster(Int_t digiNr, Int_t clusNr)
Definition: PndGemFindClustersTB.cxx:369
PndGemDigi::GetNIndices
Int_t GetNIndices()
Definition: PndGemDigi.h:102
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndGemDigi::GetSensorNr
Int_t GetSensorNr() const
Definition: PndGemDigi.h:86
refIndex
Int_t refIndex
Definition: checkmomentum.C:30
PndGemDigiPar::GetSensor
PndGemSensor * GetSensor(Int_t stationNr, Int_t sensorNr)
Definition: PndGemDigiPar.cxx:223
DigiClusterTB::digiNr
std::vector< Int_t > digiNr
Definition: PndGemFindClustersTB.h:40
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndGemSensor::GetDistance
Double_t GetDistance(Int_t iSide, Double_t chan1, Double_t chan2)
Definition: PndGemSensor.cxx:494
PndGemFindClustersTB::CompareDigiToClusters
Bool_t CompareDigiToClusters(Int_t digiNumber)
Definition: PndGemFindClustersTB.cxx:264
PndGemDigi::GetChannelNr
Double_t GetChannelNr() const
Definition: PndGemDigi.h:80
PndGemDigiPar::GetNStations
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
name
TString name
Definition: createRootGeoFile.C:82
DigiClusterTB::sigADC
std::vector< Double_t > sigADC
Definition: PndGemFindClustersTB.h:42
PndGemSensor::GetSensorPart
Int_t GetSensorPart(Int_t iSide, Int_t chan)
Definition: PndGemSensor.cxx:405
DigiClusterTB::chanNr
std::vector< Int_t > chanNr
Definition: PndGemFindClustersTB.h:41
PndGemFindClustersTB::fDigiClusters
std::vector< DigiClusterTB > fDigiClusters
Definition: PndGemFindClustersTB.h:102
PndGemDigi::GetIndex
Int_t GetIndex(int i=0) const
Definition: PndGemDigi.h:103
ClassImp
ClassImp(PndAnaContFact)
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
PndGemDigi::GetCharge
Double_t GetCharge() const
Definition: PndGemDigi.h:91
PndGemSensor::GetOuterRadius
Double_t GetOuterRadius() const
Definition: PndGemSensor.h:103
PndGemSensor::GetMeanChannel
Double_t GetMeanChannel(Int_t iSide, Double_t chan1, Double_t weight1, Double_t chan2, Double_t weight2)
Definition: PndGemSensor.cxx:764
PndGemFindClustersTB::ReInit
virtual InitStatus ReInit()
Definition: PndGemFindClustersTB.cxx:768
PndGemDigi::GetSide
Int_t GetSide() const
Definition: PndGemDigi.h:88
PndGemDigiPar::GetStation
PndGemStation * GetStation(Int_t iStation)
Definition: PndGemDigiPar.cxx:176
Generated on Wed Apr 3 2019 10:02:52 for FairRoot/PandaRoot by   doxygen 1.8.5