All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
SiStripClus.cc
Go to the documentation of this file.
1 // History:
2 //- improved clustering procedure, two different algorithms implemented: COG & analog head-tail, Z. Drasal Apr 2010
3 //- only 1 relation TrackerPulse <--> SimTrackerHit read in, Z. Drasal May 2010
4 //- output relation TrackerHit <--> SimTrackerHit not used, SimTrackerHits (without weights) accessible via TrackerHit->getRawHits() method, Z. Drasal May 2010
5 //- corrected charge units, Z. Drasal May 2010
6 //- TrackerHit -> rawHits() contains 1 SimTrackerHit, which contributed with highest weight, Z. Drasal Oct 2010
7 
8 #include "SiStripClus.h"
9 #include "Colours.h"
10 #include "PhysicalConstants.h"
11 
12 #include "SiStripGeomBuilder.h"
13 
14 #include <iomanip>
15 #include <math.h>
16 #include <stdlib.h>
17 
18 // Include LCIO header files
19 #include <lcio.h>
20 #include <EVENT/LCCollection.h>
21 #include <IMPL/LCCollectionVec.h>
22 #include <IMPL/LCFlagImpl.h>
23 #include <IMPL/LCRelationImpl.h>
24 #include <UTIL/LCRelationNavigator.h>
25 #include <IMPL/TrackerPulseImpl.h>
26 #include <IMPL/TrackerHitPlaneImpl.h>
27 #include <UTIL/CellIDDecoder.h>
28 #include <UTIL/CellIDEncoder.h>
29 #include "UTIL/LCTrackerConf.h"
30 #include <UTIL/ILDConf.h>
31 
32 
33 // Include Marlin
34 #include <streamlog/streamlog.h>
35 
36 // Namespaces
37 using namespace CLHEP;
38 
39 namespace sistrip {
40 
41 //
42 // Instantiate this object
43 //
44 SiStripClus anSiStripClus;
45 
46 //
47 // Constructor
48 //
49 SiStripClus::SiStripClus() : Processor("SiStripClus")
50 {
51  // Processor description
52  _description = "SiStripClus: Marlin processor intended for cluster finding - using digitized data worked out with SiStripDigi";
53 
54  // Define compulsory parameters
55  registerProcessorParameter( "Subdetector",
56  "Subdetector to be digitised",
57  _subdetector,
58  std::string("FTD") );
59 
60  registerProcessorParameter( "CMSnoise",
61  "Common mode subtracted noise, set in electrons",
62  _CMSnoise,
63  float(1000) );
64 
65  registerProcessorParameter( "FloatingStripsRPhi",
66  "Is every even strip floating in R-Phi?",
67  _floatStripsRPhi,
68  bool(false));
69 
70  registerProcessorParameter( "FloatingStripsZ",
71  "Is every even strip floating in Z?",
72  _floatStripsZ,
73  bool(false));
74 
75  std::vector<float> resSVDFirstInRPhi;
76  resSVDFirstInRPhi.push_back( 6.5); // 20 degrees
77  resSVDFirstInRPhi.push_back( 6.0); // 30 degrees
78  resSVDFirstInRPhi.push_back( 6.0); // 40 degrees
79 
80  resSVDFirstInRPhi.push_back( 7.0);
81  resSVDFirstInRPhi.push_back( 7.0);
82  resSVDFirstInRPhi.push_back( 6.0); // 50 , 60 , 70 degrees
83 
84  resSVDFirstInRPhi.push_back( 6.5);
85  resSVDFirstInRPhi.push_back( 6.5);
86  resSVDFirstInRPhi.push_back( 6.5); // 80 , 90 , 100 degrees
87 
88  resSVDFirstInRPhi.push_back( 6.5);
89  resSVDFirstInRPhi.push_back( 6.5);
90  resSVDFirstInRPhi.push_back( 6.5); // 110, 120, 130 degrees
91 
92  resSVDFirstInRPhi.push_back( 6.5);
93  resSVDFirstInRPhi.push_back( 6.0); // 140, 150 degrees
94 
95  registerProcessorParameter( "ResolutionOfSVDFirstInRPhi",
96  "Resolution of first SVD layer in R-Phi (in um)",
97  _resSVDFirstInRPhi,
98  resSVDFirstInRPhi);
99 
100  std::vector<float> resSVDOtherInRPhi;
101  resSVDOtherInRPhi.push_back( 7.5);
102  resSVDOtherInRPhi.push_back( 9.5);
103  resSVDOtherInRPhi.push_back(10.5); // 20 , 30 , 40 degrees
104 
105  resSVDOtherInRPhi.push_back(11.0);
106  resSVDOtherInRPhi.push_back(11.0);
107  resSVDOtherInRPhi.push_back(11.0); // 50 , 60 , 70 degrees
108 
109  resSVDOtherInRPhi.push_back(11.0);
110  resSVDOtherInRPhi.push_back(11.0);
111  resSVDOtherInRPhi.push_back(11.0); // 80 , 90 , 100 degrees
112 
113  resSVDOtherInRPhi.push_back(11.0);
114  resSVDOtherInRPhi.push_back(11.0);
115  resSVDOtherInRPhi.push_back(11.0); // 110, 120, 130 degrees
116 
117  resSVDOtherInRPhi.push_back(11.0);
118  resSVDOtherInRPhi.push_back(10.0); // 140, 150 degrees
119 
120  registerProcessorParameter( "ResolutionOfSVDOtherInRPhi",
121  "Resolution of other SVD layers in R-Phi (in um) ",
122  _resSVDOtherInRPhi,
123  resSVDOtherInRPhi);
124 
125  std::vector<float> resSVDFirstInZ;
126  resSVDFirstInZ.push_back(20.0);
127  resSVDFirstInZ.push_back(16.0);
128  resSVDFirstInZ.push_back(14.0); // 20 , 30 , 40 degrees
129 
130  resSVDFirstInZ.push_back(14.0);
131  resSVDFirstInZ.push_back(13.0);
132  resSVDFirstInZ.push_back(13.0); // 50 , 60 , 70 degrees
133 
134  resSVDFirstInZ.push_back(28.0);
135  resSVDFirstInZ.push_back(28.0);
136  resSVDFirstInZ.push_back(28.0); // 80 , 90 , 100 degrees
137 
138  resSVDFirstInZ.push_back(14.0);
139  resSVDFirstInZ.push_back(13.0);
140  resSVDFirstInZ.push_back(13.0); // 110, 120, 130 degrees
141 
142  resSVDFirstInZ.push_back(14.0);
143  resSVDFirstInZ.push_back(16.0); // 140, 150 degrees
144 
145  registerProcessorParameter( "ResolutionOfSVDFirstInZ",
146  "Resolution of first SVD layer along Z axis (in um)",
147  _resSVDFirstInZ,
148  resSVDFirstInZ);
149 
150  std::vector<float> resSVDOtherInZ;
151  resSVDOtherInZ.push_back(16.0);
152  resSVDOtherInZ.push_back(16.0);
153  resSVDOtherInZ.push_back(15.0); // 20 , 30 , 40 degrees
154 
155  resSVDOtherInZ.push_back(14.0);
156  resSVDOtherInZ.push_back(17.0);
157  resSVDOtherInZ.push_back(34.0); // 50 , 60 , 70 degrees
158 
159  resSVDOtherInZ.push_back(50.0);
160  resSVDOtherInZ.push_back(50.0);
161  resSVDOtherInZ.push_back(50.0); // 80 , 90 , 100 degrees
162 
163  resSVDOtherInZ.push_back(34.0);
164  resSVDOtherInZ.push_back(14.0);
165  resSVDOtherInZ.push_back(14.0); // 110, 120, 130 degrees
166 
167  resSVDOtherInZ.push_back(16.0);
168  resSVDOtherInZ.push_back(16.0); // 140, 150 degrees
169 
170  registerProcessorParameter( "ResolutionOfSVDOtherInZ",
171  "Resolution of other SVD layers along Z axis (in um) ",
172  _resSVDOtherInZ,
173  resSVDOtherInZ);
174 
175  registerProcessorParameter( "S/NSeedStrips",
176  "Signal to noise ratio cut for seed strips",
177  _SNseed,
178  float(5) );
179 
180  registerProcessorParameter( "S/NAdjacentStrips",
181  "Signal to noise ratio cut for adjacent strips",
182  _SNadjacent,
183  float(3) );
184 
185  registerProcessorParameter( "S/NCluster",
186  "Signal to noise ratio cut for total cluster",
187  _SNtotal,
188  float(8) );
189 
190  // registerProcessorParameter( "TanOfAvgELorentzShift",
191  // "Tangent of electrons' average Lorentz shift (0.20 for 273 K, 0.175 for 300 K)",
192  // _TanOfAvgELorentzShift,
193  // float(0.175) );
194 
195  registerProcessorParameter( "TanOfAvgHLorentzShift",
196  "Tangent of holes' average Lorentz shift (0.05 for 273 K, 0.039 for 300 K)",
197  _TanOfAvgHLorentzShift,
198  float(0.039) );
199 
200  registerProcessorParameter( "InputCollectionName",
201  "Name of TrackerPulse input collection",
202  _inColName,
203  std::string("FTDDigits") );
204 
205  registerProcessorParameter( "OutputCollectionName",
206  "Name of TrackerHitPlane output collection",
207  _outColName,
208  std::string("FTDTrackerHits") );
209 
210  registerProcessorParameter( "RelCollectionNamePlsToSim",
211  "Name of input relation collection - TrackerPulse to SimTrackerHit (if nonzero, required)",
212  _relColNamePlsToSim,
213  std::string("FTDDigitsToSimHitsRel") );
214 
215  registerProcessorParameter( "pitchFront",
216  "Pitch of the front sensor family in the center of the sensors, set in um",
217  _pitchFront,
218  double(50.0));
219 
220  registerProcessorParameter( "pitchRear",
221  "Pitch of the rear sensor family in the center of the sensors, set in um",
222  _pitchRear,
223  double(50.0));
224 }
225 
226 //
227 // Method called at the beginning of data processing
228 //
229 void SiStripClus::init()
230 {
231  // Set variables in appropriate physical units
232  _CMSnoise *= e;
233 
234  for (unsigned int i=0; i<_resSVDFirstInRPhi.size(); ++i)
235  {
236  _resSVDFirstInRPhi[i] *= um;
237  _resSVDOtherInRPhi[i] *= um;
238  _resSVDFirstInZ[i] *= um;
239  _resSVDOtherInZ[i] *= um;
240  }
241 
242  _pitchFront *= um;
243  _pitchRear *= um;
244 
245 
246  // Get geometry parameters from Gear xml file
247  _geometry = SiStripGeomBuilder::Build("FTD",_pitchFront,_pitchRear);
248  // FIXME: Provisional para compilar
249  _geometry -> initGearParams();
250 
251  // Print set parameters
252  printProcessorParams();
253 
254  //
255  // ROOT variables
256  //
257 #ifdef ROOT_OUTPUT
258  _rootFile = new TFile("BelleII_SVD_Hits.root","recreate");
259 
260  _rootFile->cd("");
261 
262  // Declare Tree
263  _rootTree = new TTree("Hits","Hit info");
264 
265  _rootTree->Branch("Layer" ,&_rootLayerID ,"Layer/I" );
266  _rootTree->Branch("Ladder" ,&_rootLadderID ,"Ladder/I" );
267  _rootTree->Branch("Sensor" ,&_rootSensorID ,"Sensor/I" );
268  _rootTree->Branch("RPhiSim" ,&_rootSimRPhi ,"SimRPhi/D" );
269  _rootTree->Branch("RPhiRec" ,&_rootRecRPhi ,"RecRPhi/D" );
270  _rootTree->Branch("ResRPhi" ,&_rootResRPhi ,"ResRPhi/D" );
271  _rootTree->Branch("ResR" ,&_rootResR ,"ResR/D" );
272  _rootTree->Branch("ResModule" ,&_rootResModule ,"Resmodule/D" );
273  _rootTree->Branch("RPhiClsSize",&_rootClsSizeRPhi ,"ClsSizeRPhi/D");
274  _rootTree->Branch("RPhiMCPDG" ,&_rootMCPDGRPhi ,"MCPDGRPhi/I" );
275  _rootTree->Branch("ZSim" ,&_rootSimZ ,"SimZ/D" );
276  _rootTree->Branch("ZRec" ,&_rootRecZ ,"RecZ/D" );
277  _rootTree->Branch("ZClsSize" ,&_rootClsSizeZ ,"ClsSizeZ/D" );
278  _rootTree->Branch("ZMCPDG" ,&_rootMCPDGZ ,"MCPDGZ/I" );
279  _rootTree->Branch("EvtNum" ,&_rootEvtNum ,"EvtNum/I" );
280 #endif
281 }
282 
283 //
284 // Method called for each run
285 //
286 void SiStripClus::processRunHeader(LCRunHeader * run)
287 {
288 // Print run number
289 // streamlog_out(MESSAGE3) << DGREEN
290 // << " Processing run: "
291 // << ENDCOLOR
292 // << (run->getRunNumber()+1)
293 // << std::endl << std::endl;
294 
295  _nRun++;
296 }
297 
298 //
299 // Method called for each event
300 //
301 void SiStripClus::processEvent(LCEvent * event)
302 {
303  // Get names of all collections saved in LCIO file
304  ConstStringVec * strVec = event->getCollectionNames();
305 
306  // Initialize relation navigators
307  _navigatorPls = NULL;
308 
309  // Initialize
310 #ifdef ROOT_OUTPUT
311  _rootEvtNum = event->getEventNumber();
312 #endif
313 
314  //
315  // Get TrackerPulse collections and relation to MCParticles collection
316  LCCollection * colOfTrkPulses = 0;
317 
318  // Print header info about collections (for 1. event)
319  if(event->getEventNumber() == 0)
320  {
321  streamlog_out(MESSAGE3) << " "
322  << DUNDERL
323  << DBLUE
324  <<"LCCollection(s) processed:"
325  << ENDCOLOR
326  << std::endl << std::endl;
327  }
328 
329  // Go through all collection names
330  for(ConstStringVec::const_iterator colName = strVec->begin();
331  colName != strVec->end(); ++colName)
332  {
333  // Tracker pulses
334  if(_inColName == (*colName))
335  {
336  // Collection must be of TrackerPulse type
337  if( (event->getCollection(*colName))->getTypeName()
338  == LCIO::TRACKERPULSE )
339  {
340  // Save collection
341  colOfTrkPulses = event->getCollection(*colName);
342  // Print info
343  if (event->getEventNumber() == 0)
344  {
345  streamlog_out(MESSAGE3) << " " << *colName
346  << std::endl;
347  }
348  }
349  // Collection NOT of TrackerPulse type
350  else
351  {
352  streamlog_out(ERROR) << "Required collection: "
353  <<_inColName
354  << " found, but NOT of TRACKERPULSE type!!!"
355  << std::endl;
356  exit(0);
357  }
358  }
359 
360  // Relation TrkPulses to SimTrackerHits
361  if( (!_relColNamePlsToSim.empty()) && (_relColNamePlsToSim == (*colName)) )
362  {
363  // Collection must be of relation type
364  if( (event->getCollection(*colName))->getTypeName()
365  == LCIO::LCRELATION )
366  {
367  // Save collection into relation navigator
368  _navigatorPls = new LCRelationNavigator(
369  event->getCollection(*colName));
370  // Print info
371  if(event->getEventNumber() == 0)
372  {
373  streamlog_out(MESSAGE3) << " " << *colName
374  << std::endl;
375  }
376  }
377  // Collection NOT of relation type
378  else
379  {
380  streamlog_out(ERROR) << "Required collection: "
381  <<_relColNamePlsToSim
382  << " found, but NOT of LCRELATION type!!!"
383  << std::endl;
384  exit(0);
385  }
386  }
387  } // For - collection names
388 
389  // Input collection NOT found
390  if(colOfTrkPulses == 0)
391  {
392  streamlog_out(WARNING) << "Required collection: "
393  <<_inColName
394  << " not found!!!"
395  << std::endl;
396  }
397 
398  if( (!_relColNamePlsToSim.empty()) && (_navigatorPls == NULL) )
399  {
400  streamlog_out(WARNING) << "Required collection: "
401  <<_relColNamePlsToSim
402  << " not found!!!"
403  << std::endl;
404  }
405 
406  //
407  // Process TrackerPulse collections
408  if(colOfTrkPulses != 0)
409  {
410  CellIDDecoder<TrackerPulse> cellIDDec(colOfTrkPulses);
411 
412  // Initialize variables
413  TrackerPulseImpl * pulse = 0;
414 
415  // Sensor map with map of hit strips (corresponding to the given sensor)
416  SensorStripMap sensorMap;
417 
418  // Get number of elements in each collection
419  int nPulses = colOfTrkPulses->getNumberOfElements();
420 
421  // Process pulses
422  for(int i=0; i<nPulses; ++i)
423  {
424  // Copy the content of collection TrackerPulse to a pulse
425  pulse = dynamic_cast<TrackerPulseImpl*>(
426  colOfTrkPulses->getElementAt(i) );
427 
428  // Update the sensor strip map with the pulse
429  updateMap(pulse, sensorMap);
430  }
431 
432  //
433  // Find clusters
434  ClsVec clsVec = findClus(sensorMap);
435 
436  // Releasing memory and clearing
437  releaseMap(sensorMap);
438 
439  //
440  // Calculate real + ghost hits from clusters - in global ref. system +
441  // create relations to MCParticles
442  IMPL::LCCollectionVec * colOfTrkHits =
443  new IMPL::LCCollectionVec(LCIO::TRACKERHITPLANE);
444 
445  // Calculate the hits
446  calcHits(clsVec, colOfTrkHits);
447 
448  //
449  // Save the collection (vector) of hits + relation to MC
450  event->addCollection(colOfTrkHits, _outColName);
451 
452  } // If (colOfTrkPulses != 0)
453  // Release memory
454  if(_navigatorPls != NULL)
455  {
456  delete _navigatorPls;
457  _navigatorPls = NULL;
458  }
459 
460  _nEvent++;
461 }
462 
463 //
464 // Method called after each event to check the data processed
465 //
466 void SiStripClus::check(LCEvent * event)
467 {
468 }
469 
470 //
471 // Method called after all data processing
472 //
473 void SiStripClus::end()
474 {
475  // Release memory
476  delete _geometry;
477  _geometry = 0;
478 
479 #ifdef ROOT_OUTPUT
480 
481  // Close file
482  _rootFile->cd("");
483  _rootFile->Write();
484  _rootFile->Close();
485 
486 #endif
487 
488  // Print message
489  streamlog_out(MESSAGE3) << std::endl
490  << " "
491  << DGREEN
492  << "Processor successfully finished!"
493  << ENDCOLOR
494  << std::endl;
495 
496 }
497 
498 // MAIN CLUSTER METHOD
499 
500 //
501 // Method searching for clusters
502 //
503 typedef std::pair<int,StripCluster*> StripClusterPair;
504 typedef std::map<int,std::map<StripType,std::vector<StripClusterPair> > > SensorStripClusterMap;
505 
506 ClsVec SiStripClus::findClus(SensorStripMap & sensorMap)
507 {
508  ClsVec clsVec;
509 
510  // Cluster vectors -
511  SensorStripClusterMap clsvectFrontRear;
512 
513  std::vector<StripType> stv;
514  stv.push_back(STRIPFRONT);
515  stv.push_back(STRIPREAR);
516  // As the sensorID was lost (see SiStripClus::updateMap method)
517  // assigning it
518  std::map<StripType,int> stSensorIDmap;
519  stSensorIDmap[STRIPFRONT] = 1;
520  stSensorIDmap[STRIPREAR] = 3;
521 
522  // Bunch of strips forming cluster
523  StripChargeMap clsStrips;
524 
525  //
526  // Search complete sensor map - find seeds & their neghbouring strips
527  for(SensorStripMap::iterator iterSMap = sensorMap.begin(); iterSMap!=sensorMap.end();
528  ++iterSMap)
529  {
530  // Save layer ID , ...
531  const int cellID = iterSMap->first;
532  std::map<std::string,int> bfmap = _geometry->decodeCellID(cellID);
533  const int layerID = bfmap["layer"];
534  const int ladderID= bfmap["module"];
535 
536  // Storing the two types of clusters
537  for(std::vector<StripType>::iterator itST = stv.begin(); itST!= stv.end();
538  ++itST)
539  {
540  StripType STRIPTYPE = *itST;
541  for(StripChargeMap::iterator iterChMap=iterSMap->second[STRIPTYPE].begin();
542  iterChMap!=iterSMap->second[STRIPTYPE].end(); iterChMap++)
543  {
544  const int sensorID= stSensorIDmap[STRIPTYPE];
545 
546  // Begin algorithm
547  // Zero: Save new candidate for seed strip + MC true info
548  const int seedStrip = iterChMap->first;
549  const double seedCharge = iterChMap->second->getCharge();
550 
551  const SimTrackerHitMap & seedSimHitMap=iterChMap->second->getSimHitMap();
552 
553  if( seedCharge < (_SNseed*_CMSnoise) )
554  {
555  continue;
556  }
557 
558  // First: New cluster and its seed strip has been found
559  clsStrips[seedStrip] = new Signal(seedCharge, 0.);
560  clsStrips[seedStrip]->updateSimHitMap(seedSimHitMap);
561 
562  // Set this charge as zero to avoid double counting
563  iterChMap->second->setCharge(0.);
564 
565  // Map ordered from lower to higher
566  // Continue searching - find left and right neighbours
567  // Taking account also the floating strips case...
568  StripChargeMap::iterator seedIt=iterSMap->second[STRIPTYPE].find(seedStrip);
569  // Second: search for left neighbours
570  //-- Left neighbours
571  StripChargeMap::reverse_iterator lschMap(seedIt);
572  clsStrips = storeHitsAdjacents<StripChargeMap::reverse_iterator>(
573  clsStrips, lschMap, iterSMap->second[STRIPTYPE].rend(),
574  iterSMap->second[STRIPTYPE]);
575 
576  // Third: search for rigth neighbours
577  //-- Right neighbours
578  StripChargeMap::iterator rschMap(++seedIt);
579  clsStrips = storeHitsAdjacents<StripChargeMap::iterator>(clsStrips,
580  rschMap, iterSMap->second[STRIPTYPE].end(),
581  iterSMap->second[STRIPTYPE]);
582 
583  // Fourth: Calculate mean position of a new cluster
584  SimTrackerHitMap clsSimHitMap;
585 
586  // Cluster: position, charge & size
587  double clsCharge = 0.0;
588 
589  double xLeftSignal = 0.0;
590  double qLeftSignal = 0.0;
591 
592  double xRightSignal = 0.0;
593  double qRightSignal = 0.0;
594 
595  double qIntermSignal = 0.0;
596 
597  int stripID = 0;
598  for(StripChargeMap::iterator iterChMap2=clsStrips.begin();
599  iterChMap2!=clsStrips.end(); iterChMap2++)
600  {
601  // Current strip ID, posZ & charge
602  stripID = iterChMap2->first;
603  const double stripPosYatz0= _geometry->getStripPosY(layerID,
604  sensorID,stripID,0.0);
605 
606  const double stripCharge = iterChMap2->second->getCharge();
607 
608  // Update info about MC particles which contributed
609  const SimTrackerHitMap & simHitMap =
610  iterChMap2->second->getSimHitMap();
611 
612  if(simHitMap.size() != 0)
613  {
614  for(SimTrackerHitMap::const_iterator iterSHM=simHitMap.begin();
615  iterSHM!=simHitMap.end(); ++iterSHM)
616  {
617  EVENT::SimTrackerHit * simHit =
618  dynamic_cast<EVENT::SimTrackerHit *>(iterSHM->first);
619  float weight = iterSHM->second;
620  if(clsSimHitMap.find(simHit)!=clsSimHitMap.end())
621  {
622  clsSimHitMap[simHit] += weight;
623  }
624  else
625  {
626  clsSimHitMap[simHit] = weight;
627  }
628  }
629  }
630 
631  // Get last iterator
632 
633  // Get leftmost signal
634  if(iterChMap2 == clsStrips.begin())
635  {
636  xLeftSignal = stripPosYatz0;
637  qLeftSignal = stripCharge;
638  }
639 
640  // Get rightmost signal
641  else if(iterChMap2 == --(clsStrips.end()))
642  {
643  xRightSignal = stripPosYatz0;
644  qRightSignal = stripCharge;
645  }
646  // Get intermediate signal
647  else
648  {
649  qIntermSignal += stripCharge;
650  }
651 
652  // Update total charge
653  clsCharge += stripCharge;
654 
655  }
656 
657  // Number of strips being part of cluster
658  int clsSize = clsStrips.size();
659  // Get average intermediate signal
660  if (clsSize > 2)
661  {
662  qIntermSignal /= (clsSize - 2.0);
663  }
664  else
665  {
666  qIntermSignal = 0.;
667  }
668 
669  // Building the cluster
670  double geomPitchAtz0 = _geometry->getSensorPitch(layerID,sensorID,0.0);
671  double readOutPitch = 0.;
672  if(_floatStripsZ) // FIXME---> CAMBIAR POR mapa _floatStrips[type]
673  {
674  readOutPitch = 2.0*geomPitchAtz0;
675  }
676  else
677  {
678  readOutPitch = geomPitchAtz0;
679  }
680 
681  double clsPosYatz0 = 0.0;
682  // Analog head-tail algorithm
683  // qIntermSignal >= 1/eps*qRighSignal || 1/eps*qLeftSignal
684  // --> if not don't use head-tail (delta electron ...) - use
685  // epsilon ~ 1.5
686  if( (clsSize>2) && (qLeftSignal<1.5*qIntermSignal)
687  && (qRightSignal<1.5*qIntermSignal) )
688  {
689  clsPosYatz0 = (xRightSignal + xLeftSignal)/2.
690  + (qRightSignal - qLeftSignal)/2./qIntermSignal * readOutPitch;
691  }
692  else // COG algorithm
693  {
694  clsPosYatz0 = (xRightSignal*qRightSignal
695  + xLeftSignal*qLeftSignal)/(qRightSignal + qLeftSignal);
696  }
697 
698 
699  //
700  // Fifth: Correct mean position to Lorentz shift
701  // As theta is Pi/2.0 --> this correction is 0: FIXME-> TO BE REMOVED
702  clsPosYatz0 += _TanOfAvgHLorentzShift
703  * _geometry->getSensorThick(layerID)/2.
704  * cos(_geometry->getLadderTheta(layerID));
705 
706  //
707  // Sixth: Save information about new cluster
708  if(clsCharge >= (_SNtotal*_CMSnoise))
709  {
710  StripCluster * pCluster = new StripCluster( layerID,
711  ladderID, sensorID,
712  Hep3Vector(_geometry->getSensorThick(layerID)/2.,
713  clsPosYatz0, 0.0),
714  Hep3Vector(_geometry->getSensorThick(layerID)/2., 0., 0.),
715  clsCharge, clsSize );
716  pCluster->updateSimHitMap(clsSimHitMap);
717 
718  clsvectFrontRear[cellID][STRIPTYPE].push_back(
719  std::pair<int,StripCluster*>(stripID,pCluster));
720  }
721 
722  // Release memory
723  for(StripChargeMap::iterator iterChMap2=clsStrips.begin();
724  iterChMap2!=clsStrips.end(); iterChMap2++)
725  {
726  delete iterChMap2->second; //Signal pointers
727  }
728 
729  // Clear content
730  clsStrips.clear();
731  }
732  } // For cluster strip type (front-rear)
733  } // For sensor map
734 
735 
736  // Stores the hit (using the STRIPFRONT as init)
737  for(SensorStripClusterMap::iterator it = clsvectFrontRear.begin();
738  it != clsvectFrontRear.end(); ++it)
739  {
740  const int cellID = it->first;
741  std::map<std::string,int> bfmap = _geometry->decodeCellID(cellID);
742  const int layerID = bfmap["layer"];
743  const int ladderID= bfmap["module"];
744 
745  // the front sensors
746  for(std::vector<StripClusterPair>::iterator
747  itFront = it->second[STRIPFRONT].begin();
748  itFront != it->second[STRIPFRONT].end(); ++itFront)
749  {
750  const int stripIDFront = itFront->first;
751  // Extracting the points in the edges
752  const double yatz0Front = _geometry->getStripPosY(layerID,1,
753  stripIDFront,0.0);
754  const double yatzLFront = _geometry->getStripPosY(layerID,1,
755  stripIDFront,_geometry->getSensorLength(layerID));
756 
757 
758  // Checking with the Rear sensors
759  for(std::vector<StripClusterPair>::iterator
760  itRear = it->second[STRIPREAR].begin();
761  itRear != it->second[STRIPREAR].end(); ++itRear)
762  {
763  const int stripIDRear = itRear->first;
764  // Extracting the points in the edges and put them in a common
765  // reference system (Front petal local ref., remember the y=0
766  // of one system is SensorWidthMax distance from the other)
767  const double yatz0Rear = _geometry->getSensorWidthMax(layerID)-
768  _geometry->getStripPosY(layerID,3,stripIDRear,0.0);
769  const double yatzLRear = _geometry->getSensorWidthMax(layerID)-
770  _geometry->getStripPosY(layerID,3,stripIDRear,
771  _geometry->getSensorLength(layerID));
772 
773  // If not crossing continue
774  if( (yatz0Front-yatz0Rear)*(yatzLFront-yatzLRear) > 0.)
775  {
776  continue;
777  }
778 
779  // There are intersection, so store the Hit
780  // Cluster in Front
781  StripCluster * pclusterFront = itFront->second;
782  // Cluster in Rear
783  StripCluster * pclusterRear = itRear->second;
784 
785  // Get the intersection point (note that the returning value
786  // is placed in the FRONT sensor and in the x_local = thick/2
787  CLHEP::Hep3Vector position=_geometry->getCrossLinePoint(layerID,
788  ladderID,stripIDFront,stripIDRear);
789 
790  // FIXME: Sigma calculation!??
791  const double sigmaY = sqrt(
792  pclusterFront->getPosY()*pclusterFront->getPosY()+
793  pclusterRear->getPosY()*pclusterRear->getPosY());
794  // The Z indetermination: putting z in the middle of the disk
795  // (i.e. middle position between front and rear sensors)
796  // FIXME---> Hay que ver que es esto ???!
797  const double sigmaZ = _geometry->getSensorThick(layerID);
798  Hep3Vector posSigma(_geometry->getSensorThick(layerID)/2.,
799  sigmaY, sigmaZ);
800 
801  double totalCharge = ( pclusterFront->getCharge() +
802  pclusterRear->getCharge())/2.0;
803 
804  // Always postioned as it was in the front
805  StripCluster * pCluster3D = new StripCluster( layerID, ladderID,
806  1, position, posSigma, totalCharge, 0);
807 
808  // Update MC true info for 3D
809  SimTrackerHitMap cls3DSimHitMap = pclusterFront->getSimHitMap();
810 
811  if(cls3DSimHitMap.size() != 0)
812  {
813  SimTrackerHitMap clsRearSimHitMap=
814  pclusterRear->getSimHitMap();
815  for(SimTrackerHitMap::iterator iterSHM=
816  clsRearSimHitMap.begin();
817  iterSHM!=clsRearSimHitMap.end(); iterSHM++)
818  {
819  EVENT::SimTrackerHit * simHit =
820  dynamic_cast<EVENT::SimTrackerHit *>(iterSHM->first);
821  float weight = iterSHM->second;
822  if(cls3DSimHitMap.find(simHit)!=cls3DSimHitMap.end())
823  {
824  cls3DSimHitMap[simHit] += weight;
825  }
826  else
827  {
828  cls3DSimHitMap[simHit] = weight;
829  }
830  }
831  pCluster3D->updateSimHitMap(cls3DSimHitMap);
832  }
833  //
834  // if defined ROOT-OUTPUT, save info
835 #ifdef ROOT_OUTPUT
836  // Set layer, ladder, sensor ID
837  _rootLayerID = layerID;
838  _rootLadderID = ladderID;
839  _rootSensorID = 0;
840 
841  // Set reconstructed position (local frame)
842  //_rootRecRPhi = pCluster3D->getPosY() / mm;
843  //_rootRecZ = pCluster3D->getPosZ() / mm;
844  // Set cluster sizes
845  _rootClsSizeRPhi = pclusterFront->getSize();
846  _rootClsSizeZ = pclusterRear->getSize();
847 
848 
849  // Set simulated position in Front & MC particle
850 
851  // Find SimTrackerHit with highest weight
852  EVENT::SimTrackerHit * simHit = 0;
853  float weight = 0;
854 
855  const SimTrackerHitMap & simHitMapFront =
856  pclusterFront->getSimHitMap();
857  for(SimTrackerHitMap::const_iterator iterSHM =
858  simHitMapFront.begin();
859  iterSHM!=simHitMapFront.end(); ++iterSHM)
860  {
861  // Find contribution with highest weight
862  if( (iterSHM->first!=0) && ((iterSHM->second)>weight) )
863  {
864  simHit = iterSHM->first;
865  weight = iterSHM->second;
866  }
867  }
868 
869  // SimHit global position
870  Hep3Vector simPosGlob;
871  Hep3Vector simPosLoc;
872 
873  simPosGlob.setX(simHit->getPosition()[0]*mm);
874  simPosGlob.setY(simHit->getPosition()[1]*mm);
875  simPosGlob.setZ(simHit->getPosition()[2]*mm);
876  simPosLoc = _geometry->transformPointToLocal(layerID,
877  ladderID, 1, simPosGlob);
878  _rootMCPDGRPhi = simHit->getMCParticle()->getPDG();
879  _rootSimRPhi = sqrt(simPosGlob.getY()*simPosGlob.getY()+
880  simPosGlob.getX()*simPosGlob.getX())/mm;
881  _rootSimZ = simPosGlob.getZ() / mm;
882 
883  // Residuals: ref. frame difined in the simHit-- rPhi, r
884  const double theta = atan2(simPosGlob.getY(),simPosGlob.getX());
885  //-- getting the reconstructed hit to the global ref.
886  CLHEP::Hep3Vector recPoint = _geometry->transformPointToGlobal(
887  layerID,ladderID,1,position);
888  _rootRecRPhi = sqrt(recPoint.getX()*recPoint.getX()
889  + recPoint.getY()*recPoint.getY())/mm;
890  _rootRecZ = recPoint.getZ()/mm;
891 
892  recPoint -= simPosGlob;
893  recPoint.rotateZ(theta); // FIXME: Esto es correcto??
894  // Extracting the residuals: In the local frame
895  // Sim position - rec position
896  _rootResRPhi = (simPosLoc.getY()-position.getY())/mm;
897  _rootResR = (simPosLoc.getZ()-position.getZ())/mm;
898  _rootResModule= sqrt(_rootResRPhi*_rootResRPhi+
899  _rootResR*_rootResR)/mm;
900 
901  //FIXME: Que hago con esta doble informacion???
902 /* // Set simulated position in Rear & MC particle
903 
904  // Find SimTrackerHit with highest weight
905  const SimTrackerHitMap & simHitMapRear =
906  pclusterRear->getSimHitMap();
907 
908  for(SimTrackerHitMap::const_iterator iterSHM=
909  simHitMapRear.begin();
910  iterSHM!=simHitMapRear.end(); ++iterSHM)
911  {
912  // Find contribution with highest weight
913  if( (iterSHM->first!=0) && ((iterSHM->second)>weight) )
914  {
915  simHit = iterSHM->first;
916  weight = iterSHM->second;
917  }
918  }
919 
920  // SimHit global position
921  simPosGlob.setX(simHit->getPosition()[0]*mm);
922  simPosGlob.setY(simHit->getPosition()[1]*mm);
923  simPosGlob.setZ(simHit->getPosition()[2]*mm);
924  simPosLoc = _geometry->transformPointToLocal(layerID, ladderID,
925  3, simPosGlob);
926  _rootMCPDGZ = simHit->getMCParticle()->getPDG();
927  _rootSimZ = simPosLoc.getZ() / mm;*/
928 
929  // Fill the tree
930  _rootFile->cd("");
931  _rootTree->Fill();
932 
933 #endif
934 
935  clsVec.push_back(pCluster3D);
936  } // for rear sensors
937  } // for front sensors
938  } // For cluster strip vector pairs
939 
940 /* for(SensorStripClusterMap::iterator it = clsvectFrontRear.begin();
941  it != clsvectFrontRear.end(); ++it)
942  {
943  std::cout << "========================================== " << std::endl;
944  std::cout << "CELLID: " << it->first << std::endl;
945  for(std::map<StripType,std::vector<StripClusterPair> >::iterator itM
946  = it->second.begin(); itM != it->second.end(); ++itM)
947  {
948  std::cout << " Strip Type: " << itM->first << std::endl;
949  for(std::vector<StripClusterPair>::iterator itP = itM->second.begin();
950  itP != itM->second.end(); ++itP)
951  {
952  std::cout << " Cluster Info " << std::endl;
953  std::cout << " - Position[mm]:"
954  << itP->second->get3Position()/mm << std::endl;
955  std::cout << " - Charge :"
956  << itP->second->getCharge() << std::endl;
957  }
958  }
959  }
960  */
961 
962  // Release memory
963  for(SensorStripClusterMap::iterator it = clsvectFrontRear.begin();
964  it != clsvectFrontRear.end(); ++it)
965  {
966  for(std::map<StripType,std::vector<StripClusterPair> >::iterator itMap =
967  it->second.begin();
968  itMap != it->second.end(); ++itMap)
969  {
970  for(std::vector<StripClusterPair>::iterator itSCl = itMap->second.begin();
971  itSCl != itMap->second.end(); ++itSCl)
972  {
973  if( itSCl->second != 0)
974  {
975  delete itSCl->second;
976  }
977  }
978  }
979  }
980 
981 
982  return clsVec;
983 }
984 
985  //-- Right neighbours
986 /* adjCharge = 0.0;
987  goNextStrip = true;
988  while( goNextStrip && LschMap != iterSMap->second[STRIPTYPE].end() )
989  {
990  adjCharge = LschMap->getCharge();
991  const SimTrackerHitMap & adjSimHitMap = LschMap->getSimHitMap();
992 std::cout << " Strip: " << LschMap->first<< " -- carga:" << adjCharge << "threshold:" <<
993  _SNadjacent*_CMSnoise<< " " ;
994  // Charge higher than threshold set
995  if( adjCharge >= (_SNadjacent*_CMSnoise) )
996  {
997  clsStrips[leftStrip] = new Signal(adjCharge, 0.);
998  clsStrips[leftStrip]->updateSimHitMap(adjSimHitMap);
999  // Set this charge as zero to avoid
1000  // double counting
1001  iterSMap->second[STRIPZ][leftStrip]->setCharge(0.);
1002  // And go on to the next strip
1003  ++LschMap;
1004 std::cout << " ---- OK!!" << std::endl;
1005  }
1006  else // Charge lower - stop searching
1007  {
1008  goNextStrip = false;
1009  }
1010  }*/
1011 
1012 
1013  //
1014  // Cluster vector iterators
1015  //---ClsVec::iterator iterClsVec, iterClsVec2;
1016 
1017  // Initiate variables - layerID, ladderID, sensorID, stripID, seed strip, seed charge
1018  /*short int layerID = 0;
1019  short int ladderID = 0;
1020  short int sensorID = 0;*/
1021 
1022  //int seedStrip = 0;
1023  //double seedCharge = 0;
1024 
1025  // Bunch of strips forming cluster
1026  //StripChargeMap clsStrips;
1027 
1028  //
1029  // Search complete sensor map - find seeds & their neighbouring strips
1030 /* for(SensorStripMap::iterator iterSMap=sensorMap.begin();
1031  iterSMap!=sensorMap.end(); iterSMap++)
1032  {
1033  // Save layer ID , ...
1034  int cellID = iterSMap->first;
1035 
1036  std::map<std::string,int> bfmap = _geometry->decodeCellID(cellID);
1037  layerID = bfmap["layer"];
1038  ladderID= bfmap["module"];
1039  sensorID= bfmap["sensor"]; // NOTA QUE YA NO TIENE SENTIDO
1040  // SI FRONT --> sensorID = 1
1041  // SI REAR --> sensorID = 3
1042 
1043  //
1044  // Clusters in Z
1045  for(StripChargeMap::iterator iterChMap=iterSMap->second[STRIPZ].begin();
1046  iterChMap!=iterSMap->second[STRIPZ].end(); iterChMap++)
1047  {
1048  sensorID = 3;
1049  //
1050  // Zero: Save new candidate for seed strip + MC true info
1051  int seedStrip = iterChMap->first;
1052  double seedCharge = iterChMap->second->getCharge();
1053 
1054  const SimTrackerHitMap & seedSimHitMap =
1055  iterChMap->second->getSimHitMap();
1056 
1057 
1058  if( seedCharge < (_SNseed*_CMSnoise) )
1059  {
1060  continue;
1061  }
1062  //
1063  // First: New cluster and its seed strip has been found
1064 
1065  clsStrips[seedStrip] = new Signal(seedCharge, 0.);
1066  clsStrips[seedStrip]->updateSimHitMap(seedSimHitMap);
1067 
1068  // Set this charge as zero to avoid double counting
1069  iterChMap->second->setCharge(0.);
1070 
1071  // Continue searching - find left and right neighbours
1072  // If floating strips exist --> go to every even strip
1073  int leftStrip = seedStrip -1;
1074  int rightStrip = seedStrip +1;
1075 
1076  if(_floatStripsZ)
1077  {
1078  leftStrip--;
1079  rightStrip--;
1080  }
1081 
1082  double adjCharge = 0;
1083 
1084  bool goLeft = true;
1085  bool goRight = true;
1086 
1087  // Control if left strip and right strip are in range
1088  if(leftStrip < 0)
1089  {
1090  goLeft = false;
1091  }
1092 // if(rightStrip >= _geometry->getSensorNStripsInZ(layerID))
1093  if(rightStrip >= _geometry->getSensorNStrips(layerID,sensorID)) //FIXME: >= o solo > ???
1094  {
1095  goRight = false;
1096  }
1097 
1098  // ---> Abajo codigo repetido
1099  // goleft, leftStrip , siguiente strip -
1100  // gorigth, rightStrip, siguiente strip +
1101  // ---> Puede hacerse tambien aprovechando que el mapa ordena
1102  // de menor a mayor (o al reves no me acuerdo..)
1103  // podemos pues ir de izq a derecha y de derecha a izq usando ++ -- en el iter y parando en rend() end()
1104  // Second: search for left neighbours
1105  while( goLeft && ((iterSMap->second[STRIPZ].find(leftStrip))
1106  !=iterSMap->second[STRIPZ].end()) )
1107  {
1108  adjCharge = iterSMap->second[STRIPZ][leftStrip]->getCharge();
1109 
1110  const SimTrackerHitMap & adjSimHitMap = iterSMap->second[STRIPZ][leftStrip]->getSimHitMap();
1111  // Charge higher than threshold set
1112  if( adjCharge >= (_SNadjacent*_CMSnoise) )
1113  {
1114  clsStrips[leftStrip] = new Signal(adjCharge, 0.);
1115  clsStrips[leftStrip]->updateSimHitMap(adjSimHitMap);
1116  // Set this charge as zero to avoid
1117  // double counting
1118  iterSMap->second[STRIPZ][leftStrip]->setCharge(0.);
1119  // Go on to the next strip
1120  if(_floatStripsZ)
1121  {
1122  leftStrip -= 2;
1123  }
1124  else
1125  {
1126  leftStrip -= 1;
1127  }
1128  }
1129  // Charge lower - stop searching
1130  else
1131  {
1132  goLeft = false;
1133  }
1134  }
1135 
1136  //
1137  // Third: search for right neighbours
1138  while(goRight && ((iterSMap->second[STRIPZ].find(rightStrip))
1139  !=iterSMap->second[STRIPZ].end()) )
1140  {
1141  adjCharge = iterSMap->second[STRIPZ][rightStrip]->getCharge();
1142  const SimTrackerHitMap & adjSimHitMap = iterSMap->second[STRIPZ][rightStrip]->getSimHitMap();
1143  // Charge higher than threshold set
1144  if( adjCharge >= (_SNadjacent*_CMSnoise) )
1145  {
1146  clsStrips[rightStrip] = new Signal(adjCharge, 0.);
1147  clsStrips[rightStrip]->updateSimHitMap(adjSimHitMap);
1148  // Set this charge as zero to avoid double counting
1149  iterSMap->second[STRIPZ][rightStrip]->setCharge(0.);
1150 
1151  // Go on to the next strip
1152  if(_floatStripsZ)
1153  {
1154  rightStrip += 2;
1155  }
1156  else
1157  {
1158  rightStrip += 1;
1159  }
1160  }
1161  // Charge lower - stop searching
1162  else
1163  {
1164  goRight = false;
1165  }
1166  }
1167 
1168  //
1169  // Fourth: Calculate mean position of a new cluster
1170  SimTrackerHitMap clsZSimHitMap;
1171  SimTrackerHitMap::const_iterator iterSHM;
1172 
1173  // Cluster: position, charge & size
1174  short int clsSizeZ = clsStrips.size();
1175  double clsPosZ = 0.;
1176  double clsChargeZ = 0.;
1177 
1178  double xLeftSignal = 0 ;
1179  double qLeftSignal = 0 ;
1180 
1181  double xRightSignal = 0 ;
1182  double qRightSignal = 0 ;
1183 
1184  double qIntermSignal = 0 ;
1185 
1186  for(StripChargeMap::iterator iterChMap2=clsStrips.begin();
1187  iterChMap2!=clsStrips.end(); iterChMap2++)
1188  {
1189  // Current strip ID, posZ & charge
1190  int stripID = iterChMap2->first;
1191  double stripPosZ = _geometry->getStripPosY(layerID, sensorID,stripID,0.0);
1192  double stripCharge = iterChMap2->second->getCharge();
1193 
1194  // Update info about MC particles which contributed
1195  const SimTrackerHitMap & simHitMap = iterChMap2->second->getSimHitMap();
1196 
1197  if(simHitMap.size() != 0)
1198  {
1199  for(iterSHM=simHitMap.begin();
1200  iterSHM!=simHitMap.end(); iterSHM++)
1201  {
1202  EVENT::SimTrackerHit * simHit = dynamic_cast<EVENT::SimTrackerHit *>(iterSHM->first);
1203  float weight = iterSHM->second;
1204 
1205  if(clsZSimHitMap.find(simHit)!=clsZSimHitMap.end())
1206  {
1207  clsZSimHitMap[simHit] += weight;
1208  }
1209  else
1210  {
1211  clsZSimHitMap[simHit] = weight;
1212  }
1213  }
1214  }
1215 
1216  // Get last iterator
1217  StripChargeMap::iterator iterChMapLast = clsStrips.end();
1218  -- iterChMapLast;
1219 
1220  // Get leftmost signal
1221  if (iterChMap2 == clsStrips.begin())
1222  {
1223  xLeftSignal = stripPosZ;
1224  qLeftSignal = stripCharge;
1225  }
1226 
1227  // Get rightmost signal
1228  else if (iterChMap2 == iterChMapLast)
1229  {
1230  xRightSignal = stripPosZ;
1231  qRightSignal = stripCharge;
1232  }
1233  // Get intermediate signal
1234  else
1235  {
1236  qIntermSignal += stripCharge;
1237  }
1238 
1239  // Update total charge
1240  clsChargeZ += stripCharge;
1241 
1242  }
1243 
1244  // Get average intermediate signal
1245  if (clsSizeZ > 2)
1246  {
1247  qIntermSignal /= (clsSizeZ - 2);
1248  }
1249  else
1250  {
1251  qIntermSignal = 0.;
1252  }
1253  // Analog head-tail algorithm
1254  // double geomPitchInZ = _geometry->getSensorPitchInZ(layerID);
1255  double geomPitchInZ = _geometry->getSensorPitch(layerID,sensorID,0);
1256  double readOutPitchInZ = 0.;
1257  if (_floatStripsZ)
1258  {
1259  readOutPitchInZ = 2*geomPitchInZ;
1260  }
1261  else
1262  {
1263  readOutPitchInZ = geomPitchInZ;
1264  }
1265 
1266  // qIntermSignal >= 1/eps*qRighSignal || 1/eps*qLeftSignal --> if not don't use head-tail (delta electron ...) - use epsilon ~ 1.5
1267  if ( (clsSizeZ>2) && (qLeftSignal<1.5*qIntermSignal) && (qRightSignal<1.5*qIntermSignal) )
1268  {
1269  clsPosZ = (xRightSignal + xLeftSignal)/2. + (qRightSignal - qLeftSignal)/2./qIntermSignal * readOutPitchInZ;
1270  }
1271  // COG algorithm
1272  else
1273  {
1274  clsPosZ = (xRightSignal*qRightSignal + xLeftSignal*qLeftSignal)/(qRightSignal + qLeftSignal);
1275  }
1276 
1277  //
1278  // Fifth: Correct cluster position to Lorentz shift (Bz is expected non-zero only --> no correction)
1279  clsPosZ += 0.;
1280 
1281  //
1282  // Sixth: Save information about new cluster (cluster in Z)
1283  if (clsChargeZ >= (_SNtotal*_CMSnoise))
1284  {
1285  StripCluster * pClusterZ = new StripCluster( layerID, ladderID, sensorID, Hep3Vector(_geometry->getSensorThick(layerID)/2., 0., clsPosZ),
1286  Hep3Vector(_geometry->getSensorThick(layerID)/2., 0., 0.), clsChargeZ, clsSizeZ );
1287  pClusterZ->updateSimHitMap(clsZSimHitMap);
1288 
1289  clsVecInZ.push_back(pClusterZ);
1290  }
1291 
1292  // Release memory
1293  for (StripChargeMap::iterator iterChMap2=clsStrips.begin(); iterChMap2!=clsStrips.end(); iterChMap2++)
1294  {
1295  delete iterChMap2->second;
1296  }
1297  // Clear content
1298  clsStrips.clear();
1299 
1300  } // For clusters in Z
1301 
1302  //
1303  // Clusters in R-Phi + 3D clusters utilizing Z (use Z position of "clusters in Z" and create 3D clusters from them)
1304  for(StripChargeMap::iterator iterChMap=iterSMap->second[STRIPRPHI].begin(); iterChMap!=iterSMap->second[STRIPRPHI].end(); iterChMap++) {
1305 
1306  sensorID = 1;
1307 
1308  //
1309  // Zero: Save new candidate for seed strip + MC true info
1310  int seedStrip = iterChMap->first;
1311  double seedCharge = iterChMap->second->getCharge();
1312 
1313  const SimTrackerHitMap & seedSimHitMap = iterChMap->second->getSimHitMap();
1314 
1315  //
1316  // First: New cluster and its seed strip has been found
1317  if ( seedCharge >= (_SNseed*_CMSnoise) ) {
1318 
1319  clsStrips[seedStrip] = new Signal(seedCharge, 0.);
1320  clsStrips[seedStrip]->updateSimHitMap(seedSimHitMap);
1321 
1322  // Set this charge as zero to avoid double counting
1323  iterChMap->second->setCharge(0.);
1324 
1325  // Continue searching - find left and right neighbours
1326  bool goLeft = true;
1327  bool goRight = true;
1328 
1329  // If floating strips exist --> go to every even strip
1330  int leftStrip = -1;
1331  int rightStrip = -1;
1332 
1333  if (_floatStripsRPhi) {
1334 
1335  leftStrip = seedStrip - 2;
1336  rightStrip = seedStrip + 2;
1337  }
1338  else {
1339 
1340  leftStrip = seedStrip - 1;
1341  rightStrip = seedStrip + 1;
1342  }
1343 
1344  double adjCharge = 0;
1345 
1346  // Control if left strip and right strip are in range
1347  if (leftStrip < 0) goLeft = false;
1348 // if (rightStrip >= _geometry->getSensorNStripsInRPhi(layerID)) goRight = false;
1349  if (rightStrip >= _geometry->getSensorNStrips(layerID,sensorID)) goRight = false;
1350 
1351  //
1352  // Second: Search for left neighbours
1353  while ( goLeft && ((iterSMap->second[STRIPRPHI].find(leftStrip))!=iterSMap->second[STRIPRPHI].end()) ) {
1354 
1355  adjCharge = iterSMap->second[STRIPRPHI][leftStrip]->getCharge();
1356 
1357  const SimTrackerHitMap & adjSimHitMap = iterSMap->second[STRIPRPHI][leftStrip]->getSimHitMap();
1358 
1359  // Charge higher than threshold set
1360  if ( adjCharge >= (_SNadjacent*_CMSnoise) ) {
1361 
1362  clsStrips[leftStrip] = new Signal(adjCharge, 0.);
1363  clsStrips[leftStrip]->updateSimHitMap(adjSimHitMap);
1364 
1365  // Set this charge as zero to avoid double counting
1366  iterSMap->second[STRIPRPHI][leftStrip]->setCharge(0.);
1367 
1368  // Go on to the next strip
1369  if (_floatStripsRPhi) leftStrip -= 2;
1370  else leftStrip -= 1;
1371  }
1372  // Charge lower - stop searching
1373  else goLeft = false;
1374  }
1375 
1376  //
1377  // Third: Search for right neighbours
1378  while (goRight && ((iterSMap->second[STRIPRPHI].find(rightStrip))!=iterSMap->second[STRIPRPHI].end()) ) {
1379 
1380  adjCharge = iterSMap->second[STRIPRPHI][rightStrip]->getCharge();
1381 
1382  const SimTrackerHitMap & adjSimHitMap = iterSMap->second[STRIPRPHI][rightStrip]->getSimHitMap();
1383 
1384  // Charge higher than threshold set
1385  if ( adjCharge >= (_SNadjacent*_CMSnoise) ) {
1386 
1387  clsStrips[rightStrip] = new Signal(adjCharge, 0.);
1388  clsStrips[rightStrip]->updateSimHitMap(adjSimHitMap);
1389 
1390  // Set this charge as zero to avoid double counting
1391  iterSMap->second[STRIPRPHI][rightStrip]->setCharge(0.);
1392 
1393  // Go on to the next strip
1394  if (_floatStripsRPhi) rightStrip += 2;
1395  else rightStrip += 1;
1396  }
1397  // Charge lower - stop searching
1398  else goRight = false;
1399  }
1400 
1401  //
1402  // Fourth: Go through all Z clusters and calculate mean position in R-Phi and save them together as a 3D cluster
1403  for (iterClsVec=clsVecInZ.begin(); iterClsVec!=clsVecInZ.end(); iterClsVec++) {
1404 
1405  // Get and save info about cluster in Z
1406  StripCluster * pClusterZ = *iterClsVec;
1407 
1408  const SimTrackerHitMap & clsZSimHitMap = pClusterZ->getSimHitMap();
1409 
1410  // Cluster: position, charge & size
1411  short int clsSizeRPhi = clsStrips.size();
1412  double clsPosRPhi = 0.;
1413  double clsChargeRPhi = 0.;
1414 
1415  double xLeftSignal = 0 ;
1416  double qLeftSignal = 0 ;
1417 
1418  double xRightSignal = 0 ;
1419  double qRightSignal = 0 ;
1420 
1421  double qIntermSignal = 0 ;
1422 
1423  SimTrackerHitMap clsRPhiSimHitMap;
1424  SimTrackerHitMap::const_iterator iterSHM;
1425 
1426  // Find clusters in R-Phi
1427  for (StripChargeMap::iterator iterChMap2=clsStrips.begin(); iterChMap2!=clsStrips.end(); iterChMap2++) {
1428 
1429  // Current strip ID, posRPhi & charge
1430  int stripID = iterChMap2->first;
1431  //double stripPosRPhi = _geometry->getStripPosInRPhi(layerID, stripID, pClusterZ->getPosZ());
1432  double stripPosRPhi = _geometry->getStripPosY(layerID, sensorID, stripID,pClusterZ->getPosZ());
1433  double stripCharge = iterChMap2->second->getCharge();
1434 
1435  // Update info about MC particles which contributed
1436  const SimTrackerHitMap & simHitMap = iterChMap2->second->getSimHitMap();
1437 
1438  if (simHitMap.size() != 0) {
1439 
1440  for (iterSHM=simHitMap.begin(); iterSHM!=simHitMap.end(); iterSHM++) {
1441 
1442  EVENT::SimTrackerHit * simHit = dynamic_cast<EVENT::SimTrackerHit *>(iterSHM->first);
1443  float weight = iterSHM->second;
1444 
1445  if (clsRPhiSimHitMap.find(simHit)!=clsRPhiSimHitMap.end()) clsRPhiSimHitMap[simHit] += weight;
1446  else clsRPhiSimHitMap[simHit] = weight;
1447  }
1448  }
1449 
1450  // Get last iterator
1451  StripChargeMap::iterator iterChMapLast = clsStrips.end();
1452  -- iterChMapLast;
1453 
1454  // Get leftmost signal
1455  if (iterChMap2 == clsStrips.begin()) {
1456 
1457  xLeftSignal = stripPosRPhi;
1458  qLeftSignal = stripCharge;
1459  }
1460 
1461  // Get rightmost signal
1462  else if (iterChMap2 == iterChMapLast) {
1463 
1464  xRightSignal = stripPosRPhi;
1465  qRightSignal = stripCharge;
1466  }
1467  // Get intermediate signal
1468  else qIntermSignal += stripCharge;
1469 
1470  // Update total charge
1471  clsChargeRPhi += stripCharge;
1472 
1473  }
1474 
1475  // Get average intermediate signal
1476  if (clsSizeRPhi > 2) qIntermSignal /= (clsSizeRPhi - 2);
1477  else qIntermSignal = 0.;
1478 
1479  // Analog head-tail algorithm
1480  //double geomPitchInRPhi = _geometry->getSensorPitchInRPhi(layerID, pClusterZ->getPosZ());
1481  double geomPitchInRPhi = _geometry->getSensorPitch(layerID,sensorID, pClusterZ->getPosZ());
1482  double readOutPitchInRPhi = 0.;
1483  if (_floatStripsRPhi) readOutPitchInRPhi = 2*geomPitchInRPhi;
1484  else readOutPitchInRPhi = geomPitchInRPhi;
1485 
1486  // qIntermSignal >= 1/eps*qRighSignal || 1/eps*qLeftSignal --> if not don't use head-tail (delta electron ...) - use epsilon ~ 1.5
1487  if ( (clsSizeRPhi>2) && (qLeftSignal<1.5*qIntermSignal) && (qRightSignal<1.5*qIntermSignal) ) {
1488 
1489  clsPosRPhi = (xRightSignal + xLeftSignal)/2. + (qRightSignal - qLeftSignal)/2./qIntermSignal * readOutPitchInRPhi;
1490  }
1491  // COG algorithm
1492  else clsPosRPhi = (xRightSignal*qRightSignal + xLeftSignal*qLeftSignal)/(qRightSignal + qLeftSignal);
1493 
1494 
1495  // Fifth: Correct mean position to Lorentz shift (Bz is expected non-zero only; correct to theta angle of each sensor)
1496  clsPosRPhi += _TanOfAvgELorentzShift * _geometry->getSensorThick(layerID)/2. * cos(_geometry->getLadderTheta(layerID));
1497 
1498  //
1499  // Sixth: Save information about new cluster (cluster in RPhi), if cluster charge higher than threshold set
1500  if (clsChargeRPhi >= (_SNtotal*_CMSnoise)) {
1501 
1502  StripCluster * pClusterRPhi = new StripCluster( layerID, ladderID, sensorID, Hep3Vector(_geometry->getSensorThick(layerID)/2., clsPosRPhi, 0.),
1503  Hep3Vector(_geometry->getSensorThick(layerID)/2., 0., 0.), clsChargeRPhi, clsSizeRPhi );
1504 
1505  pClusterRPhi->updateSimHitMap(clsRPhiSimHitMap);
1506 
1507  clsVecInRPhi.push_back(pClusterRPhi);
1508 
1509  // Create final 3D cluster
1510  Hep3Vector position( _geometry->getSensorThick(layerID)/2., pClusterRPhi->getPosY() , pClusterZ->getPosZ());
1511  Hep3Vector posSigma( _geometry->getSensorThick(layerID)/2., pClusterRPhi->getPosSigmaY(), pClusterZ->getPosSigmaZ());
1512 
1513  double totalCharge = (pClusterRPhi->getCharge() + pClusterZ->getCharge() )/2.;
1514 
1515  StripCluster * pCluster3D = new StripCluster( layerID, ladderID, sensorID, position, posSigma, totalCharge, 0);
1516 
1517  SimTrackerHitMap cls3DSimHitMap = clsRPhiSimHitMap;
1518 
1519  // Update MC true info for 3D
1520  if (cls3DSimHitMap.size() != 0) {
1521 
1522  for (iterSHM=clsZSimHitMap.begin(); iterSHM!=clsZSimHitMap.end(); iterSHM++) {
1523 
1524  EVENT::SimTrackerHit * simHit = dynamic_cast<EVENT::SimTrackerHit *>(iterSHM->first);
1525  float weight = iterSHM->second;
1526 
1527  if (cls3DSimHitMap.find(simHit)!=cls3DSimHitMap.end()) cls3DSimHitMap[simHit] += weight;
1528  else cls3DSimHitMap[simHit] = weight;
1529  }
1530 
1531  pCluster3D->updateSimHitMap(cls3DSimHitMap);
1532  }
1533 
1534  //
1535  // if defined ROOT-OUTPUT, save info
1536 #ifdef ROOT_OUTPUT
1537 
1538  // Set layer, ladder, sensor ID
1539  _rootLayerID = layerID;
1540  _rootLadderID = ladderID;
1541  _rootSensorID = sensorID;
1542 
1543  // Set reconstructed positions
1544  _rootRecRPhi = pClusterRPhi->getPosY() / mm;
1545  _rootRecZ = pClusterZ->getPosZ() / mm;
1546 
1547  // Set cluster sizes
1548  _rootClsSizeRPhi = pClusterRPhi->getSize();
1549  _rootClsSizeZ = pClusterZ->getSize();
1550 
1551  // Set simulated position in RPhi & MC particle
1552 
1553  // Find SimTrackerHit with highest weight
1554  EVENT::SimTrackerHit * simHit = 0;
1555  float weight = 0;
1556 
1557  const SimTrackerHitMap & simHitMapRPhi = pClusterRPhi->getSimHitMap();
1558 
1559  for (iterSHM=simHitMapRPhi.begin(); iterSHM!=simHitMapRPhi.end(); iterSHM++) {
1560 
1561  // Find contribution with highest weight
1562  if ( (iterSHM->first!=0) && ((iterSHM->second)>weight) ) {
1563 
1564  simHit = iterSHM->first;
1565  weight = iterSHM->second;
1566  }
1567  }
1568 
1569  // SimHit global position
1570  Hep3Vector simPosGlob;
1571  Hep3Vector simPosLoc;
1572 
1573  simPosGlob.setX(simHit->getPosition()[0]*mm);
1574  simPosGlob.setY(simHit->getPosition()[1]*mm);
1575  simPosGlob.setZ(simHit->getPosition()[2]*mm);
1576  simPosLoc = _geometry->transformPointToLocal(layerID, ladderID, sensorID, simPosGlob);
1577 
1578  _rootMCPDGRPhi = simHit->getMCParticle()->getPDG();
1579  _rootSimRPhi = simPosLoc.getY() / mm;
1580 
1581  // Set simulated position in Z & MC particle
1582 
1583  // Find SimTrackerHit with highest weight
1584  const SimTrackerHitMap & simHitMapZ = pClusterZ->getSimHitMap();
1585 
1586  for (iterSHM=simHitMapZ.begin(); iterSHM!=simHitMapZ.end(); iterSHM++) {
1587 
1588  // Find contribution with highest weight
1589  if ( (iterSHM->first!=0) && ((iterSHM->second)>weight) ) {
1590 
1591  simHit = iterSHM->first;
1592  weight = iterSHM->second;
1593  }
1594  }
1595 
1596  // SimHit global position
1597  simPosGlob.setX(simHit->getPosition()[0]*mm);
1598  simPosGlob.setY(simHit->getPosition()[1]*mm);
1599  simPosGlob.setZ(simHit->getPosition()[2]*mm);
1600  simPosLoc = _geometry->transformPointToLocal(layerID, ladderID, sensorID, simPosGlob);
1601 
1602  _rootMCPDGZ = simHit->getMCParticle()->getPDG();
1603  _rootSimZ = simPosLoc.getZ() / mm;
1604 
1605  // Fill the tree
1606  _rootFile->cd("");
1607  _rootTree->Fill();
1608 
1609 #endif
1610 
1611  clsVec.push_back(pCluster3D);
1612  }
1613 
1614  } // Go through all Z clusters
1615 
1616  // Release memory
1617  for (StripChargeMap::iterator iterChMap2=clsStrips.begin(); iterChMap2!=clsStrips.end(); iterChMap2++) delete iterChMap2->second;
1618 
1619  // Clear content
1620  clsStrips.clear();
1621 
1622  } // If found new cluster
1623 
1624  } // For clusters in R-Phi
1625 
1626  //
1627  // Release memory
1628  for (iterClsVec=clsVecInRPhi.begin(); iterClsVec!=clsVecInRPhi.end(); iterClsVec++) delete *iterClsVec;
1629  for (iterClsVec=clsVecInZ.begin(); iterClsVec!=clsVecInZ.end(); iterClsVec++) delete *iterClsVec;
1630 
1631  // Clear
1632  clsVecInRPhi.clear();
1633  clsVecInZ.clear();
1634 
1635  } // For sensor map*/
1636 
1637 // OTHER METHODS
1638 
1639 //
1640 // Calculated and stored the hits adjacents
1641 // Template for the reverse_iterator (leftStrips) and iterator (rightStrips)
1642 //
1643 
1644 template<class It>
1645 StripChargeMap & SiStripClus::storeHitsAdjacents( StripChargeMap & clsStrips,
1646  It schMap, const It & endIt, StripChargeMap & currentMap)
1647 {
1648  // Map ordered from lower to higher
1649  double adjCharge = 0.0;
1650  bool goNextStrip = true;
1651  while( goNextStrip && schMap != endIt )
1652  {
1653  adjCharge = schMap->second->getCharge();
1654  const SimTrackerHitMap & adjSimHitMap = schMap->second->getSimHitMap();
1655  // Charge higher than threshold set
1656  if( adjCharge >= (_SNadjacent*_CMSnoise) )
1657  {
1658  clsStrips[schMap->first] = new Signal(adjCharge, 0.);
1659  clsStrips[schMap->first]->updateSimHitMap(adjSimHitMap);
1660  // Set this charge as zero to avoid
1661  // double counting
1662  currentMap[schMap->first]->setCharge(0.);
1663  // And go on to the next strip
1664  ++schMap;
1665  }
1666  else // Charge lower - stop searching
1667  {
1668  goNextStrip = false;
1669  }
1670  }
1671 
1672  return clsStrips;
1673 }
1674 
1675 
1676 //
1677 // Method calculating hits from given clusters
1678 //
1679 // FIXME: Returning colOfTrkHits, clsVec reference?? --> the clsVec.clear
1680 void SiStripClus::calcHits(ClsVec & clsVec, IMPL::LCCollectionVec * colOfTrkHits)
1681 {
1682  // FIXME: Why passing clsVec by reference?
1683  // return the LCCol
1684  // Cluster - position, covariance matrix (3x3 = 6 parameters = lower triangle matrix)
1685  ClsVec::iterator iterClsVec;
1686 
1687  short int layerID;
1688  short int ladderID;
1689  short int sensorID;
1690 
1691  Hep3Vector position;
1692  Hep3Vector posSigma;
1693 
1694  double totalCharge;
1695 
1696  // Print sensor info
1697  streamlog_out(MESSAGE2) << std::endl
1698  << " Total number of reconstructed hits: "
1699  << clsVec.size()
1700  << " hit(s)"
1701  << std::endl;
1702 
1703  // Set collection flag - cellID 1 will be stored,
1704  // and LCrelations
1705  LCFlagImpl flag1(0);
1706  flag1.setBit(LCIO::RTHPBIT_ID1);
1707  colOfTrkHits->setFlag(flag1.getFlag());
1708 
1709  // CODIFICATION --->FIXME: METHOD IN GEAR (Centralizing...)
1710  CellIDEncoder<TrackerHitPlaneImpl> cellEnc(
1711  LCTrackerCellID::encoding_string()+",stripFront:11,stripRear:11",
1712  colOfTrkHits);
1713 
1714  // Go through all clusters
1715  for(iterClsVec=clsVec.begin(); iterClsVec!=clsVec.end(); iterClsVec++)
1716  {
1717  // Get cluster info
1718  StripCluster * pCluster = *iterClsVec;
1719 
1720  layerID = pCluster->getLayerID();
1721  ladderID = pCluster->getLadderID();
1722  sensorID = pCluster->getSensorID();
1723 
1724  position = pCluster->get3Position();
1725  posSigma = pCluster->get3PosSigma();
1726 
1727  totalCharge = pCluster->getCharge();
1728 
1729  CLHEP::Hep3Vector localPos = position;
1730 
1731  // Transform hit and covariance to the global ref. system
1732  position = _geometry->transformPointToGlobal(layerID, ladderID,
1733  sensorID, position);
1734  // Reasigning the z-position to be placed in the middle between the
1735  // front and rear sensor
1736  const int sign = abs(_geometry->getLayerRealID(layerID))
1737  /_geometry->getLayerRealID(layerID);
1738  const double offsetZ = _geometry->getLadderThick(layerID)/2.
1739  +_geometry->getSensorThick(layerID)/2.;
1740  position.setZ(position.getZ()+sign*offsetZ);
1741 
1742  // Save into LCIO native variables and in appropriate units
1743  double posLCIO[3] = { position.getX()/mm, position.getY()/mm,
1744  position.getZ()/mm };
1745 
1746  // Calculate resolution
1747  float * covLCIO = calcResolution(layerID, position.theta()/pi * 180.,
1748  localPos.getZ());
1749 
1750  // Create new Tracker hit
1751  TrackerHitPlaneImpl * trkHit = new TrackerHitPlaneImpl();
1752 
1753  // Set hit type: SVD type is 201 + layer ID ??
1754  trkHit->setType(layerID + 201);
1755  trkHit->setPosition(posLCIO);
1756  // U and V error, not using covMatrix anymore
1757  // trkHit->setCovMatrix(covLCIO);
1758  trkHit->setdU( covLCIO[2] );
1759  trkHit->setdV( covLCIO[5] );
1760  trkHit->setEDep(totalCharge / e); // in electrons
1761  //FIXME: Still missing EDepError
1762  trkHit->setTime(pCluster->getTime()); //FIXME: NOT IMPLEMENTED
1763 
1764  // The unit vector of the local reference frame of the front sensor
1765  // u: rphi direction, v: radial direction
1766  CLHEP::Hep3Vector Uvector = _geometry->transformVecToGlobal(layerID,
1767  ladderID,1,CLHEP::Hep3Vector(0.0,1.0,0.0));
1768  float Uf[2] = { Uvector.getX(), Uvector.getY() };
1769  CLHEP::Hep3Vector Vvector = _geometry->transformVecToGlobal(layerID,
1770  ladderID,1,CLHEP::Hep3Vector(0.0,0.0,1.0));
1771  float Vf[2] = { Vvector.getX(), Vvector.getY() };
1772  trkHit->setU( Uf );
1773  trkHit->setV( Vf );
1774 
1775  // Set simTrkHits which contributed & find hit with highest weight
1776  const SimTrackerHitMap & simHitMap = pCluster->getSimHitMap();
1777  //float weightSum = pCluster->getSimHitWeightSum();
1778  SimTrackerHit * simTrkHit = 0;
1779  float weight = 0;
1780 
1781  for(SimTrackerHitMap::const_iterator iterSHM=simHitMap.begin();
1782  iterSHM!=simHitMap.end(); iterSHM++)
1783  {
1784  // Don't save "noise hits", i.e. zero pointers
1785  // Find contribution with highest weight
1786  if( (iterSHM->first!=0) && ((iterSHM->second)>weight) )
1787  {
1788  simTrkHit = iterSHM->first;
1789  weight = iterSHM->second;
1790  }
1791  }
1792 
1793  int realLayer = _geometry->getLayerRealID(layerID);
1794  // Codification FIXME: Call FTD method?
1795  cellEnc["subdet"]=ILDDetID::FTD;
1796  cellEnc["side"]=abs(realLayer)/realLayer;
1797  cellEnc["layer"]=abs(realLayer);
1798  cellEnc["module"]=ladderID+1;
1799  cellEnc["sensor"]=1;
1800  cellEnc["stripFront"]=pCluster->getStripFront();
1801  cellEnc["stripRear"]=pCluster->getStripRear();
1802 
1803  cellEnc.setCellID(trkHit);
1804 
1805 
1806  // Save only hit with highest weight
1807  trkHit->rawHits().push_back(simTrkHit);
1808 
1809  // Save the hit to the collection
1810  colOfTrkHits->addElement(trkHit);
1811 
1812  // Print infor
1813  printHitInfo(pCluster);
1814 
1815  // Release memory
1816  delete pCluster;
1817 
1818  } // For
1819 
1820  streamlog_out(MESSAGE2) << std::endl;
1821 
1822  // Clear content
1823  clsVec.clear();
1824 }
1825 
1826 //
1827 // Method calculating hit resolution, i.e. covariance matrix
1828 // The posZ is in local frame coordinate system
1829 float * SiStripClus::calcResolution(const int & layerID,const double & hitTheta,
1830  const double & posZ)
1831 {
1832  static float covMatrix[6] = { 0., 0., 0.,
1833  0., 0., 0. };
1834  // Assuming resolution as pitch/2
1835  // we need the zposlocal in the local reference system
1836  const double halfPitch = _geometry->getSensorPitch(layerID,1,posZ);
1837 
1838  const double sigmax = halfPitch/(2.*cos(hitTheta));
1839  const double sigmay = halfPitch*sin(hitTheta)/2.;
1840 
1841  // Set covariance matrix in appropriate units*/
1842  // Assuming no correlations between the front and rear sensor
1843  covMatrix[2] = sigmax/mm * sigmax/mm; // Local frame Y (u-vector)
1844  covMatrix[5] = sigmay/mm * sigmay/mm; // Local frame Z (v-vector)
1845  covMatrix[0] = (_geometry->getLadderThick(layerID) // Local frame x
1846  +_geometry->getSensorThick(layerID))/2.0;
1847 
1848  return covMatrix;
1849  // Define array of theta angle
1850  /*static float theta[14] = {20., 30., 40., 50., 60., 70., 80., 90.,
1851  100., 110., 120., 130., 140., 150.};
1852 
1853  // Find correct bin for theta (-1 for < 20deg; 14 for > 150deg)
1854  int iTheta = -1;
1855  for(int i=0; i<14; i++)
1856  {
1857  if (hitTheta>=theta[i])
1858  {
1859  iTheta++;
1860  }
1861  }
1862 
1863  // Define tracker hit resolution
1864  float resInZ = 0.;
1865  float resInRPhi = 0.;
1866 
1867  //
1868  // Calculate resolution - interpolate between given theta angles
1869 
1870  // Theta is less than 20 degrees
1871  if (iTheta == -1)
1872  {
1873  // SVD - First layer
1874  if (layerID == 2)
1875  {
1876  resInZ = _resSVDFirstInZ[0];
1877  resInRPhi = _resSVDFirstInRPhi[0];
1878  }
1879 
1880  // SVD - Other layers
1881  else
1882  {
1883  resInZ = _resSVDOtherInZ[0];
1884  resInRPhi = _resSVDOtherInRPhi[0];
1885  }
1886  }
1887  // Theta is higher than 150 degrees
1888  else if (iTheta == 14)
1889  {
1890  // SVD - First layer
1891  if (layerID == 2)
1892  {
1893  resInZ = _resSVDFirstInZ[13];
1894  resInRPhi = _resSVDFirstInRPhi[13];
1895  }
1896 
1897  // SVD - Other layers
1898  else
1899  {
1900  resInZ = _resSVDOtherInZ[13];
1901  resInRPhi = _resSVDOtherInRPhi[13];
1902  }
1903  }
1904 
1905  // Interpolate between 20 and 150 degrees
1906  else
1907  {
1908  // SVD - First layer
1909  if (layerID == 2)
1910  {
1911  resInZ = _resSVDFirstInZ[iTheta] +
1912  (_resSVDFirstInZ[iTheta+1] - _resSVDFirstInZ[iTheta] )
1913  /(theta[iTheta+1] - theta[iTheta])*(hitTheta - theta[iTheta]);
1914  resInRPhi = _resSVDFirstInRPhi[iTheta] +
1915  (_resSVDFirstInRPhi[iTheta+1] - _resSVDFirstInRPhi[iTheta])
1916  /(theta[iTheta+1] - theta[iTheta])*(hitTheta - theta[iTheta]);
1917  }
1918 
1919  // SVD - Other layers
1920  else
1921  {
1922  resInZ = _resSVDOtherInZ[iTheta] +
1923  (_resSVDOtherInZ[iTheta+1] - _resSVDOtherInZ[iTheta] )
1924  /(theta[iTheta+1] - theta[iTheta])*(hitTheta - theta[iTheta]);
1925  resInRPhi = _resSVDOtherInRPhi[iTheta] +
1926  (_resSVDOtherInRPhi[iTheta+1] - _resSVDOtherInRPhi[iTheta])
1927  /(theta[iTheta+1] - theta[iTheta])*(hitTheta - theta[iTheta]);
1928  }
1929  }
1930 
1931  // Set covariance matrix in appropriate units
1932  covMatrix[2] = resInRPhi/mm * resInRPhi/mm;
1933  covMatrix[5] = resInZ/mm * resInZ/mm;
1934 
1935  return covMatrix;*/
1936 }
1937 
1938 //
1939 // Method to save the signal
1940 //
1941 // FTD stuff --> Need to change the sensorMap in order to store in the
1942 // same StripCharge map the two opposites sensors of the same disk-petal
1943 // as in FTD the Hits are built using two Single side sensors of the same
1944 // disk-petal
1945 //
1946 // FIXME: Returns the SensorStripMap
1947 void SiStripClus::updateMap(TrackerPulseImpl * pulse,
1948  SensorStripMap & sensorMap )
1949 {
1950  // CellID0 encodes layerID, ladderID and sensorID;
1951  // cellID1 encodes strip
1952  int cellID0 = pulse->getCellID0();
1953  // the new cellID0 is as the old one but extracted the sensor
1954  // information (cellID["sensor"] = 0
1955  cellID0 = _geometry->cellID0withSensor0(cellID0);
1956 
1957  int cellID1 = pulse->getCellID1();
1958 
1959  double charge = pulse->getCharge()*fC;
1960 
1961  // Decode stripType and stripID
1962  std::pair<StripType,int> tpIdpair = _geometry->decodeStripID(cellID1);
1963 
1964  const StripType stripType = tpIdpair.first;
1965  const int stripID = tpIdpair.second;
1966 
1967  // Controlling some errors
1968  if( stripType != STRIPFRONT && stripType != STRIPREAR )
1969  {
1970  streamlog_out(ERROR)
1971  << "cellID1 - problem to identify if strips in FRONT or REAR!!!"
1972  << std::endl;
1973  exit(0);
1974  }
1975 
1976  Signal *signal = 0;
1977 
1978  // Find if sensor already saved in map
1979  SensorStripMap::iterator iterSMap = sensorMap.find(cellID0);
1980  // If not, create map for the sensor (cellID0)
1981  if(iterSMap == sensorMap.end())
1982  {
1983  // Create map
1984  sensorMap[cellID0] = new StripChargeMap[2];
1985  // And update the iterator
1986  iterSMap = sensorMap.find(cellID0);
1987  }
1988 
1989  // Find if strip already saved in map
1990  StripChargeMap::iterator iterChMap = iterSMap->second[stripType].find(stripID);
1991  // IF not, create entry for strip, otherwise update
1992  if( iterChMap == iterSMap->second[stripType].end() )
1993  {
1994  signal = new Signal(charge,0.0);
1995  iterSMap->second[stripType][stripID] = signal;
1996  // And update the iterator
1997  iterChMap = iterSMap->second[stripType].find(stripID);
1998  }
1999  else
2000  {
2001  iterChMap->second->updateCharge(charge);
2002  }
2003 
2004  // Get MCParticles which contributed and update MCParticles
2005  if(_navigatorPls != NULL)
2006  {
2007  LCObjectVec lcObjVec = _navigatorPls->getRelatedToObjects(pulse);
2008  FloatVec floatVec = _navigatorPls->getRelatedToWeights(pulse);
2009 
2010  LCObjectVec::iterator iterLCObjVec = lcObjVec.begin();
2011  FloatVec::iterator iterFloatVec = floatVec.begin();
2012 
2013  for(iterLCObjVec=lcObjVec.begin(); iterLCObjVec!=lcObjVec.end();
2014  iterLCObjVec++, iterFloatVec++)
2015  {
2016  EVENT::SimTrackerHit * simHit = dynamic_cast<EVENT::SimTrackerHit *>(*iterLCObjVec);
2017  float weight = *iterFloatVec;
2018 
2019  iterChMap->second->updateSimHitMap(simHit, weight);
2020  }
2021  }
2022  //FIXME CONTROL DE ERRORES
2023  //return true;
2024 }
2025 
2026 void SiStripClus::releaseMap( SensorStripMap & sensorMap )
2027 {
2028  // Release memory
2029  for(SensorStripMap::iterator iterSMap=sensorMap.begin();
2030  iterSMap!=sensorMap.end(); iterSMap++)
2031  {
2032  //Array contents
2033  // Strips in Front
2034  for(StripChargeMap::iterator iterChMap=iterSMap->second[STRIPFRONT].begin();
2035  iterChMap!=iterSMap->second[STRIPFRONT].end(); iterChMap++)
2036  {
2037  delete iterChMap->second;
2038  }
2039  // Strips in Rear
2040  for(StripChargeMap::iterator iterChMap=iterSMap->second[STRIPREAR].begin();
2041  iterChMap!=iterSMap->second[STRIPREAR].end(); iterChMap++)
2042  {
2043  delete iterChMap->second;
2044  }
2045 
2046  // Release memory for array
2047  delete [] iterSMap->second;
2048  }
2049 
2050  // Clearing
2051  sensorMap.clear();
2052 }
2053 
2054 // PRINT METHODS
2055 
2056 //
2057 // Method printing processor parameters
2058 //
2059 void SiStripClus::printProcessorParams() const
2060 {
2061  streamlog_out(MESSAGE3) << std::endl
2062  << " "
2063  << DUNDERL
2064  << DBLUE
2065  << "SiStripClus parameters:"
2066  << ENDCOLOR
2067  << " "
2068  << std::endl << std::endl;
2069 
2070  streamlog_out(MESSAGE3) << std::setiosflags(std::ios::fixed | std::ios::internal )
2071  << std::setprecision(2)
2072  << " CMS noise [fC]: " << std::setw(4) << _CMSnoise/fC << std::endl
2073  << " " << std::endl
2074  << std::setprecision(1)
2075  << " S/N cut for seed strips: " << std::setw(3) << _SNseed << std::endl
2076  << " S/N cut for adjacent strips: " << std::setw(3) << _SNadjacent << std::endl
2077  << " S/N cut for total charge: " << std::setw(3) << _SNtotal << std::endl
2078  << std::resetiosflags(std::ios::showpos)
2079  << std::setprecision(0)
2080  << std::endl;
2081 
2082  if (_floatStripsRPhi)
2083  streamlog_out(MESSAGE3) << " Read-out pitch in R-Phi is 2x geom. pitch." << std::endl;
2084  if (_floatStripsZ)
2085  streamlog_out(MESSAGE3) << " Read-out pitch in Z is 2x geom. pitch." << std::endl << std::endl;
2086 
2087 }
2088 
2089 
2090 //
2091 // Method printing hit info
2092 //
2093 void SiStripClus::printHitInfo(const StripCluster * pCluster) const
2094 {
2095  short int layerID = pCluster->getLayerID();
2096  short int ladderID = pCluster->getLadderID();
2097  short int sensorID = pCluster->getSensorID();
2098 
2099  Hep3Vector position = pCluster->get3Position();
2100 
2101  // Print sensor info
2102  streamlog_out(MESSAGE2) << " Layer "
2103  << _geometry->getLayerRealID(layerID) << " "
2104  << "Ladder "
2105  << ladderID << " "
2106  << "Sensor "
2107  << sensorID << std::endl;
2108  streamlog_out(MESSAGE2) << " Hit in local ref. system: "
2109  << std::fixed
2110  << std::setprecision(3)
2111  << std::setiosflags(std::ios::showpos)
2112  << "PosX [mm]: " << position.getX()/mm << " "
2113  << "PosY [mm]: " << position.getY()/mm << " "
2114  << "PosZ [mm]: " << position.getZ()/mm
2115  << std::setprecision(0)
2116  << std::resetiosflags(std::ios::internal | std::ios::showpos)
2117  << std::endl;
2118 
2119  position = _geometry->transformPointToGlobal(layerID, ladderID, sensorID, position);
2120 
2121  streamlog_out(MESSAGE2) << " Hit in global ref. system: "
2122  << std::fixed
2123  << std::setprecision(3)
2124  << std::setiosflags(std::ios::showpos)
2125  << "PosX [mm]: " << position.getX()/mm << " "
2126  << "PosY [mm]: " << position.getY()/mm << " "
2127  << "PosZ [mm]: " << position.getZ()/mm
2128  << std::setprecision(0)
2129  << std::resetiosflags(std::ios::internal | std::ios::showpos)
2130  << std::endl;
2131 
2132 }
2133 
2134 
2135 } // Namespace
sistrip::SiStripClus::end
virtual void end()
Method called after all data processing.
Definition: SiStripClus.cc:473
sistrip::StripCluster::getTime
double getTime() const
Get time when the cluster has been created by a particle.
Definition: StripCluster.h:138
sistrip::mm
static const float mm
Definition: PhysicalConstants.h:37
sistrip::SiStripGeom::getSensorThick
virtual double getSensorThick(short int layerID) const
Get sensor thickness.
Definition: SiStripGeom.cc:421
sistrip::SiStripGeomBuilder::Build
static SiStripGeom * Build(const std::string &detector, const double &pitchFront, const double &pitchRear)
Definition: SiStripGeomBuilder.cc:28
sistrip::SiStripClus::_pitchRear
double _pitchRear
Pitch in the middle of the rear sensor.
Definition: SiStripClus.h:160
sistrip::SiStripGeom::transformPointToLocal
virtual CLHEP::Hep3Vector transformPointToLocal(short int layerID, short int ladderID, short int sensorID, const CLHEP::Hep3Vector &point)=0
Transform given point from global ref. system to local ref. system (sensor)
sistrip::StripCluster::updateSimHitMap
void updateSimHitMap(SimTrackerHitMap simHitMap)
Update MC truth information about SimTrackerHits, which contributed.
Definition: StripCluster.cc:46
sistrip::SiStripClus::calcResolution
float * calcResolution(const int &layerID, const double &hitTheta, const double &posZ)
Method calculating hit resolution, i.e. covariance matrix.
Definition: SiStripClus.cc:1829
sistrip::SiStripClus::_SNtotal
float _SNtotal
Signal to noise ratio cut for total cluster.
Definition: SiStripClus.h:138
sistrip::um
static const float um
Definition: PhysicalConstants.h:38
DBLUE
#define DBLUE
Definition: Colours.h:16
sistrip::SiStripClus::updateMap
void updateMap(TrackerPulseImpl *pulse, SensorStripMap &sensorMap)
Method to update and store the Sensor strip map.
Definition: SiStripClus.cc:1947
sistrip::StripCluster::getSensorID
short int getSensorID() const
Get cluster sensor ID.
Definition: StripCluster.h:108
sistrip::SiStripClus::check
virtual void check(LCEvent *event)
Method called after each event - used for data checking.
Definition: SiStripClus.cc:466
STRIPFRONT
#define STRIPFRONT
Definition: SiStripDigi.h:123
sistrip::SiStripClus::findClus
ClsVec findClus(SensorStripMap &sensorMap)
Method searching for clusters - first, strips above _SNseed threshold, so-called seed strips...
Definition: SiStripClus.cc:506
sistrip::SiStripClus::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Method called for each run - used for run header processing.
Definition: SiStripClus.cc:286
sistrip::StripCluster
This class holds all information about strip clusters, where the strip cluster is defined as a bunch ...
Definition: StripCluster.h:24
sistrip::SiStripClus::_TanOfAvgHLorentzShift
float _TanOfAvgHLorentzShift
Tangent of holes&#39; average Lorentz shift.
Definition: SiStripClus.h:142
sistrip::SiStripClus::_floatStripsZ
bool _floatStripsZ
Is every even strip floating in Z?
Definition: SiStripClus.h:146
sistrip::SiStripGeom::getSensorLength
virtual double getSensorLength(short int layerID) const
Get sensor length.
Definition: SiStripGeom.cc:467
sistrip::SiStripClus::_outColName
std::string _outColName
LCIO output collection name.
Definition: SiStripClus.h:131
sistrip::SiStripClus::_resSVDOtherInZ
std::vector< float > _resSVDOtherInZ
Mean strip resolution in Z - other layers; in [mm].
Definition: SiStripClus.h:153
ENDCOLOR
#define ENDCOLOR
Definition: Colours.h:21
sistrip::SiStripClus::_resSVDFirstInRPhi
std::vector< float > _resSVDFirstInRPhi
Mean strip resolution in RPhi - 1st layer; in [mm].
Definition: SiStripClus.h:150
sistrip::SiStripClus::printProcessorParams
void printProcessorParams() const
Method printing processor parameters.
Definition: SiStripClus.cc:2059
sistrip::SensorStripMap
std::map< int, StripChargeMap * > SensorStripMap
Definition: SiStripDigi.h:142
sistrip::SimTrackerHitMap
std::map< EVENT::SimTrackerHit *, float > SimTrackerHitMap
Definition: Signal.h:17
sistrip::StripCluster::getStripFront
int getStripFront() const
Get front strip ID.
Definition: StripCluster.h:150
sistrip::SiStripClus::init
virtual void init()
Method called at the beginning of data processing - used for initialization.
Definition: SiStripClus.cc:229
sistrip::SiStripClus::_relColNamePlsToSim
std::string _relColNamePlsToSim
LCIO input relation collection name - TrackerPulse &lt;-&gt; SimTrkHit.
Definition: SiStripClus.h:132
sistrip::SiStripClus::_SNadjacent
float _SNadjacent
Signal to noise ratio cut for adjacent strips.
Definition: SiStripClus.h:137
sistrip::ConstStringVec
const std::vector< std::string > ConstStringVec
Definition: SiStripClus.h:52
sistrip::SiStripGeom::getStripPosY
virtual double getStripPosY(const int &layerID, const int &sensorID, const int &stripID, const double &posZ) const =0
Get strip y-position.
sistrip::SiStripClus::storeHitsAdjacents
StripChargeMap & storeHitsAdjacents(StripChargeMap &clsStripsIn, It schmap, const It &endIt, StripChargeMap &cM)
Calculated and stored the hits adjacents.
Definition: SiStripClus.cc:1645
sistrip::Signal
Signal class holds all information that one gets from strip, pixel, etc ...
Definition: Signal.h:24
sistrip::fC
static const float fC
Definition: PhysicalConstants.h:32
sistrip::SiStripClus::_resSVDOtherInRPhi
std::vector< float > _resSVDOtherInRPhi
Mean strip resolution in RPhi - other layers; in [mm].
Definition: SiStripClus.h:151
DUNDERL
#define DUNDERL
Definition: Colours.h:20
sistrip::StripCluster::getSimHitMap
const SimTrackerHitMap & getSimHitMap() const
Get MC truth information about SimTrackerHits, which contributed.
Definition: StripCluster.h:144
sistrip::StripCluster::get3PosSigma
CLHEP::Hep3Vector get3PosSigma() const
Get cluster position Three vector.
Definition: StripCluster.h:132
sistrip::SiStripGeom::getLadderThick
virtual double getLadderThick(short int layerID) const
Get ladder thickness.
Definition: SiStripGeom.cc:256
sistrip::StripCluster::getPosY
double getPosY() const
Get cluster position Y.
Definition: StripCluster.h:114
sistrip::SiStripClus::_resSVDFirstInZ
std::vector< float > _resSVDFirstInZ
Mean strip resolution in Z - 1st layer; in [mm].
Definition: SiStripClus.h:152
sistrip::SiStripClus::processEvent
virtual void processEvent(LCEvent *event)
Method called for each event - used for event data processing.
Definition: SiStripClus.cc:301
sistrip::SiStripClus::_nRun
int _nRun
Run number.
Definition: SiStripClus.h:198
sistrip::StripCluster::getLadderID
short int getLadderID() const
Get cluster ladder ID.
Definition: StripCluster.h:105
sistrip::anSiStripClus
SiStripClus anSiStripClus
Definition: SiStripClus.cc:44
sistrip::SiStripClus::_geometry
SiStripGeom * _geometry
All geometry information from Gear xml file.
Definition: SiStripClus.h:147
sistrip::Signal::updateCharge
void updateCharge(double charge)
Update signal.
Definition: Signal.h:41
sistrip::SiStripGeom::getCrossLinePoint
virtual CLHEP::Hep3Vector getCrossLinePoint(const int &diskID, const int &sensorID, const int &stripIDFront, const int &stripIDRear) const =0
Get the point which crossed two strips.
sistrip::StripCluster::getCharge
double getCharge() const
Get cluster charge.
Definition: StripCluster.h:135
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
sistrip::SiStripClus::_CMSnoise
float _CMSnoise
Common mode subtracted noise, set in ENC.
Definition: SiStripClus.h:135
DGREEN
#define DGREEN
Definition: Colours.h:19
sistrip::StripCluster::getSize
short int getSize() const
Get cluster size.
Definition: StripCluster.h:141
sistrip::SiStripClus::_floatStripsRPhi
bool _floatStripsRPhi
Is every even strip floating in R-Phi?
Definition: SiStripClus.h:145
sistrip::SiStripClus::_nEvent
int _nEvent
Event number.
Definition: SiStripClus.h:199
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
sistrip::SiStripClus::_inColName
std::string _inColName
LCIO input collection name.
Definition: SiStripClus.h:130
sistrip::SiStripGeom::cellID0withSensor0
virtual int cellID0withSensor0(const int &cellID0) const =0
Returns the input cellID0 where the field sensor is put to 0.
sistrip::StripCluster::getStripRear
int getStripRear() const
Get rear strip ID.
Definition: StripCluster.h:153
sistrip::SiStripGeom::decodeStripID
virtual std::pair< StripType, int > decodeStripID(const int &encStripID) const =0
Decode stripID.
sistrip::StripCluster::getLayerID
short int getLayerID() const
Get cluster layer ID.
Definition: StripCluster.h:102
sistrip::ClsVec
std::vector< StripCluster * > ClsVec
Definition: SiStripClus.h:53
sistrip::StripChargeMap
std::map< int, Signal * > StripChargeMap
Definition: SiStripDigi.h:141
sistrip::SiStripGeom::getLadderTheta
virtual double getLadderTheta(short int layerID) const =0
Get ladder rotation - theta angle.
STRIPREAR
#define STRIPREAR
Definition: SiStripDigi.h:124
sistrip::SiStripGeom::getSensorPitch
virtual double getSensorPitch(const int &layerID, const int &sensorID, const double &posZ) const =0
Get sensor pitch.
sistrip::StripClusterPair
std::pair< int, StripCluster * > StripClusterPair
Definition: SiStripClus.cc:503
sistrip::SiStripGeom::transformPointToGlobal
virtual CLHEP::Hep3Vector transformPointToGlobal(short int layerID, short int ladderID, short int sensorID, const CLHEP::Hep3Vector &point)=0
Transform given point from local ref. system (sensor) to global ref. system.
sistrip::SiStripClus::calcHits
void calcHits(ClsVec &clsVec, IMPL::LCCollectionVec *colOfTrkHits)
Method calculating hits from given clusters.
Definition: SiStripClus.cc:1680
sistrip::SiStripClus::printHitInfo
void printHitInfo(const StripCluster *pCluster) const
Method printing hit info.
Definition: SiStripClus.cc:2093
sistrip::SiStripGeom::decodeCellID
virtual std::map< std::string, int > decodeCellID(const UTIL::BitField64 &cellDec) const =0
Encode cellID.
sistrip::SiStripGeom::getLayerRealID
virtual int getLayerRealID(short int layerID) const
Get layer real ID.
Definition: SiStripGeom.cc:123
sistrip::SiStripClus::_pitchFront
double _pitchFront
Pitch in the middle of the front sensor.
Definition: SiStripClus.h:159
sistrip::pi
static const double pi
Definition: PhysicalConstants.h:69
sistrip::SiStripClus::_navigatorPls
LCRelationNavigator * _navigatorPls
Definition: SiStripClus.h:156
sistrip::SiStripGeom::getSensorWidthMax
virtual double getSensorWidthMax(short int layerID) const
Get sensor width (the wider one for forward-type sensors)
Definition: SiStripGeom.cc:435
sistrip::SiStripClus::_SNseed
float _SNseed
Signal to noise ratio cut for seed strips.
Definition: SiStripClus.h:136
sistrip::SiStripGeom::transformVecToGlobal
virtual CLHEP::Hep3Vector transformVecToGlobal(short int layerID, short int ladderID, short int sensorID, const CLHEP::Hep3Vector &vec)=0
Transform given vector from local ref. system (sensor) to global ref. system.
sistrip::SensorStripClusterMap
std::map< int, std::map< StripType, std::vector< StripClusterPair > > > SensorStripClusterMap
Definition: SiStripClus.cc:504
sistrip::SiStripClus::releaseMap
void releaseMap(SensorStripMap &sensorMap)
Method to release memory of the SensorStripMap.
Definition: SiStripClus.cc:2026
sistrip::StripCluster::get3Position
CLHEP::Hep3Vector get3Position() const
Get cluster position Three vector.
Definition: StripCluster.h:120
sistrip::SiStripClus::_subdetector
std::string _subdetector
Name of the subdetector to be clusterize.
Definition: SiStripClus.h:162