Attachment 'KIT.cc'
Download 1 #include "KITdev.h"
2
3 using namespace lcio ;
4 using namespace marlin ;
5 using namespace std;
6 using namespace IMPL;
7 using namespace EVENT;
8 using namespace UTIL;
9
10
11 /*
12 Photon finding Algorythm by Predrag Krstonosic.
13 depends on Superhits, KITutil (for calibration),
14 Physics Geometry Database myPGDB (for pre-calculation needed for NN Clustering),
15 MarlinMath, B_Util
16
17 reads in CalorimeterHits, and writes out reconstructed Photons
18 developed and optimized for low energetic photons coming from the IP in LDC-00_02
19
20 Adopted to LDCPrime_02Sc by Nanda Wattimena
21 */
22
23 KITdev aKITdev ;
24 KITdev::KITdev() : Processor("KITdev") {
25
26 _description = "KITdev: an emShowerFinder for LDC" ;
27
28 registerProcessorParameter( "ECAL_Collection",
29 "Name of the CalorimeterHit Collection for ECAL ",
30 _Ecal_col,
31 std::string("ECAL"));
32 registerProcessorParameter( "HCAL_Collection",
33 "Name of the CalorimeterHit Collection for HCAL ",
34 _Hcal_col,
35 std::string("HCAL"));
36
37 registerProcessorParameter( "Core_Collection",
38 "Name of the Cluster Collection for Cores ",
39 _CoreCollection,
40 std::string("CORE"));
41
42 registerProcessorParameter( "detector_Model",
43 "chose detector geometry: LDC00, LDC0105, LDCPrime", //TODO: use this information!
44 _detectorModel,
45 std::string("LDCPrime"));
46
47 registerProcessorParameter( "Cleaning",
48 "Apply topological cleaning?",
49 _ToClean,
50 std::string("YES"));
51 registerProcessorParameter( "TopologicalCut",
52 "if so, require >= N+1 neighbors",
53 _CleanCut,
54 (int)5);
55
56 registerProcessorParameter( "NumberOfLevels",
57 "Number of levels for central loop ",
58 _Nlevels,
59 (int)10);
60 //defult values of energy per hit [mip], at which to start a new level
61 vector<float> mipstep;
62 mipstep.push_back(0.1);
63 mipstep.push_back(1.5);
64 mipstep.push_back(2.5);
65 mipstep.push_back(4.0);
66 mipstep.push_back(6.0);
67 mipstep.push_back(9.0);
68 mipstep.push_back(16.0);
69 mipstep.push_back(26.0);
70 mipstep.push_back(41.0);
71 mipstep.push_back(65.0);
72
73 registerProcessorParameter( "Levels",
74 "level thresholds (energy per hit [MIP]) ",
75 _mipstep,
76 mipstep);
77
78 registerProcessorParameter( "MinHit0",
79 "Minimal Number of hits for ground level cluster ",
80 _MinHit0,
81 (int)4);
82
83 registerProcessorParameter( "MinHitSplit",
84 "Minimal Number of hits for i-th level cluster ",
85 _MinHitSplit,
86 (int)2);
87
88 registerProcessorParameter( "rCut",
89 "Fluctuation suprresion cut",
90 _rCut,
91 (double)0.4);
92 registerProcessorParameter( "distCut",
93 "Square of distance cut [mm] for merging ",
94 _distCut,
95 (double)35.0);
96 registerProcessorParameter( "cosCut",
97 "Cosine of the angle for merging ",
98 _cosCut,
99 (double)0.95);
100
101 registerProcessorParameter( "probabilityDensityCut",
102 "cut on the probability density to assign hits to shower cores",
103 _probabilityDensityCut,
104 (double)1E-3);
105
106
107 registerProcessorParameter( "energyDeviationCut",
108 "cut on energy deviation of em shower candidates from estimated energy ( abs( (Ecluster - Eestimated)/Eestimated ) < cut )",
109 _energyDeviationCut,
110 (double)0.25);
111
112
113 /****** debugging tools: ******/
114 registerProcessorParameter( "DebugLevel",
115 "limits the amount of information written to std out (0 - none, 9 - maximal information)",
116 _debugLevel,
117 int(0) );
118 }
119
120
121
122
123 void KITdev::init() {
124
125 printParameters() ;
126 _nRun = 0 ;
127 _nEvt = 0 ;
128 _nSkippedEvt = 0 ;
129
130 // FIXME: hard coded cell id's for old Mokka (e.g. Mokka v5.4) versions)
131 CellIDDecoder<CalorimeterHit>::setDefaultEncoding("M:3,S-1:3,I:9,J:9,K-1:6");
132
133
134 _Ehit_ecal = 0.0;
135 _Ehit_hcal = 0.0;
136
137 _Et_miss = 0.0;
138 _Et_x = 0.0;
139 _Et_y = 0.0;
140
141 _NhitsEcal = 0;
142 _NhitsHcal = 0;
143 _NhitsCalos = 0;
144
145 _trueDist = 0;
146 _trueEgamma = 0;
147 _truePhi = 0;
148 _trueTheta = 0;
149
150 _trueEneutralino = 0;
151
152 }
153
154
155 void KITdev::processRunHeader( LCRunHeader* run) {
156
157 ++_nRun;
158 }
159
160 void KITdev::processEvent( LCEvent * evt ) {
161
162
163 try {
164 //open all kind of collections
165 LCCollection* ecalColl = evt->getCollection(_Ecal_col.c_str()) ;
166 LCCollection* hcalColl = evt->getCollection(_Hcal_col.c_str());
167
168 LCCollection* trueColl = evt->getCollection("TrueClusters");
169 LCCollection* mcColl = evt->getCollection("MCParticle");
170 LCCollection* LCRColl = evt->getCollection("TrueClusterToMCP");
171
172
173 std::vector<float> trueDistance;
174 std::vector<float> trueGammaEnergy;
175 std::vector<float> trueGammaStartX;
176 std::vector<float> trueGammaStartY;
177 std::vector<float> trueGammaStartZ;
178 std::vector<float> trueCogX;
179 std::vector<float> trueCogY;
180 std::vector<float> trueCogZ;
181 std::vector<float> trueGammaDirX;
182 std::vector<float> trueGammaDirY;
183 std::vector<float> trueGammaDirZ;
184 std::vector<float> trueGammaPosX;
185 std::vector<float> trueGammaPosY;
186 std::vector<float> trueGammaPosZ;
187 std::vector<float> trueNeutralinoEnergy;
188 std::vector<float> trueNeutralinoDirX;
189 std::vector<float> trueNeutralinoDirY;
190 std::vector<float> trueNeutralinoDirZ;
191 std::vector<float> trueThetaDir;
192 std::vector<float> truePhiDir;
193 std::vector<float> closestDistance;
194 std::vector<int> nClusterHits;
195 std::vector<float> phiPos;
196 std::vector<float> thetaPos;
197
198 /* get energy and vertex from photon that initiated the cluster (MCtruth) */
199 _Ntrcl = 0;
200 if ( trueColl!=0 ) {
201 _Ntrcl = trueColl->getNumberOfElements();
202
203
204 for (int i=0; i<_Ntrcl; i++) {
205
206 Cluster* cluster = dynamic_cast<Cluster*>(trueColl->getElementAt(i));
207 LCRelationNavigator* nav = new LCRelationNavigator(LCRColl);
208 const LCObjectVec& relMCParticlesToCluster = nav->getRelatedToObjects(cluster);
209
210 if ( relMCParticlesToCluster.size() > 1 )
211 std::cout << "Error: More than one MCParticle related to cluster." << std::endl;
212
213 // FIXME: here is only the 0th contribution taken into account
214 MCParticle* mcpOfCluster = dynamic_cast<MCParticle*>(relMCParticlesToCluster.at(0));
215 MCParticle* mcp = dynamic_cast<MCParticle*>(mcColl->getElementAt(0));
216
217 int mcpOfClusterID = mcpOfCluster->getPDG();
218 int parentID = mcp->getPDG();
219 int nDaughters = mcp->getDaughters().size();
220 int daughterID[nDaughters];
221
222 //get true gamma information
223 if(mcpOfClusterID == 22) {
224 _trueEgamma = mcpOfCluster->getEnergy();
225 _trueGammaPos[0] = mcpOfCluster->getVertex()[0];
226 _trueGammaPos[1] = mcpOfCluster->getVertex()[1];
227 _trueGammaPos[2] = mcpOfCluster->getVertex()[2];
228
229 _trueDist = sqrt( pow(_trueGammaPos[0],2) + pow(_trueGammaPos[1],2) + pow(_trueGammaPos[2],2) );
230
231 _trueGammaDir[0] = mcpOfCluster->getMomentum()[0];
232 _trueGammaDir[1] = mcpOfCluster->getMomentum()[1];
233 _trueGammaDir[2] = mcpOfCluster->getMomentum()[2];
234
235 double norm = sqrt( pow(_trueGammaDir[0],2) + pow(_trueGammaDir[1],2) + pow(_trueGammaDir[2],2) );
236
237 _trueGammaDir[0] /= norm;
238 _trueGammaDir[1] /= norm;
239 _trueGammaDir[2] /= norm;
240
241 _truePhi = atan2(_trueGammaDir[1],_trueGammaDir[0]) *180/M_PI;
242 _trueTheta = acos(_trueGammaDir[2]) *180/M_PI;
243 }
244 //get true neutralino information
245 if( parentID == 1000022) {
246 _trueEneutralino = mcp->getEnergy();
247
248 _trueNeutralinoDir[0] = mcp->getMomentum()[0];
249 _trueNeutralinoDir[1] = mcp->getMomentum()[1];
250 _trueNeutralinoDir[2] = mcp->getMomentum()[2];
251
252 double norm = sqrt( pow(_trueNeutralinoDir[0],2) + pow(_trueNeutralinoDir[1],2) + pow(_trueNeutralinoDir[2],2) );
253
254 _trueNeutralinoDir[0] /= norm;
255 _trueNeutralinoDir[1] /= norm;
256 _trueNeutralinoDir[2] /= norm;
257
258 if (_debugLevel >1 ) {
259 std::cout<<"/*** MC truth ***/"<<std::endl;
260 std::cout<<"found neutralino "<<std::endl;
261 std::cout<<" - with energy: "<<_trueEneutralino<<" GeV"<<std::endl;
262 std::cout<<" - with momentum: "<<_trueNeutralinoDir[0]<<", "<<_trueNeutralinoDir[1]<<", "<<_trueNeutralinoDir[2]<<std::endl;
263 std::cout<<"/*** *** *** ***/"<<std::endl;
264 }
265 }//is neutralino
266
267 else {
268 if (_debugLevel >1 )
269 std::cout<<"!!! Did not find neutralino !!!"<<std::endl;
270
271 _trueEneutralino = 0;
272 _trueNeutralinoDir[0] = 0;
273 _trueNeutralinoDir[1] = 0;
274 _trueNeutralinoDir[2] = 0;
275 }
276
277
278 if ( _debugLevel > 0 ) {
279 std::cout<<"/*** MC truth ***/"<<std::endl;
280 std::cout<<"found cluster "<<std::endl;
281 std::cout<<" - with energy: "<<_trueEgamma<<" GeV"<<std::endl;
282 std::cout<<" - at position: "<<_trueGammaPos[0]<<", "<<_trueGammaPos[1]<<", "<<_trueGammaPos[2]<<std::endl;
283 std::cout<<" - distance to IP: "<<_trueDist<<" mm"<<std::endl;
284 std::cout<<" - momentum: "<<_trueGammaDir[0]<<", "<<_trueGammaDir[1]<<", "<<_trueGammaDir[2]<<std::endl;
285 std::cout<<" - phi: "<<_truePhi<<" degree"<<std::endl;
286 std::cout<<" - theta: "<<_trueTheta<<" degree"<<std::endl;
287 std::cout<<"/*** *** *** ***/"<<std::endl;
288 }
289
290
291 trueDistance.push_back(_trueDist);
292 trueGammaEnergy.push_back(_trueEgamma);
293 trueGammaStartX.push_back(_trueGammaPos[0]);
294 trueGammaStartY.push_back(_trueGammaPos[1]);
295 trueGammaStartZ.push_back(_trueGammaPos[2]);
296 trueCogX.push_back(_trueCog[0]);
297 trueCogY.push_back(_trueCog[1]);
298 trueCogZ.push_back(_trueCog[2]);
299 trueGammaDirX.push_back(_trueGammaDir[0]);
300 trueGammaDirY.push_back(_trueGammaDir[1]);
301 trueGammaDirZ.push_back(_trueGammaDir[2]);
302 trueThetaDir.push_back(_trueTheta);
303 truePhiDir.push_back(_truePhi);
304
305 trueNeutralinoEnergy.push_back(_trueEneutralino);
306 trueNeutralinoDirX.push_back(_trueNeutralinoDir[0]);
307 trueNeutralinoDirY.push_back(_trueNeutralinoDir[1]);
308 trueNeutralinoDirZ.push_back(_trueNeutralinoDir[2]);
309 }//all true cluster
310
311 if ( _debugLevel > 0 )
312 std::cout<<"KITdev found "<<_Ntrcl<<" true cluster in event "<<_nEvt<<std::endl;
313 }//true cluster collection not empty
314
315
316 CellIDDecoder<CalorimeterHit>CDECAL = CellIDDecoder<CalorimeterHit>(ecalColl);
317 CellIDDecoder<CalorimeterHit>CDHCAL = CellIDDecoder<CalorimeterHit>(hcalColl);
318
319 LCCollectionVec * clusterColl = new LCCollectionVec(LCIO::CLUSTER);
320
321 //set everything to 0
322 _E_estimated = 0.0;
323
324 _Esum_ecal = 0.0;
325 _Esum_hcal = 0.0;
326 _Esum_calos = 0.0;
327
328 _Et_miss = 0.0;
329 _Et_x = 0.0;
330 _Et_y = 0.0;
331
332 _phiDegree = 0;
333 _thetaDegree = 0;
334 _dca = 0;
335 _nPhotons = 0;
336
337
338 unsigned int nHcalElements;
339 unsigned int nEcalElements;
340 if ( ecalColl!=0 ) {
341 if ( _debugLevel > 0 )
342 std::cout<<"KITdev found collection "<<_Ecal_col.c_str()<<std::endl;
343
344 nEcalElements = ecalColl->getNumberOfElements();
345
346 //get ECAL hits and energy in this event
347 _NhitsEcal = nEcalElements;
348
349
350 for (unsigned int i=0; i<nEcalElements; i++) {
351 CalorimeterHit* ecalHit = dynamic_cast<CalorimeterHit*>(ecalColl->getElementAt(i));
352
353 unsigned int layer = CDECAL(ecalHit)["K-1"];
354
355 _Ehit_ecal = ecalHit->getEnergy();
356 _Esum_ecal += ecalHit->getEnergy();
357 _Esum_calos += ecalHit->getEnergy();
358
359 //calculate transverse energy in ECAL
360 float normalisation = sqrt ( pow(ecalHit->getPosition()[0],2) + pow(ecalHit->getPosition()[1],2) + pow(ecalHit->getPosition()[2],2));
361 _Et_x += ecalHit->getEnergy() * ecalHit->getPosition()[0] / normalisation;
362 _Et_y += ecalHit->getEnergy() * ecalHit->getPosition()[1] / normalisation;
363
364 if ( layer <= myPGDB[myPGdb::ECAL1_BAR].max_lay ) {
365
366 _Ehit_ecal_mip = _Ehit_ecal / myPGDB[myPGdb::ECAL1_BAR].mip_whole; //mip_whole=mip_vis*e_coeff;
367
368 if ( _debugLevel > 8 )
369 std::cout<<myPGDB[myPGdb::ECAL1_BAR].mip_vis<<" * "<<myPGDB[myPGdb::ECAL1_BAR].e_coeff<<"= mip_whole SF1: "<<myPGDB[myPGdb::ECAL1_BAR].mip_whole<<std::endl;
370
371 }
372 else {
373 _Ehit_ecal_mip = _Ehit_ecal / myPGDB[myPGdb::ECAL2_BAR].mip_whole;
374
375 if ( _debugLevel > 8 )
376 std::cout<<myPGDB[myPGdb::ECAL2_BAR].mip_vis<<" * "<<myPGDB[myPGdb::ECAL2_BAR].e_coeff<<"= mip_whole SF1: "<<myPGDB[myPGdb::ECAL2_BAR].mip_whole<<std::endl;
377
378 }
379
380 }//all ECAL elem.
381 if ( _debugLevel > 0 )
382 std::cout<<" with "<<_NhitsEcal<<" hits and "<<_Esum_ecal<<"[GeV] energy"<<std::endl;
383
384 //get HCAL hits and energy in this event
385 if ( hcalColl!=0 ) {
386 if ( _debugLevel > 0 )
387 std::cout<<"KITdev found collection "<<_Hcal_col.c_str()<<std::endl;
388
389 nHcalElements = hcalColl->getNumberOfElements();
390 _NhitsHcal = nHcalElements;
391
392 for (unsigned int i=0; i<nHcalElements; i++) {
393 CalorimeterHit* hcalHit = dynamic_cast<CalorimeterHit*>(hcalColl->getElementAt(i));
394
395 _Esum_hcal += hcalHit->getEnergy();
396 _Ehit_hcal = hcalHit->getEnergy();
397 _Esum_calos += hcalHit->getEnergy();
398
399 //calculate transverse energy in HCAL
400 float normalisation = sqrt ( pow(hcalHit->getPosition()[0],2) + pow(hcalHit->getPosition()[1],2) + pow(hcalHit->getPosition()[2],2));
401 _Et_x += hcalHit->getEnergy() * hcalHit->getPosition()[0] / normalisation;
402 _Et_y += hcalHit->getEnergy() * hcalHit->getPosition()[1] / normalisation;
403 }
404 _NhitsCalos = nHcalElements + nEcalElements;
405
406 if ( _debugLevel > 0 )
407 std::cout<<" with "<<_NhitsHcal<<" hits and "<<_Esum_hcal<<"[GeV] energy"<<std::endl;
408 }//hcal col != 0
409
410 //since energy conservation should be valid, the total transverse energy should be 0 and thus Et = - Et_miss
411 _Et_miss = sqrt ( pow(_Et_x,2) + pow(_Et_y,2) );
412
413 //save parameters for later processors
414 evt->parameters().setValue("_Et_miss",(float)_Et_miss);
415 evt->parameters().setValue("_Esum_ECAL",(float)_Esum_ecal);
416 evt->parameters().setValue("_Nhits_ECAL",(int)nEcalElements);
417 evt->parameters().setValue("_Esum_HCAL",(float)_Esum_hcal);
418 evt->parameters().setValue("_Nhits_HCAL",(int)nHcalElements);
419 evt->parameters().setValue("_Esum_calos",(float)_Esum_calos);
420 evt->parameters().setValue("_Nhits_calos",(int)_NhitsCalos);
421 evt->parameters().setValue("_trueNcluster",(int)_Ntrcl);
422
423 vector<Superhit2*> calo[10]; // array with vectors of superhit pointers
424 // creating all superhits, the vector in calo[0] is filled with the superhits; applies mip calibration
425 CreateAllShits2(ecalColl,CDECAL,calo,_debugLevel);
426
427 // precalculation, shifting cells to equal distance for NN finder -> KITutil
428 TotalPrecalc2(calo,CDECAL,nEcalElements,_CleanCut,_debugLevel);
429
430
431 // setting the parameters of the alghorithm
432 vector <PROTSEED2> photonSeed;
433 CoreCut Ccut;
434 Ccut.rCut = _rCut;
435 Ccut.distCut = _distCut;
436 Ccut.cosCut = _cosCut;
437 Ccut.MinHit0 = (unsigned int) _MinHit0;
438 Ccut.MinHitSplit = (unsigned int) _MinHitSplit;
439
440 //temporary storage of cluster candidates
441 const int myNLevels = _Nlevels;
442 vector<Tmpcl2*> clusterVector[myNLevels];
443
444 // apply topological cleaning
445 // -> remove hits with < _CleanCut neighbouring hits
446 if ( _ToClean=="YES" || _ToClean=="yes" || _ToClean=="Yes") {
447 if ( _debugLevel > 0 )
448 std::cout << " topological cleaning " << std::endl;
449
450 //main loop: level wise clustering -> KITutil
451 FindCores2( &(calo[4]), clusterVector, &photonSeed, _Nlevels, _mipstep, Ccut, _debugLevel );
452 } else {
453 FindCores2( &(calo[0]), clusterVector, &photonSeed, _Nlevels, _mipstep, Ccut, _debugLevel );
454 }
455
456 // container to store information to assign hits to photon candidates
457 std::vector<ECALHitWithAttributes> ECALHitsWithAttributes;
458
459 //calibration to get from measured cluster energy to estimated photon energy
460 //is hardcoded for different levels and energies in KITutil.cc
461 std::vector<CoreCalib2> coreCalibration;
462 CreateCalibration_LDCPrime02Sc(&coreCalibration);
463 //CreateCalibration_LDC00(&coreCalibration);
464 int nClusters = 0;
465 for ( unsigned int seedCounter=0; seedCounter<photonSeed.size(); seedCounter++) {
466 if ( photonSeed[seedCounter].active==true ) {
467 if ( _debugLevel > 1 )
468 std::cout<<"KIT found photon seed"<<std::endl;
469
470 //find cluster center
471 double centerPosition[3];
472 double directionOfCluster[3];
473 centerPosition[0] = photonSeed[seedCounter].cl->getCenter()[0];
474 centerPosition[1] = photonSeed[seedCounter].cl->getCenter()[1];
475 centerPosition[2] = photonSeed[seedCounter].cl->getCenter()[2];
476
477 // already normalised cluster direction
478 photonSeed[seedCounter].cl->findInertia();
479 directionOfCluster[0] = photonSeed[seedCounter].cl->direction[0];
480 directionOfCluster[1] = photonSeed[seedCounter].cl->direction[1];
481 directionOfCluster[2] = photonSeed[seedCounter].cl->direction[2];
482
483 double coreEnergy = 0.0;
484 double startPosition[3];
485 double startPositionMinProjection[3];
486 double startPositionMaxProjection[3];
487 double minProjectionAlongMainPrincipalAxis = DBL_MAX;
488 double maxProjectionAlongMainPrincipalAxis = 0.0;
489
490 //sum over hit in this cluster
491 for ( unsigned int hitCounter=0; hitCounter<photonSeed[seedCounter].cl->hits.size(); hitCounter++ ) {
492
493 //sum up core energy
494 coreEnergy += photonSeed[seedCounter].cl->hits[hitCounter]->chit->getEnergy();
495
496 //get cluster geometry
497 double vCenterToHit[3];
498 double projectionAlongMainPrincipalAxis = 0.0;
499
500 for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
501
502 vCenterToHit[positionCounter] = centerPosition[positionCounter] //vector from center to hit
503 - photonSeed[seedCounter].cl->hits[hitCounter]->chit->getPosition()[positionCounter];
504
505 projectionAlongMainPrincipalAxis += directionOfCluster[positionCounter]*vCenterToHit[positionCounter];
506 }//all positions
507
508 if ( projectionAlongMainPrincipalAxis < minProjectionAlongMainPrincipalAxis ) {
509
510 minProjectionAlongMainPrincipalAxis = projectionAlongMainPrincipalAxis;
511
512 for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
513 startPositionMinProjection[positionCounter] = centerPosition[positionCounter]
514 - 2.0*minProjectionAlongMainPrincipalAxis*directionOfCluster[positionCounter];
515 }
516 }
517
518 if ( projectionAlongMainPrincipalAxis > maxProjectionAlongMainPrincipalAxis ) {
519
520 maxProjectionAlongMainPrincipalAxis = projectionAlongMainPrincipalAxis;
521 for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
522 startPositionMaxProjection[positionCounter] = centerPosition[positionCounter]
523 - 2.0*maxProjectionAlongMainPrincipalAxis*directionOfCluster[positionCounter];
524 }
525 }
526
527 //get cluster start
528 double rStartPositionMinProjection = sqrt( pow(startPositionMinProjection[0],2)
529 + pow(startPositionMinProjection[1],2)
530 + pow(startPositionMinProjection[2],2));
531 double rStartPositionMaxProjection = sqrt( pow(startPositionMaxProjection[0],2)
532 + pow(startPositionMaxProjection[1],2)
533 + pow(startPositionMaxProjection[2],2));
534
535 if ( rStartPositionMinProjection < rStartPositionMaxProjection ) {
536
537 startPosition[0] = startPositionMinProjection[0];
538 startPosition[1] = startPositionMinProjection[1];
539 startPosition[2] = startPositionMinProjection[2];
540 }
541 else {
542
543 startPosition[0] = startPositionMaxProjection[0];
544 startPosition[1] = startPositionMaxProjection[1];
545 startPosition[2] = startPositionMaxProjection[2];
546 }
547 }//all hits in cluster
548
549 _xStart_cluster = startPosition[0];
550 _yStart_cluster = startPosition[1];
551 _zStart_cluster = startPosition[2];
552
553 //estimate photon energy depending on the cluster energy, level and detector calibration
554
555 double photonEnergy = giveMeEEstimate2(photonSeed[seedCounter].level,photonSeed[seedCounter].cl->getEnergy(),coreCalibration);
556
557 if ( photonEnergy > 0.0 ) {
558
559 Photon2* photonFinder = new Photon2(photonEnergy,centerPosition,startPosition);
560
561 // loop over all ECAL entries, to collect remaining hits
562 for ( unsigned int ecalHitCounter = 0; ecalHitCounter < nEcalElements; ++ecalHitCounter) {
563
564 double probabilityAndDistance[2];
565 CalorimeterHit* ECALHit = dynamic_cast<CalorimeterHit*>(ecalColl->getElementAt(ecalHitCounter));
566
567 photonFinder->Prob(ECALHit,_probabilityDensityCut,probabilityAndDistance);
568
569 bool isECALHitAlreadyAssigned = false;
570 int indexOfECALHitAlreadyAssigned = 0;
571
572 for ( unsigned int k = 0; k < ECALHitsWithAttributes.size(); ++k) {
573 if ( ECALHitsWithAttributes.at(k).ECALHit == ECALHit) {
574
575 isECALHitAlreadyAssigned = true;
576 indexOfECALHitAlreadyAssigned = k;
577 break;
578 }
579 }
580 //avoid double counting
581 if ( !isECALHitAlreadyAssigned ) {
582 if ( probabilityAndDistance[0] > 0.00001 ) {
583
584 ECALHitWithAttributes hitWithAttributes;
585
586 hitWithAttributes.ECALHit = ECALHit;
587 hitWithAttributes.relatedCores.push_back(&(photonSeed[seedCounter]));
588 hitWithAttributes.probabilitiesForThisECALHit.push_back(probabilityAndDistance[0]);
589 hitWithAttributes.distancesToCoresForThisECALHit.push_back(probabilityAndDistance[1]);
590 hitWithAttributes.estimatedEnergyPerCore.push_back(photonEnergy);
591
592 ECALHitsWithAttributes.push_back(hitWithAttributes);
593 }
594 }
595 else if ( probabilityAndDistance[0] > 0.00001 ) {
596
597 ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).relatedCores.push_back(&(photonSeed[seedCounter]));
598 ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).probabilitiesForThisECALHit.push_back(probabilityAndDistance[0]);
599 ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).distancesToCoresForThisECALHit.push_back(probabilityAndDistance[1]);
600 ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).estimatedEnergyPerCore.push_back(photonEnergy);
601
602 if ( _debugLevel > 4 )
603 std::cout<<" All hits assigned "<<std::endl;
604
605 }
606 }//ecalHitCounter
607
608 estimatedE.push_back(photonEnergy);
609 _nPhotons++;
610
611 if ( _debugLevel > 0 )
612 std::cout<<" found photon with energy "<<photonEnergy<<std::endl;
613 delete photonFinder;
614 photonFinder = 0;
615 }//photon energy > 0
616 else
617 std::cout<<" energy == 0 -> no photon"<<std::endl;
618 }//photonSeed is active
619 }//all cluster
620
621 //******* Collect the remaining hits ********//
622
623 //map hit to clusters
624 std::map<PROTSEED2*,ClusterImpl*> mapCoreToCluster;
625 std::map<PROTSEED2*,double> mapCoreToEstimatedEnergy;
626
627 for ( unsigned int selectedHitsCounter = 0; selectedHitsCounter < ECALHitsWithAttributes.size(); ++selectedHitsCounter ) {
628
629 PROTSEED2* coreToAddHitTo;
630 double estimatedEnergyOfCoreToAddHitTo = 0.0;
631
632
633 if ( ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.size() == 1 ) {
634
635 coreToAddHitTo = ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.at(0);
636 estimatedEnergyOfCoreToAddHitTo = ECALHitsWithAttributes.at(selectedHitsCounter).estimatedEnergyPerCore.at(0);
637 } else {
638
639 unsigned int indexOfMaxProbability = 0;
640 double maxProbability = 0.0;
641
642 for ( unsigned int selectedCoresCounter = 0; selectedCoresCounter < ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.size(); ++selectedCoresCounter) {
643 if ( ECALHitsWithAttributes.at(selectedHitsCounter).probabilitiesForThisECALHit.at(selectedCoresCounter) > maxProbability ) {
644
645 indexOfMaxProbability = selectedCoresCounter;
646 maxProbability =
647 ECALHitsWithAttributes.at(selectedHitsCounter).probabilitiesForThisECALHit.at(selectedCoresCounter);
648 }
649 }//selected cores
650
651 coreToAddHitTo = ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.at(indexOfMaxProbability);
652 estimatedEnergyOfCoreToAddHitTo =
653 ECALHitsWithAttributes.at(selectedHitsCounter).estimatedEnergyPerCore.at(indexOfMaxProbability);
654 }
655 std::map<PROTSEED2*,ClusterImpl*>::iterator position = mapCoreToCluster.find(coreToAddHitTo);
656
657 if ( position == mapCoreToCluster.end() ) {
658
659 mapCoreToCluster[coreToAddHitTo] = new ClusterImpl();
660 mapCoreToEstimatedEnergy[coreToAddHitTo] = estimatedEnergyOfCoreToAddHitTo;
661 position = mapCoreToCluster.find(coreToAddHitTo);
662
663 if ( _debugLevel > 4 )
664 std::cout<<" Mapped Cores to Clusters "<<std::endl;
665 }
666
667 (*position).second->addHit(ECALHitsWithAttributes.at(selectedHitsCounter).ECALHit,(float)1.0);
668 }//selected ecal hits
669
670
671 if ( _debugLevel > 0 )
672 std::cout<<" Core to Cluster Map size "<<mapCoreToCluster.size()<<std::endl;
673
674 int iteratorCounter = 0;
675
676 for ( std::map<PROTSEED2*,ClusterImpl*>::iterator clusterIterator = mapCoreToCluster.begin(); clusterIterator != mapCoreToCluster.end(); clusterIterator++ ) {
677
678 ClusterImpl* cluster = (*clusterIterator).second;
679 double estimatedEnergyOfCluster = mapCoreToEstimatedEnergy[(*clusterIterator).first];
680
681 // set energy and position of cluster
682 unsigned int caloHits = cluster->getCalorimeterHits().size();
683
684 float* xPos = new float[caloHits];
685 float* yPos = new float[caloHits];
686 float* zPos = new float[caloHits];
687 float* hitEnergy = new float[caloHits];
688 float energy = 0.0;
689
690 for ( unsigned int caloHitsCounter = 0; caloHitsCounter < caloHits; ++caloHitsCounter ) {
691
692 xPos[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getPosition()[0];
693 yPos[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getPosition()[1];
694 zPos[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getPosition()[2];
695 hitEnergy[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getEnergy();
696
697 energy += cluster->getCalorimeterHits().at(caloHitsCounter)->getEnergy();
698 }
699
700 // cut of EM candidtate clusters which differ more than _energyDeviationCut from enery estimate
701 double energyDeviation = fabs( (energy-estimatedEnergyOfCluster)/estimatedEnergyOfCluster );
702 //FIXME: This accounts only for the leading cluster!
703
704 // simple cut on the starting layer of the shower
705 // FIXME: this is not working properly for the edges of the ECAL and HCAL => need pseudo layer !!!
706 const unsigned int cutOnLayer = 10; // FIXME: remove hard-coded cut
707 const unsigned int contOnNumberOfHitsInLayerCut = 4; // FIXME: remove hard-coded cut
708 unsigned int nHitsInLayerCut = 0;
709 bool hitWithinCutRangeFound = false;
710
711 // double energyDeviation = 0;
712
713
714 _E_estimated = estimatedE.at(iteratorCounter);
715 // energyDeviation = fabs( (energy - _E_estimated) / _E_estimated );
716
717 if (_debugLevel > 0) {
718 std::cout<<"energyDeviation: (" <<energy<<" - "<<_E_estimated<<") / "<<_E_estimated
719 <<" = "<<energyDeviation<<std::endl;
720 }
721
722 for ( unsigned int caloHitsCounter = 0; caloHitsCounter < cluster->getCalorimeterHits().size(); ++caloHitsCounter ) {
723
724 unsigned int layer = CDECAL(cluster->getCalorimeterHits().at(caloHitsCounter))["K-1"];
725 if ( layer <= cutOnLayer ) {
726 ++nHitsInLayerCut;
727 }
728
729 if ( nHitsInLayerCut >= contOnNumberOfHitsInLayerCut ) {//if >= 4 layers hit within the first 10 layers
730 hitWithinCutRangeFound = true;
731 break;
732 }
733
734 //COMMENT IN && HITWITHINRANGE
735 if ( energyDeviation < _energyDeviationCut ) {
736 // if ( hitWithinCutRangeFound ) {
737 // flag hits
738 for ( unsigned int j = 0; j < cluster->getCalorimeterHits().size(); ++j ) {
739 cluster->getCalorimeterHits().at(j)->ext<isPartOfEMShowerCandidate>() = 1;
740 //clusterColl->addHit(cluster->getCalorimeterHits().at(j));
741 }
742 ClusterShapes* clusterShape = new ClusterShapes(caloHits,hitEnergy,xPos,yPos,zPos);
743
744 float position[3];
745 position[0] = clusterShape->getCentreOfGravity()[0];
746 position[1] = clusterShape->getCentreOfGravity()[1];
747 position[2] = clusterShape->getCentreOfGravity()[2];
748
749 float normalisation = sqrt( pow(position[0],2) + pow(position[1],2) + pow(position[2],2) );
750 float normPos[3];
751 normPos[0] = position[0]/normalisation;
752 normPos[1] = position[1]/normalisation;
753 normPos[2] = position[2]/normalisation;
754
755 float _phiPos = atan2(normPos[1],normPos[0]);
756 float _thetaPos = acos(normPos[2]);
757
758 //get main cluster axis
759 float mainClusterAxis[3];
760 mainClusterAxis[0] = clusterShape->getEigenVecInertia()[0];
761 mainClusterAxis[1] = clusterShape->getEigenVecInertia()[1];
762 mainClusterAxis[2] = clusterShape->getEigenVecInertia()[2];
763 //get cluster angle
764 float phiCluster = atan2(clusterShape->getEigenVecInertia()[1],clusterShape->getEigenVecInertia()[0]);
765 float thetaCluster = acos(clusterShape->getEigenVecInertia()[2]);
766
767 std::vector<float> shapeParameter(3) ;
768
769 shapeParameter[0] = mainClusterAxis[0];
770 shapeParameter[1] = mainClusterAxis[1];
771 shapeParameter[2] = mainClusterAxis[2];
772
773 cluster->setPosition(position);
774 cluster->setEnergy(energy);
775 cluster->setITheta(thetaCluster);
776 cluster->setIPhi(phiCluster);
777 cluster->setShape(shapeParameter);
778 clusterColl->addElement(cluster);
779
780 nClusters++;
781
782 nClusterHits.push_back(cluster->getCalorimeterHits().size());
783 phiPos.push_back(_phiPos);
784 thetaPos.push_back(_thetaPos);
785
786 if(energy == 0.)
787 std::cout<<"WARNING: found cluster with 0 GeV energy!"<<std::endl;
788
789
790 _phiDegree = phiCluster*180/M_PI;
791 _thetaDegree = thetaCluster*180/M_PI;
792
793 //glue cluster direction to cog
794 float photonLine[3];
795 for (int axisCounter=0; axisCounter<3; axisCounter++) {
796 photonLine[axisCounter] = position[axisCounter] + mainClusterAxis[axisCounter]; //g=cog+k*EV
797 }
798
799 //calculate distance of closest approach between gamma and IP:
800 // d = | cog x EV_1 |
801 float dcaVec[3];
802 dcaVec[0] = position[1]*mainClusterAxis[2] - position[2]*mainClusterAxis[1];
803 dcaVec[1] = position[2]*mainClusterAxis[0] - position[0]*mainClusterAxis[2];
804 dcaVec[2] = position[0]*mainClusterAxis[1] - position[1]*mainClusterAxis[0];
805
806 _dca = sqrt( pow(dcaVec[0],2) + pow(dcaVec[1],2)+ pow(dcaVec[2],2) );
807
808 closestDistance.push_back(_dca);
809
810 if ( _debugLevel > 0 ) {
811 std::cout<<" ************************* "<<std::endl;
812 std::cout<<"Final cluster with energy: "<<energy
813 <<" at position "<<position[0]
814 <<" , "<<position[1]
815 <<" , "<<position[2]<<std::endl;
816 std::cout<<"Cluster direction phi: "<<_phiDegree<<std::endl;
817 std::cout<<" theta: "<<_thetaDegree<<std::endl;
818 std::cout<<" main axis: "<<mainClusterAxis[0]
819 <<" , "<<mainClusterAxis[1]
820 <<" , "<<mainClusterAxis[2]<<std::endl;
821 std::cout<<"minimal distance from IP: "<<_dca<<std::endl;
822 std::cout<<" ************************* "<<std::endl;
823 }//DEBUG
824
825
826 delete clusterShape;
827 clusterShape = 0;
828 break;
829 // }//hit within range cut
830 // else std::cout << "hit not within range cut " << std::endl;
831
832 }//energy devitation cut
833 else if (_debugLevel >0) std::cout << "energy deviation "<<energyDeviation<<" > " << _energyDeviationCut << std::endl;
834
835
836 }//all calorimeter hits
837
838 delete[] xPos;
839 xPos = 0;
840 delete[] yPos;
841 yPos = 0;
842 delete[] zPos;
843 zPos = 0;
844 delete[] hitEnergy;
845 hitEnergy = 0;
846
847 iteratorCounter++;
848 }//in cluster map
849
850 //save parameters for rootTreeWriter
851 clusterColl->parameters().setValues("_dca",closestDistance);
852 clusterColl->parameters().setValues("_nClusterHits",nClusterHits);
853 clusterColl->parameters().setValues("_phiPos",phiPos);
854 clusterColl->parameters().setValues("_thetaPos",thetaPos);
855 clusterColl->parameters().setValues("_trueDist",trueDistance);
856 clusterColl->parameters().setValues("_trueGammaEnergy",trueGammaEnergy);
857 clusterColl->parameters().setValues("_trueGammaPosX",trueGammaPosX);
858 clusterColl->parameters().setValues("_trueGammaPosY",trueGammaPosY);
859 clusterColl->parameters().setValues("_trueGammaPosZ",trueGammaPosZ);
860 clusterColl->parameters().setValues("_trueGammaDirX",trueGammaDirX);
861 clusterColl->parameters().setValues("_trueGammaDirY",trueGammaDirY);
862 clusterColl->parameters().setValues("_trueGammaDirZ",trueGammaDirZ);
863 clusterColl->parameters().setValues("_trueNeutralinoEnergy",trueNeutralinoEnergy);
864 clusterColl->parameters().setValues("_trueNeutralinoDirX",trueNeutralinoDirX);
865 clusterColl->parameters().setValues("_trueNeutralinoDirY",trueNeutralinoDirY);
866 clusterColl->parameters().setValues("_trueNeutralinoDirZ",trueNeutralinoDirZ);
867 clusterColl->parameters().setValues("_trueThetaDir",trueThetaDir);
868 clusterColl->parameters().setValues("_truePhiDir",truePhiDir);
869
870 mapCoreToCluster.clear();
871 ECALHitsWithAttributes.clear();
872 estimatedE.clear();
873
874 // for strong memory and nice dreams clean up again
875 for ( int levelCounter=0; levelCounter<_Nlevels; levelCounter++ ) {
876 if ( clusterVector[levelCounter].size()!=0 )
877 for ( unsigned int calibrationCounter=0; calibrationCounter<clusterVector[levelCounter].size(); calibrationCounter++ ) {
878 delete clusterVector[levelCounter][calibrationCounter];
879 }
880 }
881
882 for ( unsigned int caloCounter=0; caloCounter<2; caloCounter++ ) {
883 if ( calo[caloCounter].size()!=0 )
884 for ( unsigned int caloSize=0; caloSize<calo[caloCounter].size(); caloSize++ )
885 delete (calo[caloCounter][caloSize]);
886 }
887 }//ecalCollection not empty
888
889 evt->addCollection(clusterColl,_CoreCollection.c_str());
890 evt->parameters().setValue("_nGamma",(int)_nPhotons);
891
892 } catch ( DataNotAvailableException &e ) {
893 std::cout << "KITdev:: missing collection in event " << _nEvt << std::endl;
894 }
895
896 _nEvt++ ;
897
898 }//process event
899
900 void KITdev::end(){
901
902 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.