FairRoot/PandaRoot
PndGemFindClusters.cxx
Go to the documentation of this file.
1 //* $Id: */
2 
3 // -------------------------------------------------------------------------
4 // ----- PndGemFindClusters source file -----
5 // ----- Created 15/02/2009 by R. Karabowicz -----
6 // -------------------------------------------------------------------------
7 
17 #include "PndGemFindClusters.h"
18 
19 #include "PndGemCluster.h"
20 #include "PndGemDigi.h"
21 #include "PndGemDigiPar.h"
22 #include "PndGemSensor.h"
23 #include "PndGemStation.h"
24 
25 #include "PndDetectorList.h"
26 
27 #include "FairRootManager.h"
28 #include "FairRunAna.h"
29 #include "FairRuntimeDb.h"
30 
31 #include "TClonesArray.h"
32 #include "TMath.h"
33 #include "TRandom2.h"
34 
35 #include <iomanip>
36 
37 using std::cout;
38 using std::cerr;
39 using std::endl;
40 using std::flush;
41 using std::fixed;
42 using std::right;
43 using std::left;
44 using std::setw;
45 using std::setprecision;
46 using std::set;
47 using std::map;
48 
49 // ----- Default constructor ------------------------------------------
50 PndGemFindClusters::PndGemFindClusters()
51  : FairTask("GEMFindClusters", 1) {
52  fDigiPar = NULL;
53  fDigis = NULL;
54  fClusters = NULL;
55 
56  fTNofEvents = 0;
57  fTNofDigis = 0;
58  fTNofClusters = 0;
59 }
60 // -------------------------------------------------------------------------
61 
62 
63 
64 // ----- Standard constructor ------------------------------------------
65 PndGemFindClusters::PndGemFindClusters(Int_t iVerbose)
66  : FairTask("GEMFindClusters", iVerbose) {
67  fDigiPar = NULL;
68  fDigis = NULL;
69  fClusters = NULL;
70 
71  fTNofEvents = 0;
72  fTNofDigis = 0;
73  fTNofClusters = 0;
74 }
75 // -------------------------------------------------------------------------
76 
77 
78 
79 // ----- Constructor with name -----------------------------------------
80 PndGemFindClusters::PndGemFindClusters(const char* name, Int_t iVerbose)
81  : FairTask(name, iVerbose) {
82  fDigiPar = NULL;
83  fDigis = NULL;
84  fClusters = NULL;
85 
86  fTNofEvents = 0;
87  fTNofDigis = 0;
88  fTNofClusters = 0;
89 }
90 // -------------------------------------------------------------------------
91 
92 
93 
94 // ----- Destructor ----------------------------------------------------
95 PndGemFindClusters::~PndGemFindClusters() {
96  if ( fClusters ) {
97  fClusters->Delete();
98  delete fClusters;
99  }
100 }
101 // -------------------------------------------------------------------------
102 
103 
104 
105 // ----- Public method Exec --------------------------------------------
106 void PndGemFindClusters::Exec(Option_t*) {
107  if ( fVerbose > 0 ) cout <<"-I- PndGemFindClusters::Exec()"<<endl;
108 
109  fTNofEvents++;
110 
111  fTimer.Start();
112  //Bool_t warn = kFALSE; //[R.K. 01/2017] unused variable?
113 
114  // Clear output array
115  fClusters->Clear();
116 
117  //Int_t nDigis = 0; //[R.K.03/2017] unused variable
118  Int_t nClusters = 0;
119 
120  // Sort digis with respect to time
121  SortDigis();//nDigis = //[R.K.03/2017] unused variable
122 
123  //FindClusters();
124 
125  CreateClusters();
126 
127  if ( fVerbose > 1 ) {
128  Int_t nofC2W = 0;
129  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ )
130  if ( fDigiClusters[idc].cluPos > -0.5 )
131  nofC2W++;
132  cout << ">>>> created " << fDigiClusters.size() << " clusters, " << nofC2W << " to write" << endl;
133  }
134 
135  //Skip Analyze... the join is already done in the FindClusters()..
136  //AnalyzeClusters();
137 
138 // if ( fVerbose ) {
139 // Int_t nofC2W = 0;
140 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ )
141 // if ( fDigiClusters[idc].cluPos > -0.5 )
142 // nofC2W++;
143 // cout << ">>>> analyzed " << fDigiClusters.size() << " clusters, " << nofC2W << " to write" << endl;
144 // }
145 
146  //ClearClusters();
147  //ClearClusters2();//sort and print only..
148  SortClusters();
149 
150  if ( fVerbose > 1 ) {
151  PrintClusters();
152  }
153 
154 // if ( fVerbose ) {
155 // Int_t nofC2W = 0;
156 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ )
157 // if ( fDigiClusters[idc].cluPos > -0.5 )
158 // nofC2W++;
159 // cout << ">>>> cleared " << fDigiClusters.size() << " clusters, " << nofC2W << " to write" << endl;
160 // }
161 
162  nClusters = WriteClusters();
163 
164  fTimer.Stop();
165 
166  // cout << "CREATED " << nClusters << " CLUSTERS OUT OF " << nDigis << " DIGIS IN " << fTimer.RealTime() << "s" << endl;
167 
168  fTNofClusters += nClusters;
169 }
170 // -------------------------------------------------------------------------
171 
172 
173 
174 // ----- Private method SetParContainers -------------------------------
175 void PndGemFindClusters::SetParContainers() {
176 
177  // Get run and runtime database
178  FairRunAna* run = FairRunAna::Instance();
179  if ( ! run ) Fatal("SetParContainers", "No analysis run");
180 
181  FairRuntimeDb* db = run->GetRuntimeDb();
182  if ( ! db ) Fatal("SetParContainers", "No runtime database");
183 
184  // Get GEM digitisation parameter container
185  fDigiPar = (PndGemDigiPar*) db->getContainer("PndGemDetectors");
186 }
187 // -------------------------------------------------------------------------
188 
189 
190 
191 
192 // ----- Private method Init -------------------------------------------
193 InitStatus PndGemFindClusters::Init() {
194 
195  // Get input array
196  FairRootManager* ioman = FairRootManager::Instance();
197  if ( ! ioman ) Fatal("Init", "No FairRootManager");
198  fDigis = (TClonesArray*) ioman->GetObject("GEMDigi");
199 
200  // Register output array
201  //fClusters = new TClonesArray("PndGemDigi", 1000);
202  fClusters = new TClonesArray("PndGemCluster", 1000);
203  ioman->Register("GEMCluster", "Cluster in GEM", fClusters, kTRUE);
204 
205  cout << "-I- " << fName.Data() << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
206  cout << "-I- " << fName.Data() << "::Init(). Initialization succesfull." << endl;
207 
208  PndGemStation* station = (PndGemStation*)fDigiPar->GetStation(0);
209  PndGemSensor* sensor = (PndGemSensor*)station->GetSensor(0);
210  cout << "sensor out rad is " << sensor->GetOuterRadius() << endl;
211 
212  return kSUCCESS;
213 }
214 // -------------------------------------------------------------------------
215 
216 
217 
218 
219 // ----- Private method ReInit -----------------------------------------
220 InitStatus PndGemFindClusters::ReInit() {
221 
222  // Create sectorwise digi sets
223  // MakeSets();
224 
225  return kSUCCESS;
226 }
227 // -------------------------------------------------------------------------
228 
229 
230 // ----- Private method SortDigis --------------------------------------
231 Int_t PndGemFindClusters::SortDigis() {
232  if ( fVerbose ) cout << "-I- PndGemFindClusters::SortDigis()" <<endl;
233 
234  // Check input array
235  if ( ! fDigis ) {
236  cout << "-E- " << fName << "::SortDigis: No input array!" << endl;
237  return -1;
238  }
239  TRandom2* randNumber = new TRandom2();
240 
241  //PndGemDigi* digi; //[R.K. 01/2017] unused variable?
242  Int_t nDigis = fDigis->GetEntriesFast();
243 
244  fTimeOrderedDigis.clear();
245 
246  fTNofDigis += nDigis;
247 
248  for (Int_t iDigi=0; iDigi<nDigis; iDigi++) {
249  //digi = (PndGemDigi*) fDigis->At(iDigi); //[R.K. 01/2017] unused variable?
250  // Double_t time = digi->GetTimeStamp();
251  Double_t time = randNumber->Gaus(0.,5.);
252 
253  fTimeOrderedDigis[time] = iDigi;
254  }
255 
256  // cout << "have " << nDigis << " digis ( " << fTimeOrderedDigis.size() << " ) " << endl;
257 
258  map<Double_t, Int_t>::iterator mapIt;
259  Int_t nd = 0;
260  for (mapIt=fTimeOrderedDigis.begin(); mapIt!=fTimeOrderedDigis.end(); mapIt++) {
261  // cout << mapIt->first << " --> " << mapIt->second << endl;
262  nd++;
263  }
264  // cout << " ---> " << nd << endl;
265 
266  return nDigis;
267 }
268 // -------------------------------------------------------------------------
269 
270 
271 // // ----- Private method FindClusters --------------------------------------
272 // //NOT USED..
273 // void PndGemFindClusters::FindClusters() {
274 // if ( fVerbose ) {
275 // cout << "-I- PndGemFindClusters::FindClusters()" <<endl;
276 // }
277 
278 // PndGemDigi* digi;
279 
280 // map<Double_t, Int_t>::iterator mapIt;
281 
282 // fDigiClusters.clear();
283 
284 // for (mapIt=fTimeOrderedDigis.begin(); mapIt!=fTimeOrderedDigis.end(); mapIt++) {
285 
286 // if ( fVerbose ) {
287 // cout << fTNofEvents << ". fighting for " << mapIt->second << endl;
288 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
289 // cout << " cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis:" << endl;// (from " << lowChannel << " to " << topChannel << "):" << endl;
290 // }
291 // }
292 
293 // //Bool_t digiUsed = CompareDigiToClusters(mapIt->second);
294 // Bool_t digiUsed = CompareDigiToClustersDigis(mapIt->second);
295 
296 // if ( digiUsed ) continue;
297 
298 // digi = (PndGemDigi*) fDigis->At(mapIt->second);
299 
300 // if ( fVerbose > 1 )
301 // cout << "creating cluster " << fDigiClusters.size() << " from digi " << mapIt->second << " : " << digi->GetDetectorId() << " / " << digi->GetChannelNr() << endl;
302 // DigiCluster tempDC;
303 // tempDC.detId = digi->GetDetectorId();
304 // tempDC.digiNr.push_back(mapIt->second);
305 // tempDC.chanNr.push_back((Int_t)digi->GetChannelNr());
306 // tempDC.sigADC.push_back(digi->GetCharge());
307 // tempDC.cluTDC = digi->GetTimeStamp();
308 // tempDC.cluADC = digi->GetCharge();
309 // tempDC.cluPos = digi->GetChannelNr();
310 // tempDC.cluPMn = digi->GetChannelNr();
311 // tempDC.cluPMx = digi->GetChannelNr();
312 // tempDC.cluMPs = digi->GetChannelNr();
313 // tempDC.cluMVl = digi->GetCharge();
314 
315 // fDigiClusters.push_back(tempDC);
316 // }
317 
318 // }
319 // // -------------------------------------------------------------------------
320 
321 // -------------------------------------------------------------------------
322 Int_t PndGemFindClusters::CreateClusters() {
323  if ( fVerbose > 0 ) {
324  cout << "-I- PndGemFindClusters::CreateClusters()" <<endl;
325  }
326 
327  fDigiClusters.clear();
328 
329  PndGemDigi* digi = NULL;
330  for ( Int_t idigi = 0 ; idigi < fDigis->GetEntries() ; idigi++ ) {
331  digi = (PndGemDigi*) fDigis->At(idigi);
332 
333  if ( fVerbose ) {
334  // cout << "digi " << idigi << " in " << digi->GetDetectorId() << ", there are already " << fDigiClusters.size() << " clusters" << endl;
335  // cout << idigi << " ---> " << digi->GetDetectorId() << " @ " << digi->GetChannelNr() << " @ " << digi->GetTimeStamp() << endl;
336  }
337  if ( CompareDigiToClustersDigis(idigi) == kTRUE ) continue;
338 
339  DigiCluster tempDC;
340  tempDC.detId = digi->GetDetectorId();
341  tempDC.digiNr.push_back(idigi);
342  tempDC.chanNr.push_back((Int_t)digi->GetChannelNr());
343  tempDC.sigADC.push_back(digi->GetCharge());
344  tempDC.cluTDC = digi->GetTimeStamp();
345  tempDC.cluADC = digi->GetCharge();
346  tempDC.cluPos = digi->GetChannelNr();
347  tempDC.cluPMn = digi->GetChannelNr();
348  tempDC.cluPMx = digi->GetChannelNr();
349  tempDC.cluMPs = digi->GetChannelNr();
350  tempDC.cluMVl = digi->GetCharge();
351  fDigiClusters.push_back(tempDC);
352 
353  }
354  return fDigiClusters.size();
355 }
356 // -------------------------------------------------------------------------
357 
358 // // ----- Private method AnalyzeClusters --------------------------------------
359 // //CURRENTLY NOT USED...
360 // void PndGemFindClusters::AnalyzeClusters() {
361 // if ( fVerbose > 0 ) {
362 // cout << "STARTING CLUSTER ANALYSIS WITH " << fDigiClusters.size() << " CLUSTERS" << endl;
363 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
364 // cout << idc << ".th cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis:" << endl;// (from " << lowChannel << " to " << topChannel << "):" << endl;
365 // }
366 // }
367 
368 // Int_t nClusters = fDigiClusters.size();
369 
370 // PndGemSensor* sensor;
371 // Int_t side;
372 // PndGemDigi* digi;
373 // Double_t dist;
374 
375 // std::vector<Bool_t> digiVect(20000,kFALSE);
376 // std::vector<Bool_t> commVect(20000,kFALSE);
377 
378 // for ( Int_t idc1 = 0 ; idc1 < nClusters ; idc1++ ) {
379 
380 // digiVect.assign(20000,kFALSE);
381 
382 // Bool_t jointDigi = kFALSE;
383 
384 // for ( Int_t idc2 = idc1+1 ; idc2 < nClusters ; idc2++ ) {
385 // if ( fDigiClusters[idc1].detId != fDigiClusters[idc2].detId ) continue;
386 // digi = (PndGemDigi*)fDigis->At(fDigiClusters[idc1].digiNr[0]);
387 // sensor = fDigiPar->GetSensor(digi->GetStationNr(),digi->GetSensorNr());
388 // side = digi->GetSide();
389 
390 // dist = sensor->GetDistance(side,fDigiClusters[idc1].cluPos,fDigiClusters[idc2].cluPos);
391 
392 // if ( dist < 0 ) continue;
393 // if ( dist > 30 ) continue;
394 
395 // // cout << "WILL JOIN cluster " << idc1 << " and " << idc2 << endl;
396 
397 // commVect.assign(20000,kFALSE);
398 // for ( Int_t id = 0 ; id < fDigiClusters[idc2].digiNr.size() ; id++ ) {
399 // // cout << "cluster2 id = " << id << " and in fact " << fDigiClusters[idc2].digiNr[id] << endl;
400 // commVect[fDigiClusters[idc2].digiNr[id]] = kTRUE;
401 // }
402 
403 // if ( fVerbose > 1 )
404 // cout << "CLUSTER " << idc2 << " (" << fDigiClusters[idc2].digiNr.size() << " digis at " << fDigiClusters[idc2].cluPos << ") joined with cluster " << idc1 << " ---> " << flush;
405 // for ( Int_t id = 0 ; id < fDigiClusters[idc1].digiNr.size() ; id++ ) {
406 // if ( commVect[fDigiClusters[idc1].digiNr[id]] == kFALSE ) {
407 // // cout << "adding digi id " << id << " and in fact " << fDigiClusters[idc1].digiNr[id] << endl;
408 // digi = (PndGemDigi*)fDigis->At(fDigiClusters[idc1].digiNr[id]);
409 // fDigiClusters[idc2].digiNr.push_back(fDigiClusters[idc1].digiNr[id]);
410 // fDigiClusters[idc2].chanNr.push_back((Int_t)digi->GetChannelNr());
411 // fDigiClusters[idc2].sigADC.push_back(digi->GetCharge());
412 // fDigiClusters[idc2].cluTDC = digi->GetTimeStamp();
413 // fDigiClusters[idc2].cluPos = sensor->GetMeanChannel(side,fDigiClusters[idc2].cluPos,fDigiClusters[idc2].cluADC,digi->GetChannelNr(),digi->GetCharge());
414 // if ( digi->GetCharge() > fDigiClusters[idc2].cluMVl ) {
415 // fDigiClusters[idc2].cluMVl = digi->GetCharge();
416 // fDigiClusters[idc2].cluMPs = digi->GetChannelNr();
417 // }
418 
419 // fDigiClusters[idc2].cluADC = fDigiClusters[idc2].cluADC+digi->GetCharge();
420 // }
421 // }
422 // if ( fVerbose > 1 )
423 // cout << "(" << fDigiClusters[idc2].digiNr.size() << " digis at " << fDigiClusters[idc2].cluPos << ")." << endl;
424 
425 // jointDigi = kTRUE;
426 // if ( fVerbose > 1 )
427 // cout << "comparing cluster at " << fDigiClusters[idc1].cluPos << " with " << fDigiClusters[idc2].cluPos << " gives dist of " << dist << endl;
428 // }
429 
430 // if ( jointDigi == kTRUE )
431 // fDigiClusters[idc1].cluPos = -1.;
432 // }
433 
434 // if ( fVerbose ) {
435 // cout << "FINISHING CLUSTER ANALYSIS WITH " << fDigiClusters.size() << " CLUSTERS" << endl;
436 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
437 // cout << idc << ".th cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis:" << endl;// (from " << lowChannel << " to " << topChannel << "):" << endl;
438 // }
439 
440 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
441 // // if ( fDigiClusters[idc].cluADC < 105 ) continue;
442 // cout << "cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis:" << endl;// (from " << lowChannel << " to " << topChannel << "):" << endl;
443 // for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
444 // cout << idigi << " > "
445 // << fDigiClusters[idc].digiNr[idigi] << " @ "
446 // << fDigiClusters[idc].chanNr[idigi] << " h "
447 // << fDigiClusters[idc].sigADC[idigi] << endl;
448 // }
449 // }
450 // }
451 
452 // }
453 // // -------------------------------------------------------------------------
454 
455 // // ----- Private method ClearClusters --------------------------------------
456 // //CURRENTLY NOT USED..
457 // void PndGemFindClusters::ClearClusters() {
458 // if ( fVerbose )
459 // cout << "CLEAR CLUSTERS STARTING FROM " << fDigiClusters.size() << endl;
460 
461 // Int_t intEvent = -1;
462 
463 // PndGemDigi* digi;
464 // PndGemSensor* sensor;
465 
466 // Int_t nofClustersAtBegin = fDigiClusters.size();
467 // for ( Int_t idc = 0 ; idc < nofClustersAtBegin ; idc++ ) {
468 // if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
469 
470 // if ( fVerbose || fTNofEvents == intEvent )
471 // cout << "-------------------------------------------------" << endl;
472 
473 // // GET THE LOWEST AND HIGHEST CHANNEL NUMBER
474 // Int_t lowChannel = 1000000;
475 // Int_t topChannel = -1;
476 // for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
477 // if ( fVerbose ) {
478 // cout << idigi << " > "
479 // << fDigiClusters[idc].digiNr[idigi] << " @ "
480 // << fDigiClusters[idc].chanNr[idigi] << " h "
481 // << fDigiClusters[idc].sigADC[idigi] << endl;
482 // }
483 // if ( fDigiClusters[idc].chanNr[idigi] < lowChannel ) lowChannel = fDigiClusters[idc].chanNr[idigi];
484 // if ( fDigiClusters[idc].chanNr[idigi] > topChannel ) topChannel = fDigiClusters[idc].chanNr[idigi];
485 // }
486 // // ---------------------------------------------------------
487 
488 // // MAKE ORDER IN THE CLUSTER
489 // Int_t channelPosition = 0;
490 // for ( Int_t ich = lowChannel ; ich <= topChannel ; ich++ ) {
491 // for ( Int_t idigi = channelPosition ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
492 // if ( fDigiClusters[idc].chanNr[idigi] == ich ) {
493 // Int_t dN = fDigiClusters[idc].digiNr[idigi];
494 // Int_t cN = fDigiClusters[idc].chanNr[idigi];
495 // Double_t sA = fDigiClusters[idc].sigADC[idigi];
496 // fDigiClusters[idc].digiNr[idigi] = fDigiClusters[idc].digiNr[channelPosition];
497 // fDigiClusters[idc].chanNr[idigi] = fDigiClusters[idc].chanNr[channelPosition];
498 // fDigiClusters[idc].sigADC[idigi] = fDigiClusters[idc].sigADC[channelPosition];
499 // fDigiClusters[idc].digiNr[channelPosition] = dN;
500 // fDigiClusters[idc].chanNr[channelPosition] = cN;
501 // fDigiClusters[idc].sigADC[channelPosition] = sA;
502 // channelPosition++;
503 // }
504 // }
505 // }
506 // // ---------------------------------------------------------
507 
508 // // PRINT
509 // if ( fVerbose || fTNofEvents == intEvent ) {
510 // cout << fTNofEvents << " cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis (from " << lowChannel << " to " << topChannel << "):" << endl;
511 // for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
512 // cout << idigi << " > "
513 // << fDigiClusters[idc].digiNr[idigi] << " @ "
514 // << fDigiClusters[idc].chanNr[idigi] << " h "
515 // << fDigiClusters[idc].sigADC[idigi] << endl;
516 // }
517 // }
518 // // ---------------------------------------------------------
519 
520 // // FIND PEAKS AND HOLES BETWEEN TEHM
521 // Bool_t showDets = kFALSE;
522 // if ( fVerbose ) {
523 // showDets = kTRUE;
524 // cout << fTNofEvents << " cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis (from " << lowChannel << " to " << topChannel << "):" << endl;
525 // /*for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
526 // cout << idigi << " > "
527 // << fDigiClusters[idc].digiNr[idigi] << " @ "
528 // << fDigiClusters[idc].chanNr[idigi] << " h "
529 // << fDigiClusters[idc].sigADC[idigi] << endl;
530 // }*/
531 // }
532 // Int_t lastPos = fDigiClusters[idc].chanNr[0];
533 // Double_t lastADC = fDigiClusters[idc].sigADC[0];
534 // Double_t maxADC = fDigiClusters[idc].sigADC[0];
535 // Double_t gotPeak = 40.;
536 // Double_t newPeak = 0;
537 // Double_t minADC = 1000.;
538 // Int_t minPos = -1;
539 // Int_t nofPeaks = 0;
540 // std::vector<Int_t> clusterBorders;
541 // if ( showDets ) cout << "0 > "
542 // << fDigiClusters[idc].digiNr[0] << " @ "
543 // << fDigiClusters[idc].chanNr[0] << " h "
544 // << fDigiClusters[idc].sigADC[0] << endl;
545 // // for ( Int_t idigi = 1 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
546 // Int_t thisDigi = 1;
547 // Double_t thisADC = 0.;
548 // for ( Int_t ichan = lowChannel+1 ; ichan <= topChannel ; ichan++ ) {
549 // thisADC = 0.;
550 // if ( fDigiClusters[idc].chanNr[thisDigi] == ichan ) {
551 // thisADC = fDigiClusters[idc].sigADC[thisDigi];
552 // thisDigi++;
553 // }
554 // // if ( fDigiClusters[idc].chanNr[idigi] > lastPos+2 )
555 // // cout << "will probably make a break here between " << idigi-1 << " and " << idigi << endl;
556 // // if ( gotPeak < 0. &&
557 // if ( thisADC < minADC ) { // this is min
558 // minADC = thisADC;
559 // minPos = ichan;
560 // // cout << "minadc " << minADC << " at " << minPos << endl;
561 // }
562 // if ( thisADC < maxADC*.5 && gotPeak < 0. ) { // this is max
563 // gotPeak = 40.;
564 // if ( showDets ) cout << "is less than half of last max" << endl;
565 // }
566 // if ( gotPeak > 0. ) { // got peak already, adc should be going down
567 // if ( thisADC >= lastADC ) { // but adc not falling down
568 // if ( thisADC > lastADC ) newPeak+=1.0; // possibly a new peak
569 // if ( thisADC == lastADC ) {
570 // if ( thisADC != 0. ) newPeak+=0.5;
571 // if ( thisADC == 0. ) newPeak+=1.0;
572 // }
573 // if ( newPeak >= 2 ) { // it is a new peak
574 // nofPeaks++;
575 // if ( showDets || fTNofEvents == intEvent )
576 // cout << "NEW PEAK starts from " << minPos << " ( " << minADC << " )." << endl;
577 // clusterBorders.push_back(minPos);
578 // minADC = 1000.;
579 // minPos = -1;
580 // gotPeak = -1.;
581 // maxADC = 0.;
582 // }
583 // }
584 // else if ( thisADC < lastADC ) {
585 // newPeak = 0;
586 // }
587 // }
588 // lastPos = ichan;
589 // lastADC = thisADC;
590 // if ( showDets )
591 // cout << thisDigi << " > " << gotPeak << " "
592 // << fDigiClusters[idc].digiNr[thisDigi] << " @ "
593 // << fDigiClusters[idc].chanNr[thisDigi] << " h "
594 // << fDigiClusters[idc].sigADC[thisDigi] << endl;
595 // if ( lastADC > maxADC ) {
596 // maxADC = lastADC;
597 // if ( newPeak == 0 ) {
598 // minPos = -1;
599 // minADC = 1000.;
600 // }
601 // }
602 // }
603 
604 // if ( showDets ) {
605 // cout << "NOF PEAKS = " << nofPeaks << " : " << flush;
606 // for ( Int_t ipk = 0 ; ipk < nofPeaks ; ipk++ )
607 // cout << "New" << flush;
608 // cout << endl;
609 // for ( Int_t ipk = 0 ; ipk < nofPeaks ; ipk++ )
610 // cout << "---> " << clusterBorders[ipk] << endl;
611 // }
612 // // ---------------------------------------------------------
613 
614 // // ----> still do not consider strips with 0 charge...
615 
616 // if ( nofPeaks < 2 ) continue; // do not have to do anything...
617 // // DIVIDE CLUSTER IF nofPeaks > 1
618 // digi = (PndGemDigi*) fDigis->At(fDigiClusters[idc].digiNr[0]);
619 
620 // Int_t stationNr = digi->GetStationNr();
621 // Int_t sensorNr = digi->GetSensorNr();
622 // Int_t iSide = digi->GetSide();
623 
624 // sensor = fDigiPar->GetSensor(stationNr, sensorNr);
625 
626 // thisDigi = 0;
627 // Int_t ichan = lowChannel;
628 
629 // if ( fTNofEvents == intEvent ) {
630 // cout << fTNofEvents << " - - - - CLUSTER " << idc << " FROM " << lowChannel << " TO " << topChannel << " SHOULD BE DIVIDED IN " << nofPeaks << " PEAKS" << endl;
631 // for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
632 // cout << " > > > > > > > " << idigi << " > "
633 // << fDigiClusters[idc].digiNr[idigi] << " @ "
634 // << fDigiClusters[idc].chanNr[idigi] << " h "
635 // << fDigiClusters[idc].sigADC[idigi] << endl;
636 // }
637 // }
638 
639 // for ( Int_t inp = 0 ; inp < nofPeaks ; inp++ ) {
640 // if ( fTNofEvents == intEvent ) {
641 // cout << "*** p e a k " << inp << " to " << clusterBorders[inp+1] << endl;
642 // }
643 
644 // DigiCluster tempDC;
645 // tempDC.detId = fDigiClusters[idc].detId;
646 // tempDC.cluTDC = fDigiClusters[idc].cluTDC;
647 // tempDC.cluADC = 0.;
648 // tempDC.cluPos = 0.;
649 // tempDC.cluMPs = 0.;
650 // tempDC.cluMVl = 0.;
651 
652 // /* for ( ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
653 // if ( fDigiClusters[idc].chanNr[idigi] > clusterBorders[inp+1] ) break;
654 // tempDC.digiNr.push_back(fDigiClusters[idc].digiNr[idigi]);
655 // tempDC.chanNr.push_back(fDigiClusters[idc].chanNr[idigi]);
656 // tempDC.sigADC.push_back(fDigiClusters[idc].sigADC[idigi]);
657 // if ( fDigiClusters[idc].sigADC[idigi] > tempDC.cluMVl ) {
658 // tempDC.cluMVl = fDigiClusters[idc].sigADC[idigi];
659 // tempDC.cluMPs = fDigiClusters[idc].chanNr[idigi];
660 // }
661 // if ( tempDC.digiNr.size() == 1 ) {
662 // tempDC.cluPos = fDigiClusters[idc].chanNr[idigi];
663 // tempDC.cluADC = fDigiClusters[idc].sigADC[idigi];
664 // }
665 // else {
666 // tempDC.cluPos = sensor->GetMeanChannel(iSide,tempDC.cluPos,tempDC.cluADC,fDigiClusters[idc].chanNr[idigi],fDigiClusters[idc].sigADC[idigi]);
667 // tempDC.cluADC = tempDC.cluADC+fDigiClusters[idc].sigADC[idigi];
668 // }
669 // }*/
670 // for ( ; ichan <= topChannel ; ichan++ ) {
671 // if ( fTNofEvents == intEvent ) {
672 // cout << " checking channel " << ichan << endl;
673 // }
674 
675 // if ( inp < nofPeaks-1 && ichan > clusterBorders[inp+1] ) break;
676 
677 // if ( fDigiClusters[idc].chanNr[thisDigi] != ichan )
678 // continue;
679 
680 // // if ( fDigiClusters[idc].digiNr[thisDigi] > 3000 )
681 // if ( fTNofEvents == intEvent ) {
682 // cout << fTNofEvents << " > putting digiNr = " << fDigiClusters[idc].digiNr[thisDigi]
683 // << " (" << fDigiClusters[idc].chanNr[thisDigi] << ") "
684 // << " from " << thisDigi << " coming from split of cluster " << idc << endl;
685 // }
686 
687 // tempDC.digiNr.push_back(fDigiClusters[idc].digiNr[thisDigi]);
688 // tempDC.chanNr.push_back(fDigiClusters[idc].chanNr[thisDigi]);
689 // tempDC.sigADC.push_back(fDigiClusters[idc].sigADC[thisDigi]);
690 // if ( fDigiClusters[idc].sigADC[thisDigi] > tempDC.cluMVl ) {
691 // tempDC.cluMVl = fDigiClusters[idc].sigADC[thisDigi];
692 // tempDC.cluMPs = fDigiClusters[idc].chanNr[thisDigi];
693 // }
694 // if ( tempDC.digiNr.size() == 1 ) {
695 // tempDC.cluPos = fDigiClusters[idc].chanNr[thisDigi];
696 // tempDC.cluADC = fDigiClusters[idc].sigADC[thisDigi];
697 // }
698 // else {
699 // tempDC.cluPos = sensor->GetMeanChannel(iSide,tempDC.cluPos,tempDC.cluADC,fDigiClusters[idc].chanNr[thisDigi],fDigiClusters[idc].sigADC[thisDigi]);
700 // tempDC.cluADC = tempDC.cluADC+fDigiClusters[idc].sigADC[thisDigi];
701 // }
702 // thisDigi++;
703 // }
704 // if ( showDets )
705 // cout << "just to confirm, added cluster with ADC = " << tempDC.cluADC << " at " << tempDC.cluPos
706 // << ", with MAX = " << tempDC.cluMVl << " at " << tempDC.cluMPs
707 // << " with " << tempDC.digiNr.size() << " digis" << endl;
708 
709 // fDigiClusters.push_back(tempDC);
710 // }
711 // fDigiClusters[idc].cluPos = -1.;
712 
713 // // ---------------------------------------------------------
714 // }
715 // }
716 // // -------------------------------------------------------------------------
717 
718 
719 // // ----- Private method --------------------------------------
720 // // It does only sorting and print
721 // // no double peak analysis
722 // void PndGemFindClusters::ClearClusters2() {
723 // if ( fVerbose )
724 // cout << "CLEAR CLUSTERS STARTING FROM " << fDigiClusters.size() << endl;
725 
726 // Int_t intEvent = -1;
727 
728 // PndGemDigi* digi;
729 // PndGemSensor* sensor;
730 
731 // Int_t nofClustersAtBegin = fDigiClusters.size();
732 // for ( Int_t idc = 0 ; idc < nofClustersAtBegin ; idc++ ) {
733 // if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
734 
735 // if ( fVerbose || fTNofEvents == intEvent )
736 // cout << "-------------------------------------------------" << endl;
737 
738 // // GET THE LOWEST AND HIGHEST CHANNEL NUMBER
739 // Int_t lowChannel = 1000000;
740 // Int_t topChannel = -1;
741 // for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
742 // if ( fVerbose ) {
743 // cout << idigi << " > "
744 // << fDigiClusters[idc].digiNr[idigi] << " @ "
745 // << fDigiClusters[idc].chanNr[idigi] << " h "
746 // << fDigiClusters[idc].sigADC[idigi] << endl;
747 // }
748 // if ( fDigiClusters[idc].chanNr[idigi] < lowChannel ) lowChannel = fDigiClusters[idc].chanNr[idigi];
749 // if ( fDigiClusters[idc].chanNr[idigi] > topChannel ) topChannel = fDigiClusters[idc].chanNr[idigi];
750 // }
751 // // ---------------------------------------------------------
752 
753 // // MAKE ORDER IN THE CLUSTER
754 // Int_t channelPosition = 0;
755 // for ( Int_t ich = lowChannel ; ich <= topChannel ; ich++ ) {
756 // for ( Int_t idigi = channelPosition ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
757 // if ( fDigiClusters[idc].chanNr[idigi] == ich ) {
758 // Int_t dN = fDigiClusters[idc].digiNr[idigi];
759 // Int_t cN = fDigiClusters[idc].chanNr[idigi];
760 // Double_t sA = fDigiClusters[idc].sigADC[idigi];
761 // fDigiClusters[idc].digiNr[idigi] = fDigiClusters[idc].digiNr[channelPosition];
762 // fDigiClusters[idc].chanNr[idigi] = fDigiClusters[idc].chanNr[channelPosition];
763 // fDigiClusters[idc].sigADC[idigi] = fDigiClusters[idc].sigADC[channelPosition];
764 // fDigiClusters[idc].digiNr[channelPosition] = dN;
765 // fDigiClusters[idc].chanNr[channelPosition] = cN;
766 // fDigiClusters[idc].sigADC[channelPosition] = sA;
767 // channelPosition++;
768 // }
769 // }
770 // }
771 // // ---------------------------------------------------------
772 
773 // // PRINT
774 // if ( fVerbose || fTNofEvents == intEvent ) {
775 // cout << fTNofEvents << " cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis (from " << lowChannel << " to " << topChannel << "):" << endl;
776 // for ( Int_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
777 // cout << idigi << " > "
778 // << fDigiClusters[idc].digiNr[idigi] << " @ "
779 // << fDigiClusters[idc].chanNr[idigi] << " h "
780 // << fDigiClusters[idc].sigADC[idigi] << endl;
781 // }
782 // }
783 
784 // }
785 
786 // }
787 // // -------------------------------------------------------------------------
788 
789 // ----- Private method WriteClusters --------------------------------------
790 Int_t PndGemFindClusters::WriteClusters() {
791  if ( fVerbose > 0 ) cout << "-I- PndGemFindClusters::WriteClusters()" <<endl;
792 
793  Int_t nClusters = 0;
794  //PndGemDigi* digi; //[R.K. 01/2017] unused variable?
795 
796  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
797 
798  // if ( fDigiClusters[idc].cluADC < 1. ) continue;
799  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
800 
801  if ( fVerbose > 1 ) {
802  cout << "cluster at " << fDigiClusters[idc].cluPos << ", h = " << fDigiClusters[idc].cluADC << ", for detId = " << fDigiClusters[idc].detId << " from " << fDigiClusters[idc].digiNr.size() << " digis (from " << fDigiClusters[idc].cluPMn << " to " << fDigiClusters[idc].cluPMx << "): " << endl;
803 
804  for ( size_t idigi = 0 ; idigi < fDigiClusters[idc].digiNr.size() ; idigi++ ) {
805  cout <<" digiNr=" << fDigiClusters[idc].digiNr[idigi]
806  <<" chanNr=" << fDigiClusters[idc].chanNr[idigi]
807  <<" sigADC=" << fDigiClusters[idc].sigADC[idigi] << endl;
808  }
809  }
810 
811  new ((*fClusters)[nClusters]) PndGemCluster(fDigiClusters[idc].detId,
812  fDigiClusters[idc].cluPos,
813  fDigiClusters[idc].cluPMn,
814  fDigiClusters[idc].cluPMx,
815  fDigiClusters[idc].cluADC,
816  fDigiClusters[idc].cluTDC,
817  fDigiClusters[idc].digiNr);
818 
819  if ( fVerbose > 1 )
820  cout << "creating cluster at " << fDigiClusters[idc].cluPos << " with height of " << fDigiClusters[idc].cluADC << " /// compare maximum at " << fDigiClusters[idc].cluMPs << endl;
821  nClusters++;
822  }
823  return nClusters;
824 }
825 // -------------------------------------------------------------------------
826 
827 
828 // // ----- Private method ------------------------------------------
829 // //NOT USED.. keep old structure..
830 // Bool_t PndGemFindClusters::CompareDigiToClusters(Int_t digiNumber) {
831 // PndGemSensor* sensor;
832 // PndGemDigi* digi;
833 
834 // Bool_t digiUsed = kFALSE;
835 // for ( Int_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
836 // digi = (PndGemDigi*) fDigis->At(digiNumber);
837 
838 // if ( digi->GetDetectorId() != fDigiClusters[idc].detId ) continue;
839 
840 // Int_t stationNr = digi->GetStationNr();
841 // Int_t sensorNr = digi->GetSensorNr();
842 // Int_t iSide = digi->GetSide();
843 
844 // sensor = fDigiPar->GetSensor(stationNr, sensorNr);
845 
846 // Double_t digiClustDist = sensor->GetDistance(iSide,digi->GetChannelNr(),fDigiClusters[idc].cluPos);
847 // if ( digiClustDist < 0. ) continue; // there's no connection
848 // if ( digiClustDist > 50. ) continue; // they are too far
849 
850 // if ( fVerbose > 1 )
851 // cout << "digi " << digiNumber << " ( " << digi->GetDetectorId() << " / " << digi->GetChannelNr() << " ) could attach to cluster "
852 // << idc << " ( " << fDigiClusters[idc].cluPos << " ) with distance " << digiClustDist << " . " << fDigiClusters[idc].digiNr.size() << " digis in cluster" << flush;
853 
854 // fDigiClusters[idc].digiNr.push_back(digiNumber);
855 // fDigiClusters[idc].chanNr.push_back((Int_t)digi->GetChannelNr());
856 // fDigiClusters[idc].sigADC.push_back(digi->GetCharge());
857 // fDigiClusters[idc].cluTDC = digi->GetTimeStamp();
858 // fDigiClusters[idc].cluPos = sensor->GetMeanChannel(iSide,fDigiClusters[idc].cluPos,fDigiClusters[idc].cluADC,digi->GetChannelNr(),digi->GetCharge());
859 // fDigiClusters[idc].cluADC = fDigiClusters[idc].cluADC+digi->GetCharge();
860 // if ( digi->GetCharge() > fDigiClusters[idc].cluMVl ) {
861 // fDigiClusters[idc].cluMVl = digi->GetCharge();
862 // fDigiClusters[idc].cluMPs = digi->GetChannelNr();
863 // }
864 
865 // if ( fVerbose > 1 )
866 // cout << " >>> new pos = " << fDigiClusters[idc].cluPos << endl;
867 
868 // digiUsed = kTRUE;
869 // }
870 // return digiUsed;
871 
872 // }
873 // // -------------------------------------------------------------------------
874 
875 // ----- Private method ------------------------------------------
876 //compares the sample digi and all digis in clusters
877 Bool_t PndGemFindClusters::CompareDigiToClustersDigis(Int_t digiNumber) {
878  if ( fVerbose > 0 ) cout << "PndGemFindClusters::CompareDigiToClustersDigis" << endl;
879 
880  PndGemSensor* sensor;
881  PndGemDigi* digi;
882  digi = (PndGemDigi*) fDigis->At(digiNumber);
883 
884  Int_t iDigiUsed = -1;
885  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
886  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
887  if ( digi->GetDetectorId() != fDigiClusters[idc].detId ) continue;
888 
889  Int_t stationNr = digi->GetStationNr();
890  Int_t sensorNr = digi->GetSensorNr();
891  Int_t iSide = digi->GetSide();
892 
893  sensor = fDigiPar->GetSensor(stationNr, sensorNr);
894 
895  Double_t digiDigiDist=999999999, testDist;
896  for ( size_t id = 0 ; id < fDigiClusters[idc].digiNr.size() ; id++ ) {
897  testDist = sensor->GetDistance2(iSide,digi->GetChannelNr(),fDigiClusters[idc].chanNr[id]);
898  if ( testDist < digiDigiDist ) digiDigiDist=testDist;
899  }
900 
901  if ( digiDigiDist < 0. ) continue;
902  if ( digiDigiDist > 1 ) continue;//1 for next strip to accept
903 
904  if ( fVerbose > 1 )
905  cout << "digi " << digiNumber << " ( " << digi->GetDetectorId() << " / " << digi->GetChannelNr() << " ) could attach to cluster "
906  << idc << " ( " << fDigiClusters[idc].cluPos << " ) with distance " << digiDigiDist << " . " << fDigiClusters[idc].digiNr.size() << " digis in cluster" << flush;
907 
908  AddDigiToCluster(digiNumber,idc);
909 
910  if ( fVerbose > 1 )
911  cout << " >>> mean channel = " << fDigiClusters[idc].cluPos << endl;
912 
913  if ( iDigiUsed >= 0 ) {
914  //cout << "DIGI " << digiNumber << " MATCHES TO CLUSTER " << iDigiUsed << " and " << idc << endl;
915  JoinTwoClusters(idc,iDigiUsed);
916  }
917 
918  iDigiUsed = idc;
919  }
920  return (iDigiUsed==-1?kFALSE:kTRUE);
921 
922 }
923 // -------------------------------------------------------------------------
924 
925 // ----- Private method ------------------------------------------
926 void PndGemFindClusters::JoinTwoClusters(Int_t clus1, Int_t clus2) {
927  if ( fVerbose > 0 ) cout << "PndGemFindClusters::JointTwoClusters()" << endl;
928 
929  //PndGemDigi* digi; //[R.K. 01/2017] unused variable?
930 
931  for ( size_t id2 = 0 ; id2 < fDigiClusters[clus2].digiNr.size() ; id2++ ) {
932  Int_t digi2 = fDigiClusters[clus2].digiNr[id2];
933  Int_t chan2 = fDigiClusters[clus2].chanNr[id2];
934 
935  Int_t addDigi = kTRUE;
936  for ( size_t id1 = 0 ; id1 < fDigiClusters[clus1].digiNr.size() ; id1++ ) {
937  if ( chan2 == fDigiClusters[clus1].chanNr[id1] ) addDigi = kFALSE;
938  }
939 
940  if ( addDigi ) AddDigiToCluster(digi2,clus1);
941  }
942  fDigiClusters[clus2].cluPos = -1.;
943 }
944 // -------------------------------------------------------------------------
945 
946 
947 // ----- Private method AddDigiToCluster ---------------------------------
948 void PndGemFindClusters::AddDigiToCluster(Int_t digiNr, Int_t clusNr) {
949  if ( fVerbose > 0 ) cout << "PndGemFindClusters::AddDigiToCluster()" << endl;
950 
951  PndGemDigi* digi = (PndGemDigi*) fDigis->At(digiNr);
952  Int_t stationNr = digi->GetStationNr();
953  Int_t sensorNr = digi->GetSensorNr();
954  Int_t iSide = digi->GetSide();
955 
956  PndGemSensor* sensor = fDigiPar->GetSensor(stationNr, sensorNr);
957 
958  if ( fDigiClusters[clusNr].digiNr.size() == 0 ) {
959  // cout << "NEW CLUSTER " << clusNr << " CREATED WITH AddDigiToCluster(" << digiNr << "," << clusNr << ") IN DETECTOR " << digi->GetDetectorId() << endl;
960  fDigiClusters[clusNr].detId = digi->GetDetectorId();
961  fDigiClusters[clusNr].cluPos = digi->GetChannelNr();
962  fDigiClusters[clusNr].cluADC = digi->GetCharge();
963  fDigiClusters[clusNr].cluTDC = digi->GetTimeStamp();
964  fDigiClusters[clusNr].cluPMn = digi->GetChannelNr();
965  fDigiClusters[clusNr].cluPMx = digi->GetChannelNr();
966  fDigiClusters[clusNr].cluMVl = digi->GetCharge();
967  fDigiClusters[clusNr].cluMPs = digi->GetChannelNr();
968  }
969  else {
970  fDigiClusters[clusNr].cluTDC = (fDigiClusters[clusNr].cluTDC*((Double_t)(fDigiClusters[clusNr].sigADC.size()-1.))+digi->GetTimeStamp())/((Double_t)(fDigiClusters[clusNr].sigADC.size()));
971  fDigiClusters[clusNr].cluPos = sensor->GetMeanChannel(iSide,fDigiClusters[clusNr].cluPos,fDigiClusters[clusNr].cluADC,digi->GetChannelNr(),digi->GetCharge());
972  if ( fDigiClusters[clusNr].cluPMn > digi->GetChannelNr() ) fDigiClusters[clusNr].cluPMn = digi->GetChannelNr();
973  if ( fDigiClusters[clusNr].cluPMx < digi->GetChannelNr() ) fDigiClusters[clusNr].cluPMx = digi->GetChannelNr();
974  fDigiClusters[clusNr].cluADC = fDigiClusters[clusNr].cluADC+digi->GetCharge();
975  if ( digi->GetCharge() > fDigiClusters[clusNr].cluMVl ) {
976  fDigiClusters[clusNr].cluMVl = digi->GetCharge();
977  fDigiClusters[clusNr].cluMPs = digi->GetChannelNr();
978  }
979  }
980 
981  fDigiClusters[clusNr].digiNr.push_back(digiNr);
982  fDigiClusters[clusNr].chanNr.push_back((Int_t)digi->GetChannelNr());
983  fDigiClusters[clusNr].sigADC.push_back(digi->GetCharge());
984 }
985 // -------------------------------------------------------------------------
986 
987 // ----- Private method SortClusters -------------------------------
988 void PndGemFindClusters::SortClusters() {
989  if ( fVerbose > 0 ) cout << "PndGemFindClusters::SortClusters()" << endl;
990  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
991  SortCluster(idc);
992  }
993 }
994 // -------------------------------------------------------------------------
995 
996 // ----- Private method SortClusters -------------------------------
997 void PndGemFindClusters::SortCluster(Int_t clus) {
998  if ( fVerbose > 0 ) cout << "PndGemFindClusters::SortCluster()" << endl;
999  Int_t channelPosition = 0;
1000 
1001  for ( Int_t ich = fDigiClusters[clus].cluPMn ; ich <= fDigiClusters[clus].cluPMx ; ich++ ) {
1002  for ( size_t idigi = channelPosition ; idigi < fDigiClusters[clus].digiNr.size() ; idigi++ ) {
1003  if ( fDigiClusters[clus].chanNr[idigi] == ich ) {
1004  Int_t dN = fDigiClusters[clus].digiNr[idigi];
1005  Int_t cN = fDigiClusters[clus].chanNr[idigi];
1006  Double_t sA = fDigiClusters[clus].sigADC[idigi];
1007  fDigiClusters[clus].digiNr[idigi] = fDigiClusters[clus].digiNr[channelPosition];
1008  fDigiClusters[clus].chanNr[idigi] = fDigiClusters[clus].chanNr[channelPosition];
1009  fDigiClusters[clus].sigADC[idigi] = fDigiClusters[clus].sigADC[channelPosition];
1010  fDigiClusters[clus].digiNr[channelPosition] = dN;
1011  fDigiClusters[clus].chanNr[channelPosition] = cN;
1012  fDigiClusters[clus].sigADC[channelPosition] = sA;
1013  channelPosition++;
1014  }
1015  }
1016  }
1017 }
1018 // -------------------------------------------------------------------------
1019 
1020 // ----- Private method PrintClusters -------------------------------
1021 void PndGemFindClusters::PrintClusters() {
1022  if ( fVerbose > 0 ) cout << "PndGemFindClusters::PrintClusters()" << endl;
1023 
1024  for ( size_t idc = 0 ; idc < fDigiClusters.size() ; idc++ ) {
1025  // if ( fDigiClusters[idc].detId != 568329024 ) continue;
1026  if ( fDigiClusters[idc].cluPos < -0.5 ) continue;
1027  PrintCluster(idc);
1028  }
1029 }
1030 // -------------------------------------------------------------------------
1031 
1032 // ----- Private method PrintCluster -------------------------------
1033 void PndGemFindClusters::PrintCluster(Int_t clus) {
1034  if ( fVerbose > 0 ) cout << "PndGemFindClusters::PrintCluster()" << endl;
1035 
1036  cout << "CLUSTER " << clus << " IN " << fDigiClusters[clus].detId << " AT " << fDigiClusters[clus].cluPos
1037  << " CREATED FROM " << fDigiClusters[clus].digiNr.size() << " DIGIS, FROM "
1038  << fDigiClusters[clus].cluPMn << " TO " << fDigiClusters[clus].cluPMx << " WITH MAXIMUM OF "
1039  << fDigiClusters[clus].cluMVl << " AT " << fDigiClusters[clus].cluMPs << endl;
1040 
1041  for ( Int_t icol = 0 ; icol < fDigiClusters[clus].cluPMx-fDigiClusters[clus].cluPMn+3 ; icol++ ) cout << "-" << flush; cout<<endl;
1042  for ( Double_t fh = 0.99 ; fh >= 0. ; fh-=0.247475 ) {
1043  cout << "|" << flush;
1044  Int_t thisDigi = 0;
1045  Double_t thisADC = 0.;
1046  for ( Int_t ichan = fDigiClusters[clus].cluPMn ; ichan <= fDigiClusters[clus].cluPMx ; ichan++ ) {
1047  thisADC = 0.;
1048  if ( fDigiClusters[clus].chanNr[thisDigi] == ichan ) {
1049  thisADC = fDigiClusters[clus].sigADC[thisDigi];
1050  thisDigi++;
1051  }
1052  if ( thisADC > fh*fDigiClusters[clus].cluMVl ) cout << "*" << flush;
1053  else cout << " " << flush;
1054  }
1055  cout << "|" << endl;
1056  }
1057  for ( Int_t icol = 0 ; icol < fDigiClusters[clus].cluPMx-fDigiClusters[clus].cluPMn+3 ; icol++ ) cout << "-" << flush; cout<<endl;
1058 }
1059 // -------------------------------------------------------------------------
1060 
1061 
1062 // ----- Private method Finish ------------------------------------------
1063 void PndGemFindClusters::Finish() {
1064  if ( fClusters ) fClusters->Clear();
1065 
1066  cout << "-------------------- " << fName.Data() << " : Summary -----------------------" << endl;
1067  cout << " Events: " << setw(10) << fTNofEvents << endl;
1068  cout << " Digis: " << setw(10) << fTNofDigis << " ( " << (Double_t)fTNofDigis /((Double_t)fTNofEvents) << " per event )" << endl;
1069  cout << " Clusters: " << setw(10) << fTNofClusters << " ( " << (Double_t)fTNofClusters/((Double_t)fTNofEvents) << " per event )" << endl;
1070  cout << " --> ( " << (Double_t)fTNofDigis /((Double_t)fTNofClusters) << " digis per cluster )" << endl;
1071  cout << "---------------------------------------------------------------------" << endl;
1072 
1073 }
1074 // -------------------------------------------------------------------------
1075 
1076 
1077 ClassImp(PndGemFindClusters)
1078 
DigiCluster::digiNr
std::vector< Int_t > digiNr
Definition: PndGemFindClusters.h:39
PndGemFindClusters::fDigiPar
PndGemDigiPar * fDigiPar
Definition: PndGemFindClusters.h:80
fVerbose
int fVerbose
Definition: poormantracks.C:24
PndGemFindClusters::CompareDigiToClustersDigis
Bool_t CompareDigiToClustersDigis(Int_t digiNumber)
Definition: PndGemFindClusters.cxx:877
digi
TClonesArray * digi
Definition: checkdedx_helixhit.C:20
PndGemFindClusters::fTimeOrderedDigis
std::map< Double_t, Int_t > fTimeOrderedDigis
Definition: PndGemFindClusters.h:84
run
Int_t run
Definition: autocutx.C:47
PndGemFindClusters::fClusters
TClonesArray * fClusters
Definition: PndGemFindClusters.h:82
PndGemDigi::GetStationNr
Int_t GetStationNr() const
Definition: PndGemDigi.h:84
map
PndTransMap * map
Definition: sim_emc_apd.C:99
PndGemFindClusters::PrintCluster
void PrintCluster(Int_t clus)
Definition: PndGemFindClusters.cxx:1033
DigiCluster::chanNr
std::vector< Int_t > chanNr
Definition: PndGemFindClusters.h:40
PndGemFindClusters::fDigis
TClonesArray * fDigis
Definition: PndGemFindClusters.h:81
PndGemFindClusters::Init
virtual InitStatus Init()
Definition: PndGemFindClusters.cxx:193
sensor
TGeoVolume * sensor
Definition: createDummySensor.C:21
PndGemDigiPar
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
PndGemFindClusters::ReInit
virtual InitStatus ReInit()
Definition: PndGemFindClusters.cxx:220
PndGemDigi::GetDetectorId
Int_t GetDetectorId() const
Definition: PndGemDigi.h:79
PndGemStation::GetSensor
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
PndGemFindClusters::SortCluster
void SortCluster(Int_t clus)
Definition: PndGemFindClusters.cxx:997
PndGemDigi::GetSensorNr
Int_t GetSensorNr() const
Definition: PndGemDigi.h:86
PndGemDigiPar::GetSensor
PndGemSensor * GetSensor(Int_t stationNr, Int_t sensorNr)
Definition: PndGemDigiPar.cxx:223
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndGemFindClusters::fDigiClusters
std::vector< DigiCluster > fDigiClusters
Definition: PndGemFindClusters.h:85
PndGemFindClusters::AddDigiToCluster
void AddDigiToCluster(Int_t digiNr, Int_t clusNr)
Definition: PndGemFindClusters.cxx:948
PndGemDigi::GetChannelNr
Double_t GetChannelNr() const
Definition: PndGemDigi.h:80
DigiCluster::sigADC
std::vector< Double_t > sigADC
Definition: PndGemFindClusters.h:41
PndGemFindClusters::SetParContainers
virtual void SetParContainers()
Definition: PndGemFindClusters.cxx:175
PndGemDigiPar::GetNStations
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
name
TString name
Definition: createRootGeoFile.C:82
PndGemFindClusters::JoinTwoClusters
void JoinTwoClusters(Int_t clus1, Int_t clus2)
Definition: PndGemFindClusters.cxx:926
ClassImp
ClassImp(PndAnaContFact)
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
PndGemSensor::GetDistance2
Double_t GetDistance2(Int_t iSide, Double_t chan1, Double_t chan2)
Definition: PndGemSensor.cxx:653
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
PndGemDigi::GetSide
Int_t GetSide() const
Definition: PndGemDigi.h:88
PndGemFindClusters::Exec
virtual void Exec(Option_t *opt)
Definition: PndGemFindClusters.cxx:106
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