Attachment 'KITutil.cc'
Download 1 #include "KITutil_LDC00.h"
2
3
4 // Will be called by KIT.cc
5 // requires inout from myPGDB & MarlinMath
6
7 //Shits = superHits
8 void CreateAllShits2(LCCollection* colt,CellIDDecoder<CalorimeterHit>& id,vector<Superhit2*>* calo, int debug)
9 {
10 unsigned int n_elements = colt->getNumberOfElements();
11 int superHitCounterSampling1 = 0;
12 int superHitCounterSampling2 = 0;
13
14 for ( unsigned int elementCounter=0; elementCounter<n_elements; elementCounter++) {
15 CalorimeterHit* caloHit = dynamic_cast<CalorimeterHit*>(colt->getElementAt(elementCounter));
16 unsigned int layer = (unsigned int)id(caloHit)["K-1"];
17 Superhit2 *superHit;
18
19 //requires information about detector geometry from Physics Geometry database myPGDB
20 //for 2 different ECAL samplings in barrel ECAL1 & ECAL2
21 if ( layer<=myPGDB[myPGdb::ECAL1_BAR].max_lay ) {//maxlayer from gear-file
22 superHit = new Superhit2(myPGDB[myPGdb::ECAL1_BAR].mip_whole,caloHit); //mip_whole is hardcoded via e_coeff in myPGDB == energy of a mip?
23 superHitCounterSampling1++;
24 } else {
25 superHit = new Superhit2(myPGDB[myPGdb::ECAL2_BAR].mip_whole,caloHit);
26 superHitCounterSampling2++;
27 }
28 calo[0].push_back(superHit);
29 }
30 if (debug > 3 ) {
31 std::cout<<"DEBUG KITutil: * from "<< n_elements <<" calorimeter hits created "
32 << superHitCounterSampling1 <<" superhits in sampling 1 and "
33 << superHitCounterSampling2 <<" superhits in sampling 2"<<std::endl;
34 }
35 }
36
37 //precalculations does nearest neighbour clustering
38 //by shifting cells to equal distances in x, y & z
39 //called by KIT
40 void TotalPrecalc2(vector<Superhit2*>* calo,CellIDDecoder<CalorimeterHit>& id,unsigned int nelem,int Ncut,int debug)
41 {
42 if ( debug > 3 )
43 std::cout<<"DEBUG KITutil: - precalculation for "<<nelem<<" ECAL hits "<<std::endl;
44
45 vector<Superhit2*> tmpsh[9]; //fields 0 to 7 are superhits in stave 0 to 7, field 8 holds the endcaps
46
47 int barrelCounter = 0;
48 int endcapCounter = 0;
49
50 for ( unsigned int hitCounter=0; hitCounter<nelem; hitCounter++ ) {
51
52 CalorimeterHit* chh=(calo[0])[hitCounter]->chit;
53 int stave = id(chh)["S-1"]; //which octant
54 int module = id(chh)["M"]; //which module
55 int layer = id(chh)["K-1"]; //which layer
56
57 // some remaining superhit members initialised
58 (calo[0])[hitCounter]->S=stave;
59 (calo[0])[hitCounter]->M=module;
60 (calo[0])[hitCounter]->K=layer;
61 (calo[0])[hitCounter]->isCalorimeter=1; //is ecal
62
63 if ( module !=0 && module!=6 ) { //if hit not in endcaps
64 if ( layer!=0 ) { //and not in first layer
65 tmpsh[stave].push_back((calo[0])[hitCounter]);
66 } else {
67
68 //for hits in first layer, phi = angle between hit and y axis
69 //the following section determines if the first layer hit is a neighbour of the next stave.
70 //If so it is assigned to both stave vectors
71
72 double phi = atan2(-(calo[0])[hitCounter]->chit->getPosition()[0],(calo[0])[hitCounter]->chit->getPosition()[1]);
73 if ( phi< 0.0 ) {
74 if ( stave!=0 ) phi = 2.0*M_PI + phi; //for staves 1 - 7; M_PI = pi = 3.14
75 else phi = fabs(phi);
76 }
77
78 double diff = phi - M_PI_4*stave; //M_PI_4 = pi/4 = 0.79
79 // 45 degrees time is the angle between to corners of the octagon
80 // now the hit lies in stave 0 and if fabs(phi)<45 DEG, it is not overlapping, see below
81
82 //The hit is written to its original stave anyway
83 tmpsh[stave].push_back((calo[0])[hitCounter]);
84
85 if ( fabs(diff)<M_PI_4/2.0 ) { // hit in layer 0, but not in the interesting "overlap" area -> stored a 2nd time ?
86 tmpsh[stave].push_back((calo[0])[hitCounter]);
87 } else { // If its in the overlap area, put it in both stove collections
88
89 (calo[0])[hitCounter]->connect = true;
90
91 tmpsh[stave].push_back((calo[0])[hitCounter]);
92
93 int stavem1 = stave - 1;
94 // store hit in additional stave
95 if ( stavem1>=0 )
96 tmpsh[stavem1].push_back((calo[0])[hitCounter]);
97 else
98 tmpsh[7].push_back((calo[0])[hitCounter]);
99 }
100 } // special treatment for first layer (= layer 0)
101 barrelCounter++;
102
103 } // hit in barrel
104 else { //hits in endcap
105 tmpsh[8].push_back((calo[0])[hitCounter]);
106 endcapCounter++;
107 }
108 }//all elements
109 if ( debug > 3 )
110 std::cout<<"DEBUG KITutil: * found "<<endcapCounter<<" hits in ECAL endcap "<<" and "<<barrelCounter<<" hits in ECAL barrel"<<std::endl;
111
112 // Next section: find nearest neighbours
113 // CAREFULL: detector geometry dependend !!!
114
115 double celx = myPGDB[myPGdb::ECAL1_BAR].cell_size; // cell_size has to be checked
116 double distcut = 0.98*(2.0*celx)*(2.0*celx); //where does this come from?
117 // hits in barrel
118 for ( unsigned int staveCounter=0; staveCounter<8; staveCounter++ ) {
119 if(tmpsh[staveCounter].size()!=0)
120 Precalc2(tmpsh[staveCounter],myPGDB[myPGdb::ECAL1_BAR].r_inner,
121 myPGDB[myPGdb::ECAL1_CAP].z_inner,
122 myPGDB[myPGdb::ECAL1_BAR].cell_size,distcut,true,staveCounter,1,id,debug);
123 } // isCalorimeter = 1 for ecal, 2 for hcal, 0 else
124
125 if ( tmpsh[8].size()!=0 ) // hits in end caps: tmpsh[8]
126 Precalc2(tmpsh[8],myPGDB[myPGdb::ECAL1_BAR].r_inner,
127 myPGDB[myPGdb::ECAL1_CAP].z_inner,
128 myPGDB[myPGdb::ECAL1_BAR].cell_size,distcut,false,1,1,id,debug);
129
130 // rearrange calo[0]
131 for ( unsigned int ii=0; ii<calo[0].size(); ii++)
132 calo[0][ii]->top = calo[0][ii]->neighbours.size();
133
134 for ( unsigned int ii=0; ii<(calo[0]).size(); ii++) {
135 if(debug > 3)
136 std::cout<<"DEBUG KITutil: - found hit with "<<(calo[0])[ii]->top<<" neighbours"<<std::endl;
137
138 if ( (calo[0])[ii]->top !=0 ) {//if hit has neighbours
139 if ( (calo[0])[ii]->top >= 1 && (calo[0])[ii]->top < Ncut ) {
140 (calo[2]).push_back((calo[0])[ii]);// if hit with less than required neighbours -> no cluster
141 } else { //for topological cleaning
142 (calo[4]).push_back((calo[0])[ii]); // hit with required neighbours -> topological cleaning possible
143 }
144 } else { // if calo[0] hit has no neighbours at all
145 (calo[6]).push_back((calo[0])[ii]); //-> no cluster
146 }
147 }
148 }
149
150
151
152 void Precalc2 ( vector< Superhit2* >& shvec,double r, double z, double cell, double dist,bool isCalorimeter,int stave,int module,CellIDDecoder<CalorimeterHit>& id,int debug ) {
153
154 unsigned int nPoints = shvec.size();
155 ANNpointArray dataPoints;
156 ANNpoint queryPoint; //Annpoint: double[3]
157 ANNidxArray nnIndex;
158 ANNdistArray distance;
159 ANNkd_tree* kdTree;
160
161 queryPoint = annAllocPt(3); // 3d point
162 dataPoints = annAllocPts(nPoints,3);
163 nnIndex = new ANNidx[36];
164 distance = new ANNdist[36];
165 float xxx[3]; // transformed position in GridTransform
166 float Ecalradius = r;
167 float Ecalhalfz = z;
168 float Ecalcellsize = cell;
169
170 if ( debug>3 )
171 std::cout<<"DEBUG KITutil: shifting hits to new positions"<<std::endl;
172
173 for ( unsigned int hitCounter=0; hitCounter<nPoints; hitCounter++ ) {
174
175 GridTransform2(shvec[hitCounter]->chit,Ecalradius,Ecalhalfz,Ecalcellsize,xxx,isCalorimeter,stave,module,id);
176 //dataPoints are filled with transformed coordinates for search alg.
177 dataPoints[hitCounter][0] = xxx[0];
178 dataPoints[hitCounter][1] = xxx[1];
179 dataPoints[hitCounter][2] = xxx[2];
180 //point[] is the transformed position in class SuperHit
181 shvec[hitCounter]->point[0] = xxx[0];
182 shvec[hitCounter]->point[1] = xxx[1];
183 shvec[hitCounter]->point[2] = xxx[2];
184 }
185
186 //search for nearest neighbours with ANN code, assumed to work correctly see ANN.cpp
187 kdTree = new ANNkd_tree(dataPoints,nPoints,3);
188 ANNdist trd = dist;
189
190 for ( unsigned int hitCounter=0; hitCounter<nPoints; hitCounter++ ) {
191 //number of neighbours of each hit
192 unsigned int N_NearestNeighbours = (unsigned int)kdTree->annkFRSearch(dataPoints[hitCounter],trd,26,nnIndex,distance,0.0);
193 if ( N_NearestNeighbours>26 )
194 N_NearestNeighbours = 26;
195 if( isCalorimeter ) { // inside calorimeter: isCalorimeter(ecal)=1, isCalorimeter(hcal)=2
196 if( shvec[hitCounter]->connect ) { // hit in barrel and connected to two staves
197 for ( unsigned int neighbourCounter=0; neighbourCounter<N_NearestNeighbours; neighbourCounter++ ) {
198 if ( shvec[hitCounter]->neighbours.end()
199 == find(shvec[hitCounter]->neighbours.begin(),shvec[hitCounter]->neighbours.end(), shvec[nnIndex[neighbourCounter]]) ) {
200 //add neighbour to SuperHit
201 shvec[hitCounter]->neighbours.push_back(shvec[nnIndex[neighbourCounter]]);
202 shvec[hitCounter]->top++;
203
204 }
205 if ( shvec[nnIndex[neighbourCounter]]->neighbours.end()
206 == find(shvec[nnIndex[neighbourCounter]]->neighbours.begin(),shvec[nnIndex[neighbourCounter]]->neighbours.end(), shvec[hitCounter]) ) {
207 shvec[nnIndex[neighbourCounter]]->neighbours.push_back(shvec[hitCounter]);
208 shvec[nnIndex[neighbourCounter]]->top++;
209 }
210 }//all neighbours
211
212 } else if ( N_NearestNeighbours !=0 ) {
213
214 shvec[hitCounter]->top = N_NearestNeighbours;
215 for ( unsigned int neighbourCounter=0; neighbourCounter<N_NearestNeighbours ; neighbourCounter++ )
216 shvec[hitCounter]->neighbours.push_back( shvec[nnIndex[neighbourCounter]]);
217 }
218 } else if ( N_NearestNeighbours !=0 ) {
219
220 shvec[hitCounter]->top=N_NearestNeighbours;
221 for ( unsigned int neighbourCounter=0; neighbourCounter<N_NearestNeighbours; neighbourCounter++ )
222 shvec[hitCounter]->neighbours.push_back( shvec[nnIndex[neighbourCounter]]);
223 }
224 }//all hits
225
226 delete [] nnIndex;
227 delete [] distance;
228 annDeallocPts(dataPoints);
229 delete kdTree;
230 annClose();
231 }
232
233
234 //called by precalc2
235 void GridTransform2 ( CalorimeterHit* caloHit,float& radius, float& halfz, float& cellsize,float*X,bool isCalorimeter,int stave,int module,CellIDDecoder<CalorimeterHit>& id ) {
236 unsigned int layer;
237 double coss;
238 double sinn;
239
240 if ( !isCalorimeter ) { // !calorimeter -> tracker??
241 stave = id(caloHit)["S-1"];
242 module = id(caloHit)["M"]; //module is back to 6 or 1 for endcap?
243 layer = id(caloHit)["K-1"];
244
245 } else {
246
247 int stave2 = id(caloHit)["S-1"];
248 // The following calcs are done for the secondary stave the hit is assigned to, where it floats "freely"in space until now.
249 if ( stave2!=stave ) {
250 //define rotation matrix, hits are rotated to stave 0
251 coss = cos(M_PI_4*stave);
252 sinn = sin(M_PI_4*stave);
253 //rotate stave to position 0 (only y-axis)
254 float tmp_y = coss*caloHit->getPosition()[1] - sinn*caloHit->getPosition()[0];
255 //divide distance of hit from ecal face by sampling layer thickness->layer
256 layer = (unsigned int) ((tmp_y-radius)/myPGDB[myPGdb::ECAL1_BAR].sampling); // largest int number: 2.3->2, therefore n_sampl-1 in next line
257
258 if ( layer> myPGDB[myPGdb::ECAL1_BAR].n_sampl - 1 ) { // last layer of first ECAl structure
259 layer = myPGDB[myPGdb::ECAL1_BAR].n_sampl - 1
260 + (unsigned int)(( tmp_y-radius- myPGDB[myPGdb::ECAL1_BAR].n_sampl*
261 myPGDB[myPGdb::ECAL1_BAR].sampling)/myPGDB[myPGdb::ECAL2_BAR].sampling);
262 }
263
264 } else { // (stave2 == stave) original layer for hits in their original stave
265 layer = id(caloHit)["K-1"];
266 }
267 //Now we know of every hit's (new) layer
268 }
269 // out of if
270 coss = cos(M_PI_4*stave);
271 sinn = sin(M_PI_4*stave);
272
273 if ( module !=0 && module!=6 ) {
274 float tmp_x = coss*caloHit->getPosition()[0] + sinn*caloHit->getPosition()[1]; // rotate, find new y position
275 float tmp_y = radius + (layer+1)*cellsize; // center of layer position.
276 // y using the variable cellsize, the hits are equally spaced in 3d
277
278 //rotate back
279 X[0] = coss*tmp_x - sinn*tmp_y;
280 X[1] = coss*tmp_y + sinn*tmp_x;
281 X[2] = caloHit->getPosition()[2];
282
283 } else { // endcap
284 X[0] = caloHit->getPosition()[0];
285 X[1] = caloHit->getPosition()[1];
286 X[2] = (-1+module/3)*(halfz+layer*cellsize);
287 }
288 }
289
290
291 // searching for cluster cores
292 // wil be called by KIT
293 // secal1 = superhits in calo[4] with topological cleaning
294 // secal1 = superhits in calo[0] without topological cleaning
295 void FindCores2 (Shitvec2* secal1, Tmpclvec2* clusterVector , vector <PROTSEED2> * photonSeed2, unsigned int Nlevels, vector<float> mipstep, CoreCut Ccut,int debug )
296 {
297 if (debug >2 )
298 std::cout<<" ***** started main loop: searching for cluster cores in "<<secal1->size()<<" hits ***"<<std::endl;
299
300 typedef pair<int,int> test;
301 int psl_plun_global=0;
302
303 double Diagby2 = myPGDB[myPGdb::ECAL1_BAR].cell_size*sqrt(2.0)/2.0;
304 double barrel_rInner = myPGDB[myPGdb::ECAL1_BAR].r_inner; //inner ECAL radius
305 double barrel_zMax = myPGDB[myPGdb::ECAL1_CAP].z_inner; //maximum z barrel
306 double Xatf[3];
307 vector<Superhit2*> mySuperhits[Nlevels]; //lcio does not allow for more than 16 levels
308
309 for ( unsigned int hitCounter = 0; hitCounter<secal1->size(); hitCounter++ ) {
310 int myLevel = 0;
311
312 //search for cluster level thresholds: energy per hit > mipstep[level]
313 for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
314 if ( (*secal1)[hitCounter]->mip > mipstep[levelCounter] )
315 myLevel = levelCounter;
316 else
317 break;
318 }
319 //write back hits into correct levels
320 for ( int levelCounter=myLevel; levelCounter>=0; levelCounter-- ) {
321 mySuperhits[levelCounter].push_back((*secal1)[hitCounter]);
322 if ( debug >3 )
323 std::cout<<"DEBUG KITutil: * writing hits of level "<<levelCounter<<" to sorted superhits vector"<<std::endl;
324 }
325 }//all super hits
326 //temporary cluster vector
327 vector<Tmpcl2*> testvector;
328 //assigning neighbouring hits from superHitVector to temporary cluster for each level
329 for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
330 cluster5(&mySuperhits[levelCounter],&clusterVector[levelCounter],debug);
331 }
332
333 for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
334
335 if ( clusterVector[levelCounter].size()!=0 ) {
336 for ( unsigned int clusterCounter=0; clusterCounter<clusterVector[levelCounter].size(); clusterCounter++ ) {
337 clusterVector[levelCounter][clusterCounter]->calcCenter(); //energy weighted cluster center
338 if (debug >4 )
339 std::cout<<"calculated cluster center for level "<<levelCounter<<std::endl;
340 clusterVector[levelCounter][clusterCounter]->findInertia();//cluster start
341 if (debug >4 )
342 std::cout<<"calculated cluster start for level "<<levelCounter<<std::endl;
343
344 }
345 }
346 }//all levels
347
348 int myLevel = Nlevels - 1;
349 for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
350 if ( clusterVector[levelCounter].size()==0 ) {
351 myLevel = levelCounter - 1;
352 break;
353 }//no hits in this cluster level go one level down
354 myLevel = levelCounter;
355 }
356 psl_plun_global = myLevel;
357
358 vector < PROTSEED2 > photonSeed;
359
360 for (unsigned int clusterCounter=0; clusterCounter<clusterVector[0].size(); clusterCounter++ ) {
361 if ( clusterVector[0][clusterCounter]->hits.size()>Ccut.MinHit0 ) {
362 //yet another protseed
363 PROTSEED2 km;
364 km.cl = clusterVector[0][clusterCounter];
365 km.level = 0;
366 km.X[0] = 0.0;
367 km.X[1] = 0.0;
368 km.X[2] = 0.0;
369 km.active = true;
370 photonSeed.push_back(km);
371 }
372 }
373
374
375 for ( unsigned int seedCounter=0; seedCounter<photonSeed.size(); seedCounter++ ) {
376
377 Xatf[0] = 0.0;
378 Xatf[1] = 0.0;
379 Xatf[2] = 0.0;
380
381 //search for cluster closest to IP
382 LineCaloIntersectD2(photonSeed[seedCounter].cl->getCenter(), photonSeed[seedCounter].cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);//cluster center, cluster direction, inner barrel radius, inner z barrel,
383
384 if (debug >4 ) {
385 std::cout<<"calculated cluster - IP intersect for seed "<<seedCounter
386 <<" of "<<photonSeed.size()<<std::endl;
387 std::cout<<"- cluster center: "<<photonSeed[seedCounter].cl->getCenter()<<std::endl;
388 std::cout<<"- cluster direction: "<<photonSeed[seedCounter].cl->direction<<std::endl;
389 std::cout<<"- ECAL1 inner barrel: "<<barrel_rInner<<std::endl;
390 std::cout<<"- ECAL1 max z: "<<barrel_zMax<<std::endl;
391 }
392
393 //assign this cluster shifted by its depth in ECAL to photonSeed
394 photonSeed[seedCounter].X[0] = Xatf[0];
395 photonSeed[seedCounter].X[1] = Xatf[1];
396 photonSeed[seedCounter].X[2] = Xatf[2];
397 photonSeed2->push_back(photonSeed[seedCounter]);
398
399 if (debug >5 )
400 std::cout<<"assign shifted hits to new photonSeed of size"<<photonSeed2->size()<<std::endl;
401 }//all seed entries
402
403
404 //search for splitted clusters
405 for ( int levelCounter=1; levelCounter<=myLevel; levelCounter++) {
406
407 for ( unsigned int seed2Counter=0; seed2Counter<photonSeed2->size(); seed2Counter++) {
408 if ( (*photonSeed2)[seed2Counter].active ) {
409
410 Tmpclvec2 tempCluster; //clusters inside clusters
411 ClusterInCluster2((*photonSeed2)[seed2Counter].cl,clusterVector[levelCounter],tempCluster, debug);
412
413 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
414 //FIXME: CAREFULL photonSeed.Size increases 1 for each call of ClusterInCluster2 ! -> Infinite loop !
415 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
416 if (debug >4 )
417 std::cout<<"after ClusterInCluster2: photonSeed size = "<<photonSeed2->size()<<std::endl;
418
419 if ( tempCluster.size()<=1 ) {
420 if (debug >4 )
421 std::cout<<"did not find cluster in cluster "<<std::endl;
422 // 1-in-1 nothing to do
423 } else { //found cluster tempCluster in cluster photonSeed2
424 // Split the split in two categories
425 // clear one (larger distance) and not so clear
426
427 if (debug >4 )
428 std::cout<<"level "<<levelCounter<<": found cluster in cluster "<<std::endl;
429
430 vector < int > largerDistanceCluster;
431 vector < int > smallerDistanceCluster;
432 //check included clusters
433 for ( unsigned int clusterCounter1=0; clusterCounter1<tempCluster.size(); clusterCounter1++ ) {
434 for ( unsigned int clusterCounter2=0; clusterCounter2<tempCluster.size(); clusterCounter2++ ) {
435 if ( clusterCounter1!=clusterCounter2 ) {
436 if ( tempCluster[clusterCounter1]->hits.size()>=Ccut.MinHitSplit //more than minimum numgber of hits in included cluster
437 && tempCluster[clusterCounter2]->hits.size()>=Ccut.MinHitSplit ) {
438
439 double XIP[3];
440 XIP[0] = 0.0;
441 XIP[1] = 0.0;
442 XIP[2] = 0.0;
443
444
445 if ( (LinePointDistance2(XIP,(*photonSeed2)[seed2Counter].cl->center,tempCluster[clusterCounter2]->center) > Diagby2
446 || LinePointDistance2(XIP,(*photonSeed2)[seed2Counter].cl->center,tempCluster[clusterCounter1]->center) > Diagby2 ) ) {//distance of tempCluster to photonSeed2 > cellSize *sqrt(2) / 2 ->clear cluster split
447 if ( find(largerDistanceCluster.begin(),largerDistanceCluster.end(),(int)clusterCounter2)==largerDistanceCluster.end() )
448
449 largerDistanceCluster.push_back(clusterCounter2);
450
451 if ( find(largerDistanceCluster.begin(),largerDistanceCluster.end(),(int)clusterCounter1)==largerDistanceCluster.end() )
452
453 largerDistanceCluster.push_back(clusterCounter1);
454 } else { //small distance between clusters -> split unclear
455 if (debug >5 )
456 std::cout<<"not so clearly seperated small distant cluster -> Split unclear"<<std::endl;
457
458 if (find(smallerDistanceCluster.begin(),smallerDistanceCluster.end(),(int)clusterCounter2)==smallerDistanceCluster.end() )
459
460 smallerDistanceCluster.push_back(clusterCounter2);
461
462 if (find(smallerDistanceCluster.begin(),smallerDistanceCluster.end(),(int)clusterCounter1)==smallerDistanceCluster.end() )
463
464 smallerDistanceCluster.push_back(clusterCounter1);
465 }
466 }//found cluster in cluster
467 }//clusterCounters don't agree
468 }//all tempCluster entries, first loop
469 }//all tempCluster entries, second loop
470
471 if ( largerDistanceCluster.size()!=0 ) { //if clearly seperated cluster exists
472
473 if (debug >5 )
474 std::cout<<"clearly seperated large distant cluster "<<std::endl;
475
476 if (debug >6 )
477 std::cout<<"for clearly seperated large distant clusters: photonSeed size = "<<photonSeed2->size()<<std::endl;
478 int clusterSize1 = 0;
479 int clusterSize3 = (*photonSeed2)[seed2Counter].cl->hits.size();
480 for ( unsigned int vectorEntry=0; vectorEntry<largerDistanceCluster.size(); vectorEntry++ )
481 clusterSize1 += tempCluster[largerDistanceCluster[vectorEntry]]->hits.size();
482
483 double ratio = ((double)clusterSize1) / ((double)clusterSize3);
484 (*photonSeed2)[seed2Counter].active = false;
485
486 if ( ratio > Ccut.rCut ) { //require minimum ratio of cluster entries to supress fluctuations
487 for ( unsigned int clusterEntry=0; clusterEntry<smallerDistanceCluster.size(); clusterEntry++ ) { //have a look at the not so clear splitted clusters
488
489 if (debug >5 )
490 std::cout<<"cluster ratio larger than"<<Ccut.rCut<<" hits, loop over small distance clusters"<<std::endl;
491
492 Xatf[0] = 0.0;
493 Xatf[1] = 0.0;
494 Xatf[2] = 0.0;
495
496 PROTSEED2 clusterSeed;
497
498 clusterSeed.cl = tempCluster[smallerDistanceCluster[clusterEntry]];
499 clusterSeed.level = levelCounter;
500 LineCaloIntersectD2(tempCluster[smallerDistanceCluster[clusterEntry]]->getCenter(),tempCluster[smallerDistanceCluster[clusterEntry]]->direction,barrel_rInner,barrel_zMax,Xatf,debug);
501 clusterSeed.X[0] = Xatf[0];
502 clusterSeed.X[1] = Xatf[1];
503 clusterSeed.X[2] = Xatf[2];
504 clusterSeed.active = true;
505 bool split = true;
506
507 for ( unsigned int seed2Counter=0; seed2Counter<photonSeed2->size(); seed2Counter++ ) {
508 if ( (*photonSeed2)[seed2Counter].active ) {
509 if ( (*photonSeed2)[seed2Counter].cl==clusterSeed.cl )
510 split=false;
511 if (debug >5 )
512 std::cout<<"do not split cluster"<<std::endl;
513 }
514 }
515 if (split) {
516 photonSeed2->push_back(clusterSeed);
517 if (debug >5 )
518 std::cout<<"splitting cluster"<<std::endl;
519 }
520 }
521 }//ratio > ratio cut
522
523
524 for ( unsigned int clusterEntry=0; clusterEntry<largerDistanceCluster.size(); clusterEntry++ ) {
525
526 Xatf[0] = 0.0;
527 Xatf[1] = 0.0;
528 Xatf[2] = 0.0;
529
530 PROTSEED2 clusterSeed;
531
532 if (debug >4 )
533 std::cout<<"cluster larger than "<<Ccut.rCut<<" hits, loop over large distance clusters"<<std::endl;
534
535 clusterSeed.cl = tempCluster[largerDistanceCluster[clusterEntry]];
536 clusterSeed.level = myLevel;
537 LineCaloIntersectD2(tempCluster[largerDistanceCluster[clusterEntry]]->getCenter(),tempCluster[largerDistanceCluster[clusterEntry]]->direction,barrel_rInner,barrel_zMax,Xatf,debug);
538 clusterSeed.X[0] = Xatf[0];
539 clusterSeed.X[1] = Xatf[1];
540 clusterSeed.X[2] = Xatf[2];
541 clusterSeed.active = true;
542 bool split = true;
543 for ( unsigned int seed2Counter=0; seed2Counter<photonSeed2->size(); seed2Counter++ ) {
544 if ( (*photonSeed2)[seed2Counter].active ) {
545 if ( (*photonSeed2)[seed2Counter].cl==clusterSeed.cl )
546 split = false;
547 }
548 }
549 if ( split )
550 photonSeed2->push_back(clusterSeed);
551 }
552 }//found split clusterSeed at large distance
553 if ( debug >6 ) std::cout<<"exit large distance cluster"<<std::endl;
554 }//found cluster in cluster
555 if ( debug >6 ) std::cout<<"exit cluster in cluster"<<std::endl;
556 }// if active
557 if ( debug >6 ) std::cout<<"exit photonseed.active"<<std::endl;
558 if ( debug >6 ) std::cout<<"for seed "<<seed2Counter<<" of "<< photonSeed2->size()<<std::endl;
559 }//in photonSeed2
560 if ( debug >6 ) std::cout<<"exit loop over photonseed"<<std::endl;
561 }//all hit energy threshold levels
562 if ( debug >6 ) std::cout<<"exit loop over threhold limits"<<std::endl;
563
564 vector <test> mergeVector;
565
566 for ( unsigned int hitCounter=0; hitCounter<photonSeed2->size(); hitCounter++ ) {
567 if( (*photonSeed2)[hitCounter].active ) {
568
569 double Xa[3];
570 ModuleNormal2((*photonSeed2)[hitCounter].cl->center,barrel_zMax,Xa); //shift hits in barrel
571 if (debug >3 )
572 std::cout<<" start merging: shifting barrel hits"<<std::endl;
573
574 //double loop over all hits?
575 for ( unsigned int hitCounter2=0; hitCounter2<photonSeed2->size(); hitCounter2++) {
576
577 //if hitcounters do not agree, get difference between these hits
578 if ( hitCounter2!=hitCounter && (*photonSeed2)[hitCounter2].active ) {
579
580 double dir_diff[3];
581 dir_diff[0] = (*photonSeed2)[hitCounter].cl->center[0]
582 - (*photonSeed2)[hitCounter2].cl->center[0];
583 dir_diff[1] = (*photonSeed2)[hitCounter].cl->center[1]
584 - (*photonSeed2)[hitCounter2].cl->center[1];
585 dir_diff[2] = (*photonSeed2)[hitCounter].cl->center[2]
586 - (*photonSeed2)[hitCounter2].cl->center[2];
587
588 //MERGE CONDITION for clusters:
589 //D_cl_cl2 returns smallest distance between cluster1 & cluster2
590 //Dot2 returns angle between clusters center & split candidate ??
591 if ( D_cl_cl2((*photonSeed2)[hitCounter].cl,(*photonSeed2)[hitCounter2].cl) < Ccut.distCut //minimum distance between clusters < distance cut
592 && fabs( Dot2((*photonSeed2)[hitCounter2].cl->center,dir_diff)) > Ccut.cosCut) { //minimum angle between clusters
593
594 test cluster1,cluster2; //pair of cluster1 & cluster2
595 cluster1.first = hitCounter;
596 cluster1.second = hitCounter2;
597 cluster2.first = hitCounter2;
598 cluster2.second = hitCounter;
599 //write result in test vector
600 if ( find(mergeVector.begin(),mergeVector.end(),cluster1)==mergeVector.end()
601 && find(mergeVector.begin(),mergeVector.end(),cluster2)==mergeVector.end() ) {
602 mergeVector.push_back(cluster1);
603
604 if (debug >3 )
605 std::cout<<" merging clusters"<<std::endl;
606 }
607 }//if merge condition
608 }//hitCounters do not agree
609 }//all photonSeed entries
610 }//photonSeed active
611 }//all photonSeed entries
612
613
614 // if(testVector.size()!=0) //if testvector not empty -> differences in hits
615 if ( mergeVector.size()!=0) {
616
617 int benefit[mergeVector.size()];
618 for ( unsigned int vectorCounter=0; vectorCounter<mergeVector.size(); vectorCounter++) {
619 benefit[vectorCounter] = 0;
620 }
621 for ( unsigned int vectorCounter=0; vectorCounter<mergeVector.size(); vectorCounter++ ) {
622 if ( benefit[vectorCounter]==0 ) {
623
624 vector <int> trinityVector; //holds entries from merged clusters
625 trinityVector.push_back(mergeVector[vectorCounter].first);
626 trinityVector.push_back(mergeVector[vectorCounter].second);
627 benefit[vectorCounter] = 1;
628 for(unsigned int vectorCounter2=0; vectorCounter2<mergeVector.size(); vectorCounter2++ ) {
629
630 if ( (find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].first) != trinityVector.end()
631 && find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].second) == trinityVector.end() )
632 || (find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].second) != trinityVector.end()
633 && find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].first) == trinityVector.end() ) ) {
634
635 if ( find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].first) == trinityVector.end() ) {
636 if ( benefit[vectorCounter2] !=1 ) { //->vectorCounter1 != vectorCounter
637
638 benefit[vectorCounter2] = 1;
639 if ( (*photonSeed2)[mergeVector[vectorCounter2].first].active )
640 trinityVector.push_back(mergeVector[vectorCounter2].first);
641 }
642
643 } else if ( benefit[vectorCounter2] != 1 ) {
644 if( (*photonSeed2)[mergeVector[vectorCounter2].second].active)
645 trinityVector.push_back(mergeVector[vectorCounter2].second);
646 benefit[vectorCounter2] = 1;
647 }
648 }
649 }
650
651
652 int levelCheck = 1;
653 for ( unsigned int vectorEntry=1; vectorEntry<trinityVector.size(); vectorEntry++) {
654 if( (*photonSeed2)[trinityVector[0]].level != (*photonSeed2)[trinityVector[vectorEntry]].level )
655 levelCheck++;
656 }
657
658 if ( levelCheck==1 ) { //all cluster entries in same level
659
660 Tmpcl2* cl = new Tmpcl2();
661 for ( unsigned int vectorEntry=0; vectorEntry<trinityVector.size(); vectorEntry++ ) {
662 if ( (*photonSeed2)[trinityVector[vectorEntry]].active ) {
663 for ( unsigned int hitCounter=0; hitCounter<(*photonSeed2)[trinityVector[vectorEntry]].cl->hits.size(); hitCounter++ ) {
664 cl->hits.push_back((*photonSeed2)[trinityVector[vectorEntry]].cl->hits[hitCounter]);
665 }
666 (*photonSeed2)[trinityVector[vectorEntry]].active = false;
667 }
668 }
669 if ( cl->hits.size()!=0 ) {//found merged cluster
670
671 cl->calcCenter();
672 if (debug >4 )
673 std::cout<<" calculate center of merged cluster"<<std::endl;
674 cl->findInertia();
675 if (debug >4 )
676 std::cout<<" calculate start of merged cluster"<<std::endl;
677 clusterVector[(*photonSeed2)[trinityVector[0]].level].push_back(cl);
678
679 Xatf[0] = 0.0;
680 Xatf[1] = 0.0;
681 Xatf[2] = 0.0;
682
683 PROTSEED2 clusterSeed;
684 clusterSeed.cl = cl;
685 clusterSeed.level = (*photonSeed2)[trinityVector[0]].level;
686 LineCaloIntersectD2(cl->getCenter(),cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);
687 clusterSeed.X[0] = Xatf[0];
688 clusterSeed.X[1] = Xatf[1];
689 clusterSeed.X[2] = Xatf[2];
690 clusterSeed.active = true;
691
692 bool split = true;
693 for ( unsigned int seedEntry=0; seedEntry<photonSeed2->size(); seedEntry++ ) {
694 if ( (*photonSeed2)[seedEntry].active ) {
695 if( (*photonSeed2)[seedEntry].cl==clusterSeed.cl)
696 split = false;
697 }
698 }
699 if(split)
700 photonSeed2->push_back(clusterSeed);
701
702 } else {//no merged clusters
703 delete cl;
704 }
705 } else {//levels do not agree
706
707 if (debug >5 )
708 std::cout<<" clusters to merge have different levels -> create new vector"<<std::endl;
709 vector <int > newVector;
710 for ( unsigned int j=0; j<trinityVector.size(); j++ ) {//write trinityVector into newVector
711 if ( find(newVector.begin(),newVector.end(),(*photonSeed2)[trinityVector[j]].level)==newVector.end() )
712 newVector.push_back((*photonSeed2)[trinityVector[j]].level);
713 }
714
715 if ( trinityVector.size()==2 ) {//should always be, else trinityVector empty??
716 if ( (*photonSeed2)[trinityVector[0]].cl->getEnergy() >=(*photonSeed2)[trinityVector[1]].cl->getEnergy() ) {
717 (*photonSeed2)[trinityVector[1]].active = false;
718 } else {
719 (*photonSeed2)[trinityVector[0]].active = false;
720 }
721 } else {//trinityVector.size !=2
722
723 vector <int> selectedHits;
724 int sumedUp[trinityVector.size()];
725 for ( unsigned int trinityEntry=0; trinityEntry<trinityVector.size(); trinityEntry++ ) {
726 sumedUp[trinityEntry] = 0;
727 }
728 for ( unsigned int newVectorEntry=0; newVectorEntry<newVector.size(); newVectorEntry++ ) {
729 vector<int> sumUpVector;
730 for ( unsigned int trinityEntry=0; trinityEntry<trinityVector.size(); trinityEntry++) {
731 if ( sumedUp[trinityEntry]==0 ) {
732 if ( newVector[newVectorEntry]==(*photonSeed2)[trinityVector[trinityEntry]].level) {
733 sumUpVector.push_back(trinityVector[trinityEntry]);
734 sumedUp[trinityEntry] = 1;
735 }
736 }
737 }
738
739 if ( sumUpVector.size()>1 ) {
740
741 Tmpcl2* cl = new Tmpcl2();
742 for ( unsigned int sumUpVectorEntry=0; sumUpVectorEntry<sumUpVector.size(); sumUpVectorEntry++ ) {
743 (*photonSeed2)[sumUpVector[sumUpVectorEntry]].active = false;
744
745 for ( unsigned int hitCounter=0; hitCounter<(*photonSeed2)[sumUpVector[sumUpVectorEntry]].cl->hits.size(); hitCounter++ )
746 cl->hits.push_back((*photonSeed2)[sumUpVector[sumUpVectorEntry]].cl->hits[hitCounter]);
747 if (debug >5 )
748 std::cout<<" writing new, merged cluster"<<std::endl;
749 }
750
751 cl->calcCenter();
752 cl->findInertia();
753 clusterVector[newVector[newVectorEntry]].push_back(cl);
754
755
756 Xatf[0] = 0.0;
757 Xatf[1] = 0.0;
758 Xatf[2] = 0.0;
759
760 PROTSEED2 clusterSeed;
761
762 clusterSeed.cl = cl;
763 clusterSeed.level = (*photonSeed2)[newVector[newVectorEntry]].level;
764 LineCaloIntersectD2(cl->getCenter(),cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);
765 clusterSeed.X[0] = Xatf[0];
766 clusterSeed.X[1] = Xatf[1];
767 clusterSeed.X[2] = Xatf[2];
768 clusterSeed.active = true;
769
770 bool split = true;
771 for ( unsigned int im=0; im<photonSeed2->size(); im++ ) {
772 if ( (*photonSeed2)[im].active ) {
773 if ( (*photonSeed2)[im].cl==clusterSeed.cl )
774 split=false;
775 }
776 }
777 if ( split ) {
778 photonSeed2->push_back(clusterSeed);
779 selectedHits.push_back(photonSeed2->size()-1);
780 } else {
781 delete cl;
782 }
783 } else {
784 selectedHits.push_back(sumUpVector[0]);
785 }
786 }//all newVector entries
787
788 if ( selectedHits.size()>1 ) {
789 if (debug >5 )
790 std::cout<<" pairing merged cluster"<<std::endl;
791 typedef pair <double,int> DpairI;
792 vector <DpairI> selectedPairsVector;
793 for ( unsigned int j=0; j<selectedHits.size(); j++ ) {
794
795 DpairI tmp;
796 tmp.first = (*photonSeed2)[selectedHits[j]].cl->getEnergy();
797 tmp.second = selectedHits[j];
798 selectedPairsVector.push_back(tmp);
799 }
800
801 sort(selectedPairsVector.begin(),selectedPairsVector.end());
802
803 for ( unsigned int j=0; j<selectedHits.size()-1; j++ ) {
804 (*photonSeed2)[selectedPairsVector[j].second].active=false;
805 }
806 }//found selected Hits
807 }//trinityVector.size !=2
808 }//levels do not agree
809 }//benefit==0
810 } // mergeVector entries
811 }//merge clusters
812 if (debug >2 )
813 std::cout<<" ***** FindCores2 is done ***"<<std::endl;
814 }//FindCores2
815
816
817
818 //called by FindCores2, asigning neighbouring hits to cluster
819 void cluster5 ( vector<Superhit2*>* superHitVector, vector<Tmpcl2*>* clusterVector, int debug )
820 {
821 if ( (*superHitVector).size()!=0 ) {
822 if ( debug > 5 )
823 std::cout<<"DEBUG KITutil: * clustering neighbouring hits of same level "<<std::endl;
824
825 unsigned int superHitVectorSize = (*superHitVector).size();
826
827 for ( unsigned int hitCounter=0; hitCounter<superHitVectorSize; hitCounter++ ) {
828 (*superHitVector)[hitCounter]->is_assigned = false;
829 //fill superHitVector into temporary clusterVector
830 if ( (*superHitVector)[hitCounter]->top==0 ) {
831 Tmpcl2* tempCluster = new Tmpcl2();
832 tempCluster->hits.push_back((*superHitVector)[hitCounter]);
833 (*superHitVector)[hitCounter]->is_assigned = true;
834 clusterVector->push_back(tempCluster);
835 continue;
836
837 } else if ( (*superHitVector)[hitCounter]->is_assigned==false ) {
838 vector<Superhit2*> sshv;
839 stack <Superhit2*> superhitStack;
840 sshv.push_back((*superHitVector)[hitCounter]);
841
842 //assigning neighbouring hits
843 for ( unsigned int neighbourCounter=0; neighbourCounter<(*superHitVector)[hitCounter]->neighbours.size(); neighbourCounter++ ) {
844 if ( !((*superHitVector)[hitCounter]->neighbours[neighbourCounter]->is_assigned)
845 && superHitVector->end() != find(superHitVector->begin(),superHitVector->end(),(*superHitVector)[hitCounter]->neighbours[neighbourCounter]) ) {
846
847 sshv.push_back((*superHitVector)[hitCounter]->neighbours[neighbourCounter]);
848 (*superHitVector)[hitCounter]->neighbours[neighbourCounter]->is_assigned = true;
849 superhitStack.push((*superHitVector)[hitCounter]->neighbours[neighbourCounter]);
850 }
851 }//all neighbours
852
853 while ( !superhitStack.empty() ) {
854
855 Superhit2* superHit = superhitStack.top();
856 superhitStack.pop();
857
858 for ( unsigned int neighbourCounter=0; neighbourCounter<superHit->neighbours.size(); neighbourCounter++ ) {
859 if ( superHit->neighbours[neighbourCounter]!=0 ) {
860 if ( ! superHit->neighbours[neighbourCounter]->is_assigned &&
861 superHitVector->end() != find(superHitVector->begin(),superHitVector->end(),superHit->neighbours[neighbourCounter]) )
862 if ( sshv.end() == find(sshv.begin(),sshv.end(),superHit->neighbours[neighbourCounter]) ) {
863 sshv.push_back(superHit->neighbours[neighbourCounter]);
864 superHit->neighbours[neighbourCounter]->is_assigned=true;
865
866 superhitStack.push(superHit->neighbours[neighbourCounter]);
867 }
868 } // neighbours not 0
869 } // all neighbours
870 if ( debug >5 )
871 std::cout<<"DEBUG KITutil: * found neighbours "<<superHit->neighbours.size()<<std::endl;
872 } // stack not empty
873
874 if ( sshv.size()!=0 ) {
875 Tmpcl2* tempCluster = new Tmpcl2();
876 for ( unsigned int sshvCounter=0; sshvCounter<sshv.size(); sshvCounter++ ) {
877 tempCluster->hits.push_back(sshv[sshvCounter]);
878 }
879 clusterVector->push_back(tempCluster);
880 }
881 } // no superhit assigned
882 } // over all hits
883
884 for ( unsigned int superHitVectorCounter=0; superHitVectorCounter<superHitVector->size(); superHitVectorCounter++ )
885 if ( (*superHitVector)[superHitVectorCounter]->is_assigned == false ) {
886 Tmpcl2* tempCluster = new Tmpcl2();
887 tempCluster->hits.push_back((*superHitVector)[superHitVectorCounter]);
888 clusterVector->push_back(tempCluster);
889 if ( debug >5 )
890 std::cout<<"DEBUG KITutil: * found cluster "<<tempCluster<<std::endl;
891 }
892 }//superhit vector size not 0
893 }//cluster5
894
895
896 Superhit2::~Superhit2 () {}
897
898 Superhit2::Superhit2 ( float E, CalorimeterHit* chit_in ) {
899 chit = chit_in;
900 for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
901 point[positionCounter] = 0.0;
902 pointt[positionCounter] = chit_in->getPosition()[positionCounter];
903 }
904
905 mip = chit->getEnergy() / E; //mip calibration
906 connect = false;
907 is_assigned = false;
908 mipE = E;
909 top = 0;
910 cl = 0;
911
912 S = 0; //stave
913 M = 0; //module
914 K = 0; //layer
915 isCalorimeter = 0; // ecal=1, hcal=2, else=0
916
917 }
918
919 mySuperhit::~mySuperhit () {}
920
921 mySuperhit::mySuperhit ( float E, CalorimeterHit* chit_in ) {
922 chit = chit_in;
923 for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
924 point[positionCounter] = 0.0;
925 pointt[positionCounter] = chit_in->getPosition()[positionCounter];
926 }
927
928 mip = chit->getEnergy() / E; //mip calibration
929 connect = false;
930 is_assigned = false;
931 mipE = E;
932 top = 0;
933 cl = 0;
934
935 S = 0; //stave
936 M = 0; //module
937 K = 0; //layer
938 }
939
940
941 double* myTmpcl:: getCenter()
942 {
943 return center;
944 }
945
946
947 myTmpcl::myTmpcl(){
948 energy=0.0;
949 }
950 myTmpcl::~myTmpcl(){}
951
952 double myTmpcl::getEnergy(){
953
954 energy =0.0;
955 for ( unsigned int i=0; i<hits.size(); i++ ) {
956 energy+=hits[i]->chit->getEnergy();
957 }
958 return energy;
959 }
960
961 //calculating energy weighted center of cluster
962 void myTmpcl::calcCenter()
963 {
964 center[0] = 0.0;
965 center[1] = 0.0;
966 center[2] = 0.0;
967 double energySum = 0.0;
968 double energy = 0.0;
969 for ( unsigned int hitCounter=0; hitCounter<hits.size(); hitCounter++ ) {
970 energy = hits[hitCounter]->chit->getEnergy();
971
972 center[0] += hits[hitCounter]->pointt[0] * energy;
973 center[1] += hits[hitCounter]->pointt[1] * energy;
974 center[2] += hits[hitCounter]->pointt[2] * energy;
975 energySum += energy;
976 }
977 center[0] = center[0] / energySum;
978 center[1] = center[1] / energySum;
979 center[2] = center[2] / energySum;
980 }
981
982
983 //finding cluster start
984 void myTmpcl::findInertia()
985 {
986
987 double aIne[3][3];
988 for ( int xPosCounter(0); xPosCounter<3; ++xPosCounter )
989 for ( int yPosCounter(0); yPosCounter < 3; ++yPosCounter ) {
990 aIne[xPosCounter][yPosCounter] = 0.0;
991 }
992
993
994 for ( unsigned int hitCounter=0; hitCounter<hits.size(); hitCounter++) {
995
996 //distance hit - cluster center
997 double distX = hits[hitCounter]->pointt[0] - center[0];
998 double distY = hits[hitCounter]->pointt[1] - center[1];
999 double distZ = hits[hitCounter]->pointt[2] - center[2];
1000
1001 double eHit = hits[hitCounter]->chit->getEnergy();
1002 eHit = 1.0;
1003
1004 aIne[0][0] += eHit * (distY*distY + distZ*distZ);
1005 aIne[1][1] += eHit * (distX*distX + distZ*distZ);
1006 aIne[2][2] += eHit * (distX*distX + distY*distY);
1007
1008 aIne[0][1] -= eHit * distX*distY;
1009 aIne[0][2] -= eHit * distX*distZ;
1010 aIne[1][2] -= eHit * distY*distZ;
1011 }
1012
1013 for (int i(0); i<2; ++i ) {
1014 for (int j=i+1; j<3; ++j ) {
1015 aIne[j][i] = aIne[i][j];
1016 }
1017 }
1018
1019
1020 gsl_matrix_view aMatrix = gsl_matrix_view_array((double*)aIne,3,3);
1021 gsl_vector* aVector = gsl_vector_alloc(3);
1022 gsl_matrix* aEigenVec = gsl_matrix_alloc(3,3);
1023 gsl_eigen_symmv_workspace* wa = gsl_eigen_symmv_alloc(3);
1024 gsl_eigen_symmv(&aMatrix.matrix,aVector,aEigenVec,wa);
1025 gsl_eigen_symmv_free(wa);
1026 gsl_eigen_symmv_sort(aVector,aEigenVec,GSL_EIGEN_SORT_ABS_ASC);
1027
1028 double norm = 0.0;
1029 for ( int i(0); i < 3; i++ ) {
1030 inteigen[i] = gsl_vector_get(aVector,i);
1031 norm += inteigen[i] * inteigen[i];
1032 for (int j(0); j<3; j++ ) {
1033 inteigenvec[i+3*j] = gsl_matrix_get(aEigenVec,i,j);
1034 }
1035 }
1036 norm = sqrt(norm);
1037
1038
1039 for ( int i(0); i<3; i++ )
1040 inteigen[i] = inteigen[i]/norm;
1041
1042 direction[0] = inteigenvec[0];
1043 direction[1] = inteigenvec[1];
1044 direction[2] = inteigenvec[2];
1045
1046 gsl_vector_free(aVector);
1047 gsl_matrix_free(aEigenVec);
1048 }//find inertia
1049
1050
1051 double* Tmpcl2:: getCenter()
1052 {
1053 return center;
1054 }
1055
1056
1057 Tmpcl2::Tmpcl2(){
1058 energy=0.0;
1059 }
1060 Tmpcl2::~Tmpcl2(){}
1061
1062 double Tmpcl2::getEnergy(){
1063
1064 energy =0.0;
1065 for ( unsigned int i=0; i<hits.size(); i++ ) {
1066 energy+=hits[i]->chit->getEnergy();
1067 }
1068 return energy;
1069 }
1070
1071 //calculating energy weighted center of cluster
1072 void Tmpcl2::calcCenter()
1073 {
1074 center[0] = 0.0;
1075 center[1] = 0.0;
1076 center[2] = 0.0;
1077 double energySum = 0.0;
1078 double energy = 0.0;
1079 for ( unsigned int hitCounter=0; hitCounter<hits.size(); hitCounter++ ) {
1080 energy = hits[hitCounter]->chit->getEnergy();
1081
1082 center[0] += hits[hitCounter]->pointt[0] * energy;
1083 center[1] += hits[hitCounter]->pointt[1] * energy;
1084 center[2] += hits[hitCounter]->pointt[2] * energy;
1085 energySum += energy;
1086 }
1087 center[0] = center[0] / energySum;
1088 center[1] = center[1] / energySum;
1089 center[2] = center[2] / energySum;
1090 }
1091
1092
1093 //finding cluster start
1094 void Tmpcl2::findInertia()
1095 {
1096
1097 double aIne[3][3];
1098 for ( int xPosCounter(0); xPosCounter<3; ++xPosCounter )
1099 for ( int yPosCounter(0); yPosCounter < 3; ++yPosCounter ) {
1100 aIne[xPosCounter][yPosCounter] = 0.0;
1101 }
1102
1103
1104 for ( unsigned int hitCounter=0; hitCounter<hits.size(); hitCounter++) {
1105
1106 //distance hit - cluster center
1107 double distX = hits[hitCounter]->pointt[0] - center[0];
1108 double distY = hits[hitCounter]->pointt[1] - center[1];
1109 double distZ = hits[hitCounter]->pointt[2] - center[2];
1110
1111 double eHit = hits[hitCounter]->chit->getEnergy();
1112 eHit = 1.0;
1113
1114 aIne[0][0] += eHit * (distY*distY + distZ*distZ);
1115 aIne[1][1] += eHit * (distX*distX + distZ*distZ);
1116 aIne[2][2] += eHit * (distX*distX + distY*distY);
1117
1118 aIne[0][1] -= eHit * distX*distY;
1119 aIne[0][2] -= eHit * distX*distZ;
1120 aIne[1][2] -= eHit * distY*distZ;
1121 }
1122
1123 for (int i(0); i<2; ++i ) {
1124 for (int j=i+1; j<3; ++j ) {
1125 aIne[j][i] = aIne[i][j];
1126 }
1127 }
1128
1129
1130 gsl_matrix_view aMatrix = gsl_matrix_view_array((double*)aIne,3,3);
1131 gsl_vector* aVector = gsl_vector_alloc(3);
1132 gsl_matrix* aEigenVec = gsl_matrix_alloc(3,3);
1133 gsl_eigen_symmv_workspace* wa = gsl_eigen_symmv_alloc(3);
1134 gsl_eigen_symmv(&aMatrix.matrix,aVector,aEigenVec,wa);
1135 gsl_eigen_symmv_free(wa);
1136 gsl_eigen_symmv_sort(aVector,aEigenVec,GSL_EIGEN_SORT_ABS_ASC);
1137
1138 double norm = 0.0;
1139 for ( int i(0); i < 3; i++ ) {
1140 inteigen[i] = gsl_vector_get(aVector,i);
1141 norm += inteigen[i] * inteigen[i];
1142 for (int j(0); j<3; j++ ) {
1143 inteigenvec[i+3*j] = gsl_matrix_get(aEigenVec,i,j);
1144 }
1145 }
1146 norm = sqrt(norm);
1147
1148
1149 for ( int i(0); i<3; i++ )
1150 inteigen[i] = inteigen[i]/norm;
1151
1152 direction[0] = inteigenvec[0];
1153 direction[1] = inteigenvec[1];
1154 direction[2] = inteigenvec[2];
1155
1156 gsl_vector_free(aVector);
1157 gsl_matrix_free(aEigenVec);
1158 }//find inertia
1159
1160
1161 //find shower start in first ECAL layer?
1162 //where does 400 come from?
1163 void LineCaloIntersectD2( double* X1,double* dir,double& r_inner,double& zmax,double* X,int debug)
1164 { //X1: cluster center; dir: cluster direction; r_inner: inner barrel radius, zmax: maximum barrel z-value; X: 0
1165
1166 if(debug >4)
1167 std::cout<<"---- FindCores2 called LineCaloIntersect -----"<<std::endl;
1168
1169 double X2[3]; //cluster center + direction
1170 //distinguish between upper and lower barrel half
1171 X2[0] = X1[0] + dir[0] * 400.0;
1172 X2[1] = X1[1] + dir[1] * 400.0;
1173 X2[2] = X1[2] + dir[2] * 400.0;
1174
1175 LineCaloIntersect2(X2,X1,r_inner,zmax,X,debug);
1176
1177 if ( X[0]==0.0 ) {
1178
1179 X2[0] = X1[0] - dir[0] * 400.0;
1180 X2[1] = X1[1] - dir[1] * 400.0;
1181 X2[2] = X1[2] - dir[2] * 400.0;
1182
1183 LineCaloIntersect2(X2,X1,r_inner,zmax,X,debug);
1184
1185 }
1186 if(debug >4)
1187 std::cout<<"---- Finished LineCaloIntersect ----"<<std::endl;
1188
1189 }
1190
1191
1192 // CAREFULL hardcoded numbers!!
1193 // written for LDC00_02 -> assumes 8fold detector geometry + 2 endcaps!
1194 //searches for cluster closest to IP
1195 void LineCaloIntersect2(double* X1, double* X2,double& r_inner,double& zmax,double* X, int debug)
1196 {
1197 double n[10][3];
1198
1199 double sqrt2_half = M_SQRT2/2.0; //= sqrt(2)/2 = 0.71
1200
1201 //n[stave][xyz] ->stave geometry
1202 n[0][0] = 0.0;
1203 n[0][1] = 1.0;
1204 n[0][2] = 0.0;
1205
1206 n[1][0] = sqrt2_half;
1207 n[1][1] = sqrt2_half;
1208 n[1][2] = 0.0;
1209
1210 n[2][0] = 1.0;
1211 n[2][1] = 0.0;
1212 n[2][2] = 0.0;
1213
1214 n[3][0] = -1.0;
1215 n[3][1] = 0.0;
1216 n[3][2] = 0.0;
1217
1218 n[4][0] = 0.0;
1219 n[4][1] = -1.0;
1220 n[4][2] = 0.0;
1221
1222 n[5][0] = -sqrt2_half;
1223 n[5][1] = -sqrt2_half;
1224 n[5][2] = 0.0;
1225
1226 n[6][0] = sqrt2_half;
1227 n[6][1] = -sqrt2_half;
1228 n[6][2] = 0.0;
1229
1230 n[7][0] = -sqrt2_half;
1231 n[7][1] = sqrt2_half;
1232 n[7][2] = 0.0;
1233
1234 //endcaps
1235 n[8][0] = 0.0;
1236 n[8][1] = 0.0;
1237 n[8][2] = 1.0;
1238
1239 n[9][0] = 0.0;
1240 n[9][1] = 0.0;
1241 n[9][2] = -1.0;
1242
1243
1244
1245 if ( abs(X2[2])<zmax ) { //hit in barrel:
1246
1247
1248 if(debug >5)
1249 std::cout<<"---- LineCaloIntersect2 found hit in barrel ----"<<std::endl;
1250
1251 vector < vector <double> > test;
1252 for ( unsigned int staveCounter=0; staveCounter<8; staveCounter++ ) {
1253
1254 double tmp = - (n[staveCounter][0]*n[staveCounter][0] + n[staveCounter][1]*n[staveCounter][1] + n[staveCounter][2]*n[staveCounter][2]) * r_inner
1255 + n[staveCounter][0]*X1[0] + n[staveCounter][1]*X1[1] + n[staveCounter][2]*X1[2];
1256
1257 double tmp2 = n[staveCounter][0]*(X1[0]-X2[0])
1258 + n[staveCounter][1]*(X1[1]-X2[1])
1259 + n[staveCounter][2]*(X1[2]-X2[2]);
1260
1261 double t = -2.0;
1262 if ( tmp2!=0.0 ) {
1263 t = tmp/tmp2;
1264 //for stave[0]: t = (-r_inner+clusterCenterY)/(400*clusterDirectionY)
1265 }
1266
1267 if (t>0.0 && t<=1.0 ) {
1268 //for stave[0]: if (clusterCenterY>r_inner && -r_inner+clusterCenterY<400*clusterDirectionY) ->hit inside ECAL
1269 vector <double> tmpi;
1270 tmpi.push_back(X1[0]+(X2[0]-X1[0])*t); //2*clusterCenterX-r_inner ->shifts cluster center 2*depper into the ECAL stave
1271 tmpi.push_back(X1[1]+(X2[1]-X1[1])*t);
1272 tmpi.push_back(X1[2]+(X2[2]-X1[2])*t);
1273 test.push_back(tmpi);
1274 }
1275 }//all staves
1276
1277 if (test.size()==1 ) {
1278
1279 X[0] = test[0][0];
1280 X[1] = test[0][1];
1281 X[2] = test[0][2];
1282 double rts = sqrt(X[0]*X[0] + X[1]*X[1]);
1283 if ( rts< ((r_inner+2.0)*1.0823922002924) ) { //why +2 ?? why *1.08...??
1284 return; //check wether hit still in ECAL??
1285 } else {
1286 X[0] = 0.0;
1287 X[1] = 0.0;
1288 X[2] = 0.0;
1289 return;//if so, seteverything back
1290 }
1291 } //if test positive for exactly one stave
1292 else {
1293
1294 double rmin = 1e+22;
1295 unsigned int imin = 999;
1296
1297 //loop over all positive staves
1298 for ( unsigned int i=0; i<test.size(); i++ ) {
1299 double Xa[3];
1300 Xa[0] = test[i][0];
1301 Xa[1] = test[i][1];
1302 Xa[2] = test[i][2];
1303
1304 double tmpr = sqrt(Xa[0]*Xa[0] + Xa[1]*Xa[1]);
1305 if( tmpr<rmin ) {
1306 rmin = tmpr; //search smallest cluster distance to IP
1307 imin = i;
1308 }
1309 }
1310
1311 if ( imin < test.size() ) {
1312 X[0] = test[imin][0]; //return shifted cluster position of stave with cluster closest to ECAL surface
1313 X[1] = test[imin][1];
1314 X[2] = test[imin][2];
1315 return;
1316 } else { //did not find stave smallest cluster distance
1317 X[0] = 0.0;
1318 X[1] = 0.0;
1319 X[2] = 0.0;
1320 return;
1321 }
1322 }// testsize != 1
1323
1324 } else { // endcap start
1325 if(debug >5)
1326 std::cout<<"---- LineCaloIntersect2 found hit in endcap ----"<<std::endl;
1327
1328 double n[8][3];
1329 if ( X2[2]>0.0 ) {
1330 n[0][0] = 0.0;
1331 n[0][1] = 0.0;
1332 n[0][2] = 1.0;
1333 } else {
1334 n[0][0] = 0.0;
1335 n[0][1] = 0.0;
1336 n[0][2] = -1.0;
1337 }
1338
1339
1340 double tmp = -(n[0][0]*n[0][0] + n[0][1]*n[0][1] + n[0][2]*n[0][2]) * zmax
1341 + n[0][0]*X1[0] + n[0][1]*X1[1] + n[0][2]*X1[2];
1342 double tmp2 = n[0][0]*(X1[0]-X2[0])
1343 + n[0][1]*(X1[1]-X2[1])
1344 + n[0][2]*(X1[2]-X2[2]);
1345 double t = 0.0;
1346
1347 if ( tmp2!=0.0 ) t=tmp/tmp2;
1348
1349 if ( t>0.0 && t<=1.0 ) {
1350 X[0] = X1[0]+(X2[0]-X1[0])*t ;
1351 X[1] = X1[1]+(X2[1]-X1[1])*t ;
1352 X[2] = X1[2]+(X2[2]-X1[2])*t ;
1353 return;
1354 } else {
1355 X[0] = 0.0;
1356 X[1] = 0.0;
1357 X[2] = 0.0;
1358 return;
1359 }
1360 }//endcap end
1361 }
1362 void ClusterInCluster(myTmpcl* cl, vector<myTmpcl*>& clv,vector<myTmpcl*>& clout)
1363 {
1364 sort(cl->hits.begin(),cl->hits.end());//sort hits in first cluster
1365 for (unsigned int i=0; i<clv.size(); i++) {//loop over 2nd cluster
1366 sort(clv[i]->hits.begin(),clv[i]->hits.end());//sort hits in 2nd cluster
1367 if ( includes(cl->hits.begin(),cl->hits.end(),clv[i]->hits.begin(),clv[i]->hits.end()) )
1368 //every element of clusterVector[i] included in cl
1369 clout.push_back(clv[i]);//return included cluster
1370 }
1371 // std::cout<<"ClusterInCluster loop over "<<clv.size() <<" hits of clusterVector included in photonSeed2 -> filled tempCluster"<<std::endl;
1372 }
1373
1374 void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clusterVector, int debug)
1375 {
1376
1377 sort(cl->hits.begin(),cl->hits.end());
1378
1379 for ( unsigned int vectorCounter=0; vectorCounter<clusterVector.size(); vectorCounter++ ) {
1380
1381 sort(clusterVector[vectorCounter]->hits.begin(),clusterVector[vectorCounter]->hits.end());
1382
1383 if ( includes(cl->hits.begin(),cl->hits.end(),clusterVector[vectorCounter]->hits.begin(),clusterVector[vectorCounter]->hits.end()) ) {
1384 cl->daughters.push_back(clusterVector[vectorCounter]);
1385 clusterVector[vectorCounter]->parents.push_back(cl);
1386 }
1387 }
1388 }
1389
1390
1391 //called by FindCores2
1392 //returns distance of included cluster hit to surrounding cluster hit ??
1393 double LinePointDistance2( double* X1, double* X2, double* X0){
1394
1395 double tmp1 = X1[1]*X2[0] - X1[0]*X2[1] - X1[1]*X0[0] + X2[1]*X0[0] + X1[0]*X0[1] - X2[0]*X0[1];
1396 double tmp2 = X1[2]*X2[0] - X1[0]*X2[2] - X1[2]*X0[0] + X2[2]*X0[0] + X1[0]*X0[2] - X2[0]*X0[2];
1397 double tmp3 = X1[2]*X2[1] - X1[1]*X2[2] - X1[2]*X0[1] + X2[2]*X0[1] + X1[1]*X0[2] - X2[1]*X0[2];
1398
1399 double tmp4 = (X1[0]-X2[0])*(X1[0]-X2[0]) + (X1[1]-X2[1])*(X1[1]-X2[1]) + (X1[2]-X2[2])*(X1[2]-X2[2]);
1400
1401 double tmp5 = tmp1*tmp1 + tmp2*tmp2 + tmp3*tmp3;
1402 tmp5 = tmp5/tmp4;
1403 return sqrt(tmp5);
1404 }
1405
1406 //called by findCores2
1407 void ModuleNormal2(double* X1,double& zmax, double* X0)
1408 {
1409 if ( abs(X1[2])<zmax ) { //in barrel
1410
1411 //for hits in first layer, phi = angle between hit and y axis
1412 //the following section determines if the first layer hit is a neighbour of the next stave.
1413 //If so it is assigned to both stave vectors -> very similar to precalculation
1414
1415 double phi = atan2(X1[0],X1[1]);
1416 double phi0 = M_PI_4/2.0; //pi/4/2 = 0.39
1417
1418 int i;
1419 if ( phi>=0. )
1420 i = (int)((phi+phi0) / M_PI_4);
1421 else
1422 i = (int)((phi-phi0) / M_PI_4);
1423
1424 double phi_p = i * M_PI_4;
1425 if ( i!=0 ) {
1426
1427 X0[0] = sin(phi_p);
1428 X0[1] = cos(phi_p);
1429 X0[2] = 0.0;
1430 } else {
1431
1432 X0[0] = 0.0;
1433 X0[1] = 1.0;
1434 X0[2] = 0.0;
1435 }
1436 } else {
1437 if ( X1[2]>0. ) {
1438
1439 X0[0] = 0.0;
1440 X0[1] = 0.0;
1441 X0[2] = 1.0;
1442 } else {
1443
1444 X0[0] = 0.0;
1445 X0[1] = 0.0;
1446 X0[2] = -1.0;
1447 }
1448 }
1449 }
1450
1451 double D_cl_cl(myTmpcl* cl1,myTmpcl* cl2)
1452 {
1453
1454 vector <double> dist;
1455
1456 for ( unsigned int cluster1entry=0; cluster1entry<cl1->hits.size(); cluster1entry++ ) {
1457
1458 double x1 = cl1->hits[cluster1entry]->pointt[0];
1459 double y1 = cl1->hits[cluster1entry]->pointt[1];
1460 double z1 = cl1->hits[cluster1entry]->pointt[2];
1461
1462 for ( unsigned int cluster2entry=0; cluster2entry<cl2->hits.size(); cluster2entry++ ) {
1463
1464 double tmp = ( (x1 - cl2->hits[cluster2entry]->pointt[0]) * (x1 - cl2->hits[cluster2entry]->pointt[0])
1465 + (y1 - cl2->hits[cluster2entry]->pointt[1]) * (y1 - cl2->hits[cluster2entry]->pointt[1])
1466 + (z1 - cl2->hits[cluster2entry]->pointt[2]) * (z1 - cl2->hits[cluster2entry]->pointt[2]) ) ;
1467 dist.push_back(tmp);
1468 }
1469 }
1470 sort(dist.begin(),dist.end());
1471 return sqrt(dist[0]); //return smallest distance between cluster1 & cluster2
1472 }
1473
1474 //calculate distance between 2 clusters called by FindCores
1475 double D_cl_cl2(Tmpcl2* cl1,Tmpcl2* cl2)
1476 {
1477
1478 vector <double> dist;
1479
1480 for ( unsigned int cluster1entry=0; cluster1entry<cl1->hits.size(); cluster1entry++ ) {
1481
1482 double x1 = cl1->hits[cluster1entry]->pointt[0];
1483 double y1 = cl1->hits[cluster1entry]->pointt[1];
1484 double z1 = cl1->hits[cluster1entry]->pointt[2];
1485
1486 for ( unsigned int cluster2entry=0; cluster2entry<cl2->hits.size(); cluster2entry++ ) {
1487
1488 double tmp = ( (x1 - cl2->hits[cluster2entry]->pointt[0]) * (x1 - cl2->hits[cluster2entry]->pointt[0])
1489 + (y1 - cl2->hits[cluster2entry]->pointt[1]) * (y1 - cl2->hits[cluster2entry]->pointt[1])
1490 + (z1 - cl2->hits[cluster2entry]->pointt[2]) * (z1 - cl2->hits[cluster2entry]->pointt[2]) ) ;
1491 dist.push_back(tmp);
1492 }
1493 }
1494 sort(dist.begin(),dist.end());
1495 return sqrt(dist[0]); //return smallest distance between cluster1 & cluster2
1496 }
1497
1498 // called by FindCores
1499 inline double Dot2(double* X1,double* X2)
1500 {
1501 double n1 = sqrt(X1[0]*X1[0] + X1[1]*X1[1] + X1[2]*X1[2]);
1502 double n2 = sqrt(X2[0]*X2[0] + X2[1]*X2[1] + X2[2]*X2[2]);
1503
1504 return (X1[0]*X2[0] + X1[1]*X2[1] + X1[2]*X2[2]) / (n1*n2);
1505 }
1506
1507
1508 void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clusterVector, vector<Tmpcl2*>& clout, int debug)
1509 {
1510 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
1511 //FIXME: The main problem seems to be this sort routine, but why?
1512 //This should not touch the content, only the order of the cluster!
1513 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
1514
1515 //sort(cl->hits.begin(),cl->hits.end());//sort hits in first cluster
1516 for (unsigned int i=0; i<clusterVector.size(); i++) {//loop over 2nd cluster
1517 sort(clusterVector[i]->hits.begin(),clusterVector[i]->hits.end());//sort hits in 2nd cluster
1518 if ( includes(cl->hits.begin(),cl->hits.end(),clusterVector[i]->hits.begin(),clusterVector[i]->hits.end()) )
1519 //every element of clusterVector[i] included in cl
1520 clout.push_back(clusterVector[i]);//return included cluster
1521 }
1522 if ( debug >2)
1523 std::cout<<"ClusterInCluster loop over "<<clusterVector.size() <<" hits of clusterVector included in photonSeed2 -> filled tempCluster"<<std::endl;
1524
1525 }
1526
1527
1528 Photon2::~Photon2(){}
1529
1530 Photon2:: Photon2(double E_in, double* direction, double* beginning)
1531 {
1532 //photon energy
1533 double E_e = E_in;
1534
1535 dir[0] = direction[0];
1536 dir[1] = direction[1];
1537 dir[2] = direction[2];
1538
1539 start[0] = beginning[0];
1540 start[1] = beginning[1];
1541 start[2] = beginning[2];
1542
1543 //my Particle Geometrical database finds hits in ECAL barrel
1544 myPGdb::ZONE zones ;
1545 zones = myPGdb::ECAL1_BAR;
1546 //get ECAL parameters
1547 Z = myPGDB[zones].Zeff;
1548 x0eff = myPGDB[zones].x0eff;
1549 sampling = myPGDB[zones].sampling;
1550 Eceff = myPGDB[zones].Eceff;
1551 eprime = myPGDB[zones].eprime;
1552 Rm = myPGDB[zones].Rmeff;
1553
1554 z1 = 0.0251 + 0.00319 * log(E_e*1000.0);
1555 z2 = 0.1162 - 0.000381 * Z;
1556 k1 = 0.659 - 0.00309 * Z;
1557 k2 = 0.645;
1558 k3 = -2.59;
1559 k4 = 0.3585 + 0.0421 * log(E_e*1000.0);
1560 p1 = 2.632 - 0.00094 * Z;
1561 p2 = 0.401 + 0.00187 * Z;
1562 p3 = 1.313 - 0.0686 * log(E_e*1000.0);
1563 y = E_e * 1000.0 / Eceff;
1564
1565 Fs = x0eff/sampling;
1566
1567 Thom = log(y) - 0.858;
1568 Tsam = Thom - 0.59/Fs - 0.53*(1.0-eprime);
1569 alfahom = 0.21 + (0.492+2.38/Z) * log(y);
1570 alfasam = alfahom - 0.444/Fs;
1571 betasam = 0.5;
1572
1573 }
1574
1575
1576 void Photon2::Prob(CalorimeterHit* ch,double cut,double* out)
1577 {
1578
1579 double X[3];
1580 PointOnLine22(start, dir,ch->getPosition(),X);
1581
1582 double TP[3];
1583 TP[0] = ch->getPosition()[0] - start[0];
1584 TP[1] = ch->getPosition()[1] - start[1];
1585 TP[2] = ch->getPosition()[2] - start[2];
1586
1587 // distance must have direction
1588 if ( Dot2(TP,dir)<0.0 ) {
1589 out[0] = 0.0;
1590 out[1] = 0.0;
1591 return ;
1592 }
1593
1594 double pos_X0 = sqrt( (start[0]-X[0])*(start[0]-X[0])+
1595 (start[1]-X[1])*(start[1]-X[1])+
1596 (start[2]-X[2])*(start[2]-X[2]))/x0eff;
1597 double tau = pos_X0/Tsam;
1598
1599 double pos[3];
1600 pos[0] = ch->getPosition()[0];
1601 pos[1] = ch->getPosition()[1];
1602 pos[2] = ch->getPosition()[2];
1603
1604 double radius = LinePointDistance2(start,dir,pos);
1605
1606 // average radial profiles
1607 double RChom = z1 + z2*tau;
1608 double RThom = k1 * (exp(k3*(tau-k2)) + exp(k4*(tau-k2)));
1609 double phom = p1 * exp((p2-tau)/p3 - exp((p2-tau)/p3));
1610
1611 // average radial profiles sampling
1612 double RCsam = RChom - 0.0203*(1.0-eprime) + 0.0397*exp(-tau)/Fs;
1613 double RTsam = RThom - 0.14*(1.0-eprime) - 0.495*exp(-tau)/Fs;
1614 double psam = phom + (1.0-eprime)*(0.348-0.642*exp(-pow(tau-1.0,2.0))/Fs);
1615
1616 double rb = radius/Rm;
1617
1618 if ( rb<20.0 && pos_X0<35.0 ) {
1619 double tmp = pow(betasam*pos_X0,alfasam-1.0)*betasam*exp(-betasam*pos_X0)/gsl_sf_gamma(alfasam);
1620 tmp = Ee*tmp*(psam*2.0*rb*RCsam*RCsam/pow(rb*rb+RCsam*RCsam,2.0) +
1621 (1-psam)*(2.0*rb*RTsam*RTsam/pow(rb*rb+RTsam*RTsam,2.0)))/rb;
1622 if ( tmp>cut ) {
1623 out[0] = tmp;
1624 out[1] = radius;
1625
1626 // std::cout<<"set probability to "<<tmp<<std::endl;
1627
1628 } else {
1629 out[0] = 0.0;
1630 out[1] = 0.0;
1631
1632 // std::cout<<"probability "<<tmp<<" < "<<cut<<" -> set to 0."<<std::endl;
1633 }
1634 } else {
1635 out[0] = 0.0;
1636 out[1] = 0.0;
1637
1638 // std::cout<<"radius "<<rb<<" < 20 || pos_X0 "<<pos_X0<<" < 35 -> set to 0."<<std::endl;
1639 }
1640
1641 }
1642
1643
1644
1645 void PointOnLine3(const double* X1,const double* X2,const float* X0,double* Xline)
1646 {
1647
1648 double tmp1 = (X1[0]-X0[0]) * (X2[0]-X1[0]) + (X1[1]-X0[1]) * (X2[1]-X1[1]) + (X1[2]-X0[2]) * (X2[2]-X1[2]);
1649 double tmp4 = (X1[0]-X2[0]) * (X1[0]-X2[0]) + (X1[1]-X2[1]) * (X1[1]-X2[1]) + (X1[2]-X2[2]) * (X1[2]-X2[2]);
1650 double t = -tmp1/tmp4;
1651
1652
1653 Xline[0] = X1[0] + (X2[0]-X1[0])*t;
1654 Xline[1] = X1[1] + (X2[1]-X1[1])*t;
1655 Xline[2] = X1[2] + (X2[2]-X1[2])*t;
1656
1657 }
1658 void PointOnLine22(const double* Xstart,const double* dir,const float* X0,double* Xline)
1659 {
1660 double X1[3];
1661 double X2[3];
1662 X1[0] = Xstart[0] - dir[0]*100.0;
1663 X1[1] = Xstart[1] - dir[1]*100.0;
1664 X1[2] = Xstart[2] - dir[2]*100.0;
1665
1666 X2[0] = Xstart[0] + dir[0]*100.0;
1667 X2[1] = Xstart[1] + dir[1]*100.0;
1668 X2[2] = Xstart[2] + dir[2]*100.0;
1669
1670
1671 PointOnLine3(X1,X2,X0,Xline);
1672
1673 }
1674
1675 double giveMeEEstimate2(int level,double Ecore, vector<CoreCalib2> cc)
1676 {
1677 for ( unsigned int i=0; i<cc.size(); i++ ) {
1678 // std::cout<<cc[i].level<<" == "<<level<<" ? && "<<cc[i].Emax<<" > "<<Ecore<<" ?"<<std::endl;
1679 if ( cc[i].level==level && cc[i].Emax>Ecore ) { //FIXME: does find more than one appropriate calibration!
1680 /*
1681 std::cout << "For level " << level
1682 << " core with energy: " << Ecore
1683 << " : estimated photon energy = cc[i].b: " << cc[i].b
1684 << " + cc[i].a: " << cc[i].a
1685 << " * Ecore = " << cc[i].b + cc[i].a*Ecore << std::endl;
1686 */
1687 if ( cc[i].b==0.0 || cc[i].a==0.0 )
1688 return Ecore; //if no valid calibration, return uncalibrated value
1689 else
1690 return cc[i].b + cc[i].a*Ecore;
1691 }
1692 else
1693 return Ecore; //if no valid calibration, return uncalibrated value
1694 }
1695 return 0.0;
1696 }
1697
1698
1699 // !!! has to be adopted for different detector models !!!
1700 // relates measured energy in cluster to estimated true energy of the cluster,
1701 // depending on the cluster level
1702 // numbers = energy cut in the enery per hit spectrum per level, for various energies
1703 // E_true = aa*E_cluster + b; E_min < E_cluster < E_max
1704 // can be obtained with KALIBRATOR.cc
1705
1706 void CreateCalibration_LDC00(vector<CoreCalib2>* cc)
1707 {
1708
1709
1710 double aEnom[9] = {0.5,1.0,2.0,3.0,5.0,8.0,12.0,20.0,40.0};
1711 double aa[9][10] = {{0.283,0.357,0.337,0.398,0.000,0.000,0.000,0.000,0.000,0.000},
1712 {0.703,0.346,0.814,0.785,0.939,0.000,0.000,0.000,0.000,0.000},
1713 {0.505,0.821,1.000,1.112,1.527,1.860,0.000,0.000,0.000,0.000},
1714 {0.809,0.892,1.999,1.554,2.170,2.870,0.000,0.000,0.000,0.000},
1715 {1.075,0.827,1.290,2.350,2.123,3.958,4.814,0.000,0.000,0.000},
1716 {1.140,1.310,1.765,3.652,4.030,5.600,6.670,0.000,0.000,0.000},
1717 {0.760,1.845,3.444,2.850,5.150,6.760,10.72,0.000,0.000,0.000},
1718 {3.770,4.460,2.270,5.930,8.730,11.00,17.40,17.77,0.000,0.000},
1719 {1.890,3.560,5.120,7.320,9.650,11.71,21.74,0.000,0.000,0.000}};
1720
1721 double bb[9][10] = {{0.730,0.595,0.772,0.683,0.000,0.000,0.000,0.000,0.000,0.000},
1722 {0.493,1.079,0.458,0.592,0.447,0.000,0.000,0.000,0.000,0.000},
1723 {0.941,0.813,0.766,0.782,0.559,0.405,0.000,0.000,0.000,0.000},
1724 {0.892,0.910,0.550,0.809,0.593,0.347,0.000,0.000,0.000,0.000},
1725 {0.918,1.015,0.967,0.800,0.954,0.510,0.356,0.000,0.000,0.000},
1726 {0.974,0.986,0.970,0.778,0.784,0.593,0.509,0.000,0.000,0.000},
1727 {1.032,0.971,0.869,0.986,0.835,0.742,0.339,0.000,0.000,0.000},
1728 {0.891,0.888,1.025,0.884,0.782,0.620,0.314,0.361,0.000,0.000},
1729 {1.020,0.992,0.977,0.962,0.945,0.943,0.724,0.000,0.000,0.000}};
1730
1731 double aEmin[9][10] = {{0.250,0.200,0.150,0.060,0.000,0.000,0.000,0.000,0.000,0.000},
1732 {0.600,0.500,0.350,0.200,0.100,0.000,0.000,0.000,0.000,0.000},
1733 {1.500,1.200,1.100,0.900,0.650,0.300,0.000,0.000,0.000,0.000},
1734 {2.400,2.100,1.900,1.600,1.300,0.600,0.000,0.000,0.000,0.000},
1735 {4.000,4.100,3.900,3.200,2.900,2.400,1.100,0.000,0.000,0.000},
1736 {7.000,6.800,6.200,5.500,5.000,4.200,2.600,0.000,0.000,0.000},
1737 {10.50,10.00,9.500,8.600,8.400,8.300,5.200,0.000,0.000,0.000},
1738 {18.50,18.20,17.40,16.00,15.50,13.00,11.50,9.000,0.000,0.000},
1739 {36.00,35.00,34.00,33.50,32.50,31.00,28.00,0.000,0.000,0.000}};
1740
1741 double aEmax[9][10]={{0.600,0.500,0.350,0.300,0.000,0.000,0.000,0.000,0.000,0.000},
1742 {1.200,1.100,1.000,0.900,0.700,0.000,0.000,0.000,0.000,0.000},
1743 {2.500,2.200,1.900,1.800,1.400,1.200,0.000,0.000,0.000,0.000},
1744 {3.300,3.100,3.000,2.700,2.400,2.000,0.000,0.000,0.000,0.000},
1745 {5.600,5.500,5.200,4.800,4.300,3.900,2.900,0.000,0.000,0.000},
1746 {9.000,8.800,8.300,7.600,7.200,6.400,5.200,0.000,0.000,0.000},
1747 {13.50,13.00,12.50,11.70,11.10,11.00,8.300,0.000,0.000,0.000},
1748 {22.20,22.80,21.50,20.50,19.00,16.30,15.50,13.50,0.000,0.000},
1749 {44.00,43.50,42.50,41.50,39.20,37.50,33.50,0.000,0.000,0.000}};
1750
1751 for ( unsigned int levelCounter=0; levelCounter<10; levelCounter++ ) {
1752 for ( unsigned int energyCounter=0; energyCounter<12; energyCounter++ ) {
1753 if ( aEmin[energyCounter][levelCounter]!=0.0 && aEmax[energyCounter][levelCounter]!=0.0 ) {
1754
1755 CoreCalib2 cck;
1756 cck.level = levelCounter;
1757 cck.Enom = aEnom[energyCounter];
1758 cck.a = aa[energyCounter][levelCounter];
1759 cck.b = bb[energyCounter][levelCounter];
1760 cck.Emin = aEmin[energyCounter][levelCounter];
1761 cck.Emax = aEmax[energyCounter][levelCounter];
1762
1763 cc->push_back(cck);
1764 }
1765 }
1766 }
1767
1768 }
1769
1770
1771 void CreateCalibration_LDCPrime02Sc(vector<CoreCalib2>* cc)
1772 {
1773 // fit straight line: E_true = aa*E_reco + bb for each energy and level
1774
1775 double aEnom[13] = {0.5,1.0,2.0,3.0,5.0,8.0,12.0,20.0,40.0,70.0,100.0,150.0,200.0}; //nominal photon energy
1776
1777 double aa[13][10] = {{0.390,0.386,0.875,0.327,0.000,0.000,0.000,0.000,0.000,0.000},
1778 {0.574,0.437,0.365,0.380,0.378,0.000,0.000,0.000,0.000,0.000},
1779 {0.693,0.644,0.735,0.701,0.663,0.386,0.000,0.000,0.000,0.000},
1780 {0.992,1.003,0.703,0.920,0.705,0.666,0.000,0.000,0.000,0.000},
1781 {0.942,0.894,0.857,0.685,0.851,0.719,0.000,0.000,0.000,0.000},
1782 {0.850,0.876,1.036,0.806,0.746,0.532,0.000,0.000,0.000,0.000},
1783 {0.911,0.901,0.851,0.847,0.796,0.693,0.468,0.464,0.000,0.000},
1784 {0.832,1.005,1.006,0.920,0.975,0.865,0.655,0.441,0.472,0.000},
1785 {0.969,0.899,0.854,0.939,0.858,0.762,0.716,0.378,0.417,0.000},
1786 {1.056,1.073,1.046,0.894,0.879,1.127,1.140,1.199,0.308,0.650},
1787 {0.738,1.080,2.044,1.089,1.032,0.862,0.885,0.758,1.294,0.422},
1788 {0.879,0.879,0.879,0.879,0.879,0.879,0.879,0.879,0.000,0.000},
1789 {0.836,0.836,0.836,0.836,0.836,0.836,0.386,0.836,0.000,0.000}};
1790
1791 double bb[13][10] = {{0.412,0.436,0.383,0.471,0.000,0.000,0.000,0.000,0.000,0.000},
1792 {0.690,0.784,0.870,0.886,0.926,0.000,0.000,0.000,0.000,0.000},
1793 {1.095,1.209,1.155,1.281,1.422,1.748,0.000,0.000,0.000,0.000},
1794 {0.894,0.965,1.675,1.454,1.986,2.251,0.000,0.000,0.000,0.000},
1795 {1.466,1.804,2.109,2.900,2.760,3.442,0.000,0.000,0.000,0.000},
1796 {2.658,2.722,2.135,3.841,4.589,5.984,0.000,0.000,0.000,0.000},
1797 {3.158,3.588,4.480,5.106,6.226,7.719,9.918,10.47,0.000,0.000},
1798 {6.274,3.616,4.306,6.645,7.095,9.971,14.18,17.14,17.78,0.000},
1799 {5.938,9.247,11.98,10.65,14.84,20.11,25.15,33.79,34.32,0.000},
1800 {3.345,3.711,7.038,18.97,22.79,13.86,21.53,-6.09,62.73,54.05},
1801 {33.77,3.409,-81.3,8.968,17.70,36.57,43.36,60.02,35.92,85.20},
1802 {29.30,29.30,29.30,29.30,29.30,29.30,29.30,29.30,0.000,0.000},
1803 {46.13,46.13,46.13,46.13,46.13,46.13,46.13,46.13,0.000,0.000}};
1804
1805 double aEmin[13][10] = {{0.250,0.200,0.150,0.060,0.000,0.000,0.000,0.000,0.000,0.000},
1806 {0.600,0.500,0.350,0.200,0.100,0.000,0.000,0.000,0.000,0.000},
1807 {1.500,1.200,1.100,0.900,0.650,0.300,0.000,0.000,0.000,0.000},
1808 {2.400,2.100,1.900,1.600,1.300,0.600,0.000,0.000,0.000,0.000},
1809 {4.000,4.100,3.900,3.200,2.900,2.400,0.000,0.000,0.000,0.000},
1810 {4.279,3.848,3.501,4.618,2.439,1.695,0.000,0.000,0.000,0.000},
1811 {9.189,6.657,6.144,5.337,4.811,3.780,2.749,0.954,0.000,0.000},
1812 {14.29,13.64,12.69,11.56,10.76,8.683,6.213,3.714,3.040,0.000},
1813 {31.50,30.82,28.61,28.23,25.75,22.29,18.14,13.61,9.857,0.000},
1814 {57.10,56.41,51.97,51.54,48.47,32.96,39.48,28.08,25.26,23.19},
1815 {62.30,84.58,71.61,68.98,59.43,62.63,53.12,33.60,35.83,34.79},
1816 {134.9,134.9,134.9,194.9,134.9,134.9,134.9,134.9,0.000,0.000},
1817 {182.7,182.7,182.7,182.7,182.7,165.0,135.0,112.0,0.000,0.000}};
1818
1819 double aEmax[13][10]={{0.600,0.500,0.350,0.300,0.000,0.000,0.000,0.000,0.000,0.000},
1820 {1.200,1.100,1.000,0.900,0.700,0.000,0.000,0.000,0.000,0.000},
1821 {2.500,2.200,1.900,1.800,1.400,1.200,0.000,0.000,0.000,0.000},
1822 {3.300,3.100,3.000,2.700,2.400,2.000,0.000,0.000,0.000,0.000},
1823 {5.600,5.500,5.200,4.800,4.300,3.900,0.000,0.000,0.000,0.000},
1824 {8.464,8.442,8.051,6.774,6.945,6.120,0.000,0.000,0.000,0.000},
1825 {12.79,12.45,12.01,11.35,10.21,9.137,7.234,6.474,0.000,0.000},
1826 {20.40,19.61,19.38,18.38,16.55,15.58,12.83,11.91,9.463,0.000},
1827 {40.94,39.63,39.58,36.69,34.23,32.34,27.23,24.58,22.26,0.000},
1828 {71.53,70.61,73.92,66.78,63.47,71.53,49.70,50.04,42.17,34.28},
1829 {136.1,100.0,106.0,139.5,140.0,92.21,85.62,91.04,72.17,132.5},
1830 {152.6,152.6,152.6,152.6,152.6,152.6,152.6,152.6,0.000,0.000},
1831 {204.1,204.1,204.1,204.1,204.1,204.1,204.1,204.1,0.000,0.000}};
1832
1833 for ( unsigned int levelCounter=0; levelCounter<10; levelCounter++ ) {
1834 for ( unsigned int energyCounter=0; energyCounter<12; energyCounter++ ) {
1835 if ( aEmin[energyCounter][levelCounter]!=0.0 && aEmax[energyCounter][levelCounter]!=0.0 ) {
1836
1837 CoreCalib2 cck;
1838 cck.level = levelCounter;
1839 cck.Enom = aEnom[energyCounter];
1840 cck.a = aa[energyCounter][levelCounter];
1841 cck.b = bb[energyCounter][levelCounter];
1842 cck.Emin = aEmin[energyCounter][levelCounter];
1843 cck.Emax = aEmax[energyCounter][levelCounter];
1844
1845 cc->push_back(cck);
1846 }
1847 }
1848 }
1849
1850 }//creatCalibration_LDCPrime02Sc
1851
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.