1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | #include "CCAL/DCCALShower_factory.h" |
10 | |
11 | |
12 | static mutex CCAL_MUTEX; |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | DCCALShower_factory::DCCALShower_factory() |
23 | { |
24 | |
25 | |
26 | MIN_CLUSTER_BLOCK_COUNT = 2; |
27 | MAX_HITS_FOR_CLUSTERING = 80; |
28 | MIN_CLUSTER_SEED_ENERGY = 0.035; |
29 | MIN_CLUSTER_ENERGY = 0.05; |
30 | MAX_CLUSTER_ENERGY = 15.9; |
31 | TIME_CUT = 15.0; |
32 | |
33 | SHOWER_DEBUG = 0; |
34 | DO_NONLINEAR_CORRECTION = 1; |
35 | |
36 | CCAL_RADIATION_LENGTH = 0.86; |
37 | CCAL_CRITICAL_ENERGY = 1.1e-3; |
38 | |
39 | LOG_POS_CONST = 4.2; |
40 | |
41 | |
42 | gPARMS->SetDefaultParameter( "CCAL:SHOWER_DEBUG", SHOWER_DEBUG ); |
43 | gPARMS->SetDefaultParameter( "CCAL:MIN_CLUSTER_BLOCK_COUNT", MIN_CLUSTER_BLOCK_COUNT, |
44 | "minimum number of blocks to form a cluster" ); |
45 | gPARMS->SetDefaultParameter( "CCAL:MIN_CLUSTER_SEED_ENERGY", MIN_CLUSTER_SEED_ENERGY, |
46 | "minimum energy for a block to be considered as a seed for a cluster" ); |
47 | gPARMS->SetDefaultParameter( "CCAL:MIN_CLUSTER_ENERGY", MIN_CLUSTER_ENERGY, |
48 | "minimum allowed cluster energy" ); |
49 | gPARMS->SetDefaultParameter( "CCAL:MAX_CLUSTER_ENERGY", MAX_CLUSTER_ENERGY, |
50 | "maximum allowed cluster energy" ); |
51 | gPARMS->SetDefaultParameter( "CCAL:MAX_HITS_FOR_CLUSTERING", MAX_HITS_FOR_CLUSTERING, |
52 | "maximum hits allowed to call clustering algorithm" ); |
53 | gPARMS->SetDefaultParameter( "CCAL:TIME_CUT", TIME_CUT, |
54 | "time cut for associating CCAL hits together into a cluster" ); |
55 | gPARMS->SetDefaultParameter( "CCAL:DO_NONLINEAR_CORRECTION", DO_NONLINEAR_CORRECTION, |
56 | "set this to zero when no nonlinear correction is desired" ); |
57 | gPARMS->SetDefaultParameter( "CCAL:CCAL_RADIATION_LENGTH", CCAL_RADIATION_LENGTH ); |
58 | gPARMS->SetDefaultParameter( "CCAL:CCAL_CRITICAL_ENERGY", CCAL_CRITICAL_ENERGY ); |
59 | gPARMS->SetDefaultParameter( "CCAL:LOG_POS_CONST", LOG_POS_CONST ); |
60 | |
61 | VERBOSE = 0; |
62 | gPARMS->SetDefaultParameter("DFCALShower:VERBOSE", VERBOSE, "Verbosity level for DFCALShower objects and factories"); |
63 | } |
64 | |
65 | |
66 | |
67 | |
68 | |
69 | |
70 | |
71 | |
72 | |
73 | |
74 | jerror_t DCCALShower_factory::brun(JEventLoop *locEventLoop, int32_t runnumber) |
75 | { |
76 | |
77 | static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
78 | static set<int> runs_announced; |
79 | pthread_mutex_lock(&print_mutex); |
80 | bool print_messages = false; |
81 | if(runs_announced.find(runnumber) == runs_announced.end()){ |
82 | print_messages = true; |
83 | runs_announced.insert(runnumber); |
84 | } |
85 | pthread_mutex_unlock(&print_mutex); |
86 | |
87 | |
88 | DApplication *dapp = dynamic_cast<DApplication*>(eventLoop->GetJApplication()); |
89 | const DGeometry *geom = dapp->GetDGeometry(runnumber); |
90 | |
91 | if (geom) { |
92 | geom->GetTargetZ(m_zTarget); |
93 | geom->GetCCALPosition(m_CCALdX,m_CCALdY,m_CCALfront); |
94 | } |
95 | else{ |
96 | jerr << "No geometry accessible." << endl; |
97 | return RESOURCE_UNAVAILABLE; |
98 | } |
99 | |
100 | |
101 | |
102 | |
103 | |
104 | std::unique_lock<std::mutex> lck(CCAL_MUTEX); |
105 | |
106 | string ccal_profile_file; |
107 | gPARMS->SetDefaultParameter("CCAL_PROFILE_FILE", ccal_profile_file, |
108 | "CCAL profile data file name"); |
109 | |
110 | |
111 | |
112 | map< string,string > profile_file_name; |
113 | JCalibration *jcalib = dapp->GetJCalibration(runnumber); |
114 | |
115 | if( jcalib->GetCalib("/CCAL/profile_data/profile_data_map", profile_file_name) ) |
116 | { |
117 | jout << "Can't find requested /CCAL/profile_data/profile_data_map in CCDB for this run!" |
118 | << endl; |
119 | |
120 | } else if( profile_file_name.find("map_name") != profile_file_name.end() |
121 | && profile_file_name["map_name"] != "None" ) |
122 | { |
123 | JResourceManager *jresman = dapp->GetJResourceManager(runnumber); |
124 | ccal_profile_file = jresman->GetResource(profile_file_name["map_name"]); |
125 | } |
126 | |
127 | if(print_messages) |
128 | jout<<"Reading CCAL profile data from "<<ccal_profile_file<<" ..."<<endl; |
129 | |
130 | |
131 | if(ccal_profile_file.empty()) |
132 | { |
133 | if(print_messages) |
134 | jerr << "Empty file..." << endl; |
135 | return RESOURCE_UNAVAILABLE; |
136 | } |
137 | |
138 | ifstream ccal_profile(ccal_profile_file.c_str()); |
139 | for(int i=0; i<=500; i++) { |
140 | for(int j=0; j<=i; j++) { |
141 | int id1, id2; |
142 | double fcell_hyc, fd2c; |
143 | |
144 | ccal_profile >> id1 >> id2 >> fcell_hyc >> fd2c; |
145 | |
146 | acell[id1][id2] = fcell_hyc; |
147 | acell[id2][id1] = fcell_hyc; |
148 | ad2c[id1][id2] = fd2c; |
149 | ad2c[id2][id1] = fd2c; |
150 | } |
151 | } |
152 | ccal_profile.close(); |
153 | |
154 | lck.unlock(); |
155 | |
156 | |
157 | |
158 | |
159 | |
160 | |
161 | |
162 | for(int icol = 0; icol < MCOL12; ++icol) { |
163 | for(int irow = 0; irow < MROW12; ++irow) { |
164 | if(icol>=5 && icol<=6 && irow>=5 && irow<=6) { stat_ch[irow][icol] = -1; } |
165 | else { stat_ch[irow][icol] = 0; } |
166 | } |
167 | } |
168 | |
169 | |
170 | |
171 | |
172 | |
173 | |
174 | |
175 | vector< vector<double> > timewalk_params; |
176 | if( eventLoop->GetCalib("/CCAL/shower_timewalk_correction",timewalk_params) ) |
177 | jout << "Error loading /CCAL/shower_timewalk_correction !" << endl; |
178 | else { |
179 | if( (int)timewalk_params.size() != 2 ) { |
180 | cout << "DCCALShower_factory: Wrong number of entries to timewalk correction table (should be 144)." << endl; |
181 | for( int ii = 0; ii < 2; ++ii ) { |
182 | timewalk_p0.push_back(0.0); |
183 | timewalk_p1.push_back(0.0); |
184 | timewalk_p2.push_back(0.0); |
185 | timewalk_p3.push_back(0.0); |
186 | } |
187 | } else { |
188 | for( vector< vector<double> >::const_iterator iter = timewalk_params.begin(); |
189 | iter != timewalk_params.end(); ++iter ) { |
190 | if( iter->size() != 4 ) { |
191 | cout << "DCCALShower_factory: Wrong number of values in timewalk correction table (should be 4)" << endl; |
192 | continue; |
193 | } |
194 | |
195 | timewalk_p0.push_back( (*iter)[0] ); |
196 | timewalk_p1.push_back( (*iter)[1] ); |
197 | timewalk_p2.push_back( (*iter)[2] ); |
198 | timewalk_p3.push_back( (*iter)[3] ); |
199 | |
200 | } |
201 | } |
202 | } |
203 | |
204 | |
205 | |
206 | |
207 | |
208 | |
209 | |
210 | vector< vector<double> > nonlin_params; |
211 | if( eventLoop->GetCalib("/CCAL/nonlinear_energy_correction",nonlin_params) ) |
212 | jout << "Error loading /CCAL/nonlinear_energy_correction !" << endl; |
213 | else { |
214 | if( (int)nonlin_params.size() != CCAL_CHANS ) { |
215 | cout << "DCCALShower_factory: Wrong number of entries to nonlinear energy correction table (should be 144)." << endl; |
216 | for( int ii = 0; ii < CCAL_CHANS; ++ii ) { |
217 | Nonlin_p0.push_back(0.0); |
218 | Nonlin_p1.push_back(0.0); |
219 | Nonlin_p2.push_back(0.0); |
220 | Nonlin_p3.push_back(0.0); |
221 | } |
222 | } else { |
223 | for( vector< vector<double> >::const_iterator iter = nonlin_params.begin(); |
224 | iter != nonlin_params.end(); ++iter ) { |
225 | if( iter->size() != 4 ) { |
226 | cout << "DCCALShower_factory: Wrong number of values in nonlinear energy correction table (should be 4)" << endl; |
227 | continue; |
228 | } |
229 | |
230 | Nonlin_p0.push_back( (*iter)[0] ); |
231 | Nonlin_p1.push_back( (*iter)[1] ); |
232 | Nonlin_p2.push_back( (*iter)[2] ); |
233 | Nonlin_p3.push_back( (*iter)[3] ); |
234 | |
235 | } |
236 | } |
237 | } |
238 | |
239 | if( SHOWER_DEBUG ) { |
240 | cout << "\n\nNONLIN_P0 NONLIN_P1 NONLIN_P2 NONLIN_P3" << endl; |
241 | for(int ii = 0; ii < (int)Nonlin_p0.size(); ii++) { |
242 | cout << Nonlin_p0[ii] << " " << Nonlin_p1[ii] << " " << Nonlin_p2[ii] << " " << |
243 | Nonlin_p3[ii] << endl; |
244 | } |
245 | cout << "\n\n"; |
246 | } |
247 | |
248 | |
249 | |
250 | |
251 | return NOERROR; |
252 | } |
253 | |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | jerror_t DCCALShower_factory::evnt(JEventLoop *locEventLoop, uint64_t eventnumber) |
264 | { |
265 | |
266 | |
267 | vector< const DCCALGeometry* > ccalGeomVect; |
268 | locEventLoop->Get( ccalGeomVect ); |
269 | if (ccalGeomVect.size() < 1) |
270 | return OBJECT_NOT_AVAILABLE; |
271 | const DCCALGeometry& ccalGeom = *(ccalGeomVect[0]); |
272 | |
273 | |
274 | |
275 | |
276 | |
277 | |
278 | |
279 | vector< const DCCALHit* > ccalhits; |
280 | locEventLoop->Get( ccalhits ); |
281 | |
282 | int n_hits = static_cast<int>( ccalhits.size() ); |
283 | if( n_hits < 1 || n_hits > MAX_HITS_FOR_CLUSTERING ) return NOERROR; |
284 | |
285 | |
286 | vector< vector< const DCCALHit* > > hitPatterns; |
287 | getHitPatterns( ccalhits, hitPatterns ); |
288 | |
289 | int n_patterns = static_cast<int>( hitPatterns.size() ); |
290 | if( n_patterns < 1 ) return NOERROR; |
291 | |
292 | vector< ccalcluster_t > ccalClusters; |
293 | vector< cluster_t > clusterStorage; |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | for( int ipat = 0; ipat < n_patterns; ipat++ ) |
302 | { |
303 | |
304 | |
305 | |
306 | |
307 | vector< const DCCALHit* > rawHitPattern = hitPatterns[ipat]; |
308 | vector< const DCCALHit* > locHitPattern; |
309 | cleanHitPattern( rawHitPattern, locHitPattern ); |
310 | |
311 | int n_hits = static_cast<int>( locHitPattern.size() ); |
312 | |
313 | vector< int > ia; |
314 | vector< int > id; |
315 | |
316 | for( int ih = 0; ih < n_hits; ih++ ) { |
317 | const DCCALHit *ccalhit = locHitPattern[ih]; |
318 | int row = 12 - ccalhit->row; |
319 | int col = 12 - ccalhit->column; |
320 | int ie = static_cast<int>( ccalhit->E*10. + 0.5 ); |
321 | if( ie > 0 ) { |
322 | int address = 100*col + row; |
323 | ia.push_back( address ); |
324 | id.push_back( ie ); |
325 | } |
326 | } |
327 | |
328 | |
329 | |
330 | |
331 | vector< gamma_t > gammas; gammas.clear(); |
332 | main_island( ia, id, gammas ); |
333 | |
334 | |
335 | |
336 | |
337 | int init_clusters = static_cast<int>( gammas.size() ); |
338 | if( !init_clusters ) continue; |
339 | |
340 | processShowers( gammas, ccalGeom, locHitPattern, eventnumber, |
341 | ccalClusters, clusterStorage ); |
342 | |
343 | |
344 | } |
345 | |
346 | |
347 | |
348 | |
349 | |
350 | |
351 | int n_clusters = static_cast<int>( ccalClusters.size() ); |
352 | |
353 | for( int k = 0; k < n_clusters; k++ ) { |
354 | |
355 | DCCALShower *shower = new DCCALShower; |
356 | |
357 | shower->E = ccalClusters[k].E; |
358 | shower->Esum = ccalClusters[k].Esum; |
359 | |
360 | shower->x = ccalClusters[k].x+m_CCALdX; |
361 | shower->y = ccalClusters[k].y+m_CCALdY; |
362 | shower->x1 = ccalClusters[k].x1+m_CCALdX; |
363 | shower->y1 = ccalClusters[k].y1+m_CCALdY; |
364 | shower->z = m_CCALfront; |
365 | |
366 | shower->chi2 = ccalClusters[k].chi2; |
367 | shower->sigma_E = ccalClusters[k].sigma_E; |
368 | shower->Emax = ccalClusters[k].emax; |
369 | shower->time = ccalClusters[k].time; |
370 | |
371 | shower->dime = ccalClusters[k].nhits; |
372 | shower->idmax = ccalClusters[k].idmax; |
373 | shower->id = ccalClusters[k].id; |
374 | |
375 | shower->type = ccalClusters[k].type; |
376 | |
377 | int showTypeVal = (ccalClusters[k].type)/10; |
378 | shower->ClusterType = static_cast<ClusterType_t>( showTypeVal/2 ); |
379 | shower->PeakType = static_cast<PeakType_t>( showTypeVal%2 ); |
380 | |
381 | float xerr = 0.1016; |
382 | float yerr = 0.1016; |
383 | float zerr = 2.5; |
384 | float terr = 0.2; |
385 | float eerr = 0.5; |
386 | |
387 | shower->ExyztCovariance.ResizeTo(5,5); |
388 | for (int col = 0; col < 5; col ++) { |
389 | for (int row = 0; row < 5; row ++) { |
390 | if (col != row) shower->ExyztCovariance[col][row] = 0; |
391 | } |
392 | } |
393 | shower->ExyztCovariance[0][0] = pow(xerr, 2); |
394 | shower->ExyztCovariance[1][1] = pow(yerr, 2); |
395 | shower->ExyztCovariance[2][2] = pow(zerr, 2); |
396 | shower->ExyztCovariance[3][3] = pow(terr, 2); |
397 | shower->ExyztCovariance[4][4] = pow(eerr, 2); |
398 | |
399 | if (VERBOSE>2) {printf("(E,x,y,z,t) "); shower->ExyztCovariance.Print(); } |
400 | |
401 | shower->hitsInCluster.clear(); |
402 | for( int icell = 0; icell < ccalClusters[k].nhits; icell++ ) { |
403 | |
404 | int hitID = clusterStorage[k].id[icell]; |
405 | float hitX = ccalGeom.positionOnFace( hitID ).X(); |
406 | float hitY = ccalGeom.positionOnFace( hitID ).Y(); |
407 | int hitROW = ccalGeom.row( hitY ); |
408 | int hitCOL = ccalGeom.column( hitX ); |
409 | |
410 | DCCALHit clusHit; |
411 | clusHit.row = hitROW; |
412 | clusHit.column = hitCOL; |
413 | clusHit.x = hitX; |
414 | clusHit.y = hitY; |
415 | clusHit.E = static_cast<float>( 1000.*clusterStorage[k].E[icell] ); |
416 | clusHit.t = static_cast<float>( clusterStorage[k].t[icell] ); |
417 | |
418 | shower->hitsInCluster.push_back( clusHit ); |
419 | |
420 | } |
421 | |
422 | _data.push_back( shower ); |
423 | } |
424 | |
425 | |
426 | return NOERROR; |
427 | |
428 | } |
429 | |
430 | |
431 | |
432 | |
433 | |
434 | |
435 | |
436 | |
437 | |
438 | |
439 | void DCCALShower_factory::getHitPatterns( vector< const DCCALHit* > hitarray, |
440 | vector< vector< const DCCALHit* > > &hitPatterns ) |
441 | { |
442 | |
443 | |
444 | |
445 | |
446 | |
447 | |
448 | |
449 | |
450 | |
451 | |
452 | |
453 | |
454 | |
455 | |
456 | |
457 | |
458 | int n_hits = static_cast<int>( hitarray.size() ); |
459 | |
460 | if( n_hits < 1 ) return; |
461 | if( n_hits < 2 ) { |
462 | vector< const DCCALHit* > hitVec; |
463 | hitVec.push_back( hitarray[0] ); |
464 | hitPatterns.push_back( hitVec ); |
465 | return; |
466 | } |
467 | |
468 | |
469 | vector< const DCCALHit* > clonedHitArray = hitarray; |
470 | |
471 | while( clonedHitArray.size() ) |
472 | { |
473 | |
474 | vector< const DCCALHit* > locHitVec; |
475 | |
476 | float maxE = -1.; |
477 | float maxT = 1.e6; |
478 | |
479 | for( unsigned int ih = 0; ih < clonedHitArray.size(); ih++ ) { |
480 | float trialE = clonedHitArray[ih]->E; |
481 | if( trialE > maxE ) { maxE = trialE; maxT = clonedHitArray[ih]->t; } |
482 | } |
483 | |
484 | if( maxE < 0. ) break; |
485 | |
486 | |
487 | sortByTime( clonedHitArray, maxT ); |
488 | |
489 | |
490 | n_hits = static_cast<int>( clonedHitArray.size() ); |
491 | |
492 | int n_good_hits = 0; |
493 | for( int ih = 0; ih < n_hits; ih++ ) { |
494 | |
495 | const DCCALHit *locHit = clonedHitArray[ih]; |
496 | float timeDiff = fabs( locHit->t - maxT ); |
497 | |
498 | if( timeDiff < TIME_CUT ) { |
499 | locHitVec.push_back( locHit ); |
500 | n_good_hits++; |
501 | } else { break; } |
502 | |
503 | } |
504 | |
505 | if( locHitVec.size() ) hitPatterns.push_back( locHitVec ); |
506 | |
507 | clonedHitArray.erase( clonedHitArray.begin(), clonedHitArray.begin() + n_good_hits ); |
508 | |
509 | } |
510 | |
511 | |
512 | return; |
513 | |
514 | } |
515 | |
516 | |
517 | |
518 | |
519 | |
520 | |
521 | |
522 | |
523 | |
524 | |
525 | void DCCALShower_factory::sortByTime( vector< const DCCALHit* > &hitarray, float hitTime ) |
526 | { |
527 | |
528 | int nhits = static_cast<int>( hitarray.size() ); |
529 | |
530 | if( nhits < 2 ) return; |
531 | |
532 | for( int ih = 1; ih < nhits; ih++ ) |
533 | { |
534 | |
535 | float timeDiff = fabs( hitarray[ih]->t - hitTime ); |
536 | float lastTimeDiff = fabs( hitarray[ih-1]->t - hitTime ); |
537 | |
538 | if( timeDiff <= lastTimeDiff ) |
539 | { |
540 | |
541 | const DCCALHit *Hit = hitarray[ih]; |
542 | |
543 | for( int ii = ih-1; ii >= -1; ii-- ) { |
544 | |
545 | if( ii >= 0 ) { |
546 | |
547 | const DCCALHit *locHit = hitarray[ii]; |
548 | float locTimeDiff = fabs( locHit->t - hitTime ); |
549 | |
550 | if( timeDiff < locTimeDiff ) { |
551 | hitarray[ii+1] = locHit; |
552 | } else { |
553 | hitarray[ii+1] = Hit; |
554 | break; |
555 | } |
556 | } else { |
557 | hitarray[0] = Hit; |
558 | } |
559 | } |
560 | |
561 | } |
562 | |
563 | } |
564 | |
565 | return; |
566 | } |
567 | |
568 | |
569 | |
570 | |
571 | |
572 | |
573 | |
574 | |
575 | |
576 | |
577 | void DCCALShower_factory::cleanHitPattern( vector< const DCCALHit* > hitarray, |
578 | vector< const DCCALHit* > &hitarrayClean ) |
579 | { |
580 | |
581 | for( vector< const DCCALHit* >::const_iterator iHit = hitarray.begin(); |
582 | iHit != hitarray.end(); ++iHit ) { |
583 | |
584 | int id12 = ((*iHit)->row)*12 + (*iHit)->column; |
585 | int findVal = -1; |
586 | for( vector<const DCCALHit*>::size_type ii = 0; ii != hitarrayClean.size(); ++ii ) { |
587 | int id = 12*(hitarrayClean[ii]->row) + hitarrayClean[ii]->column; |
588 | if( id == id12 ) { findVal = (int)ii; break; } |
589 | } |
590 | if( findVal >= 0 ) { |
591 | if( (*iHit)->E > hitarrayClean[findVal]->E ) { |
592 | hitarrayClean.erase( hitarrayClean.begin()+findVal ); |
593 | hitarrayClean.push_back( (*iHit) ); |
594 | } |
595 | } else { |
596 | hitarrayClean.push_back( (*iHit) ); |
597 | } |
598 | |
599 | } |
600 | |
601 | return; |
602 | |
603 | } |
604 | |
605 | |
606 | |
607 | |
608 | |
609 | |
610 | |
611 | |
612 | |
613 | |
614 | void DCCALShower_factory::processShowers( vector< gamma_t > gammas, DCCALGeometry ccalGeom, |
615 | vector< const DCCALHit* > locHitPattern, int eventnumber, |
616 | vector< ccalcluster_t > &ccalClusters, vector< cluster_t > &clusterStorage ) |
617 | { |
618 | |
619 | |
620 | |
621 | |
622 | |
623 | int n_clusters = 0; |
624 | int n_hits = static_cast<int>( locHitPattern.size() ); |
625 | |
626 | int init_clusters = static_cast<int>( gammas.size() ); |
627 | for( int k = 0; k < init_clusters; k++ ) { |
628 | |
629 | ccalcluster_t locCluster; |
630 | cluster_t locClusterStorage; |
631 | |
632 | int type = gammas[k].type; |
633 | int dime = gammas[k].dime; |
634 | int id = gammas[k].id; |
635 | double chi2 = gammas[k].chi2; |
636 | double e = gammas[k].energy; |
637 | double x = gammas[k].x; |
638 | double y = gammas[k].y; |
639 | double xc = gammas[k].xc; |
640 | double yc = gammas[k].yc; |
641 | |
642 | |
643 | |
644 | |
645 | if( dime < MIN_CLUSTER_BLOCK_COUNT ) { continue; } |
646 | if( e < MIN_CLUSTER_ENERGY || e > MAX_CLUSTER_ENERGY ) { continue; } |
647 | |
648 | n_clusters++; |
649 | |
650 | |
651 | |
652 | |
653 | double ecellmax = -1; int idmax = -1; |
654 | double e1 = 0.0; |
655 | for( int j = 0; j < (dime>MAX_CC60 ? MAX_CC60 : dime); j++ ) { |
656 | |
657 | |
658 | double ecell = 1.e-4*static_cast<double>( gammas[k].icl_en[j] ); |
659 | |
660 | int ccal_id = gammas[k].icl_in[j]; |
661 | int kx = 12 - (ccal_id/100), ky = 12 - (ccal_id%100); |
662 | |
663 | ccal_id = ky*12 + kx; |
664 | |
665 | e1 += ecell; |
666 | if( ecell > ecellmax ) { |
667 | ecellmax = ecell; |
668 | idmax = ccal_id; |
669 | } |
670 | |
671 | } |
672 | |
673 | double xmax = ccalGeom.positionOnFace(idmax).X(); |
674 | double ymax = ccalGeom.positionOnFace(idmax).Y(); |
675 | |
676 | |
677 | |
678 | |
679 | double sW = 0.0; |
680 | double xpos = 0.0; |
681 | double ypos = 0.0; |
682 | double W; |
683 | |
684 | for( int j = 0; j < (dime>MAX_CC60 ? MAX_CC60 : dime); j++ ) { |
685 | |
686 | int ccal_id = gammas[k].icl_in[j]; |
687 | int kx = 12 - (ccal_id/100), ky = 12 - (ccal_id%100); |
688 | ccal_id = ky*12 + kx; |
689 | |
690 | double ecell = 1.e-4*static_cast<double>( gammas[k].icl_en[j] ); |
691 | double xcell = ccalGeom.positionOnFace( ccal_id ).X(); |
692 | double ycell = ccalGeom.positionOnFace( ccal_id ).Y(); |
693 | |
694 | if(id%10 == 1 || id%10 == 2) { |
695 | xcell += xc; |
696 | ycell += yc; |
697 | } |
698 | |
699 | double hittime = 0.; |
700 | for( int ihit = 0; ihit < n_hits; ihit++ ) { |
701 | int trialid = 12*( locHitPattern[ihit]->row) + locHitPattern[ihit]->column; |
702 | if( trialid == ccal_id ) { |
703 | hittime = locHitPattern[ihit]->t; |
704 | break; |
705 | } |
706 | } |
707 | |
708 | locClusterStorage.id[j] = ccal_id; |
709 | locClusterStorage.E[j] = ecell; |
710 | locClusterStorage.x[j] = xcell; |
711 | locClusterStorage.y[j] = ycell; |
712 | locClusterStorage.t[j] = hittime; |
713 | |
714 | |
715 | |
716 | |
717 | if( ecell > 0.009 && fabs(xcell-xmax) < 6. && fabs(ycell-ymax) < 6.) { |
718 | W = LOG_POS_CONST + log(ecell/e); |
719 | if( W > 0 ) { |
720 | sW += W; |
721 | xpos += xcell*W; |
722 | ypos += ycell*W; |
723 | } |
724 | } |
725 | |
726 | } |
727 | |
728 | for( int j = dime; j < MAX_CC60; j++ ) |
729 | locClusterStorage.id[j] = -1; |
730 | |
731 | double weightedTime = getEnergyWeightedTime( locClusterStorage, dime ); |
732 | double showerTime = getCorrectedTime( weightedTime, e ); |
733 | |
734 | |
735 | |
736 | |
737 | |
738 | DVector3 vertex(0.0, 0.0, m_zTarget); |
739 | |
740 | double x1, y1; |
741 | double zV = vertex.Z(); |
742 | double z0 = m_CCALfront - zV; |
743 | if(sW) { |
744 | double dz = getShowerDepth( e ); |
745 | double zk = 1. / (1. + dz/z0); |
746 | x1 = zk*xpos/sW; |
747 | y1 = zk*ypos/sW; |
748 | } else { |
749 | printf("WRN bad cluster log. coord at event %i: center id = %i, energy = %f\n", |
750 | eventnumber, idmax, e); |
751 | x1 = 0.0; |
752 | y1 = 0.0; |
753 | } |
754 | |
755 | |
756 | |
757 | |
758 | |
759 | |
760 | |
761 | if( idmax < 0 ) { |
762 | printf("WRN negative idmax recorded at event %i; energy = %f\n", eventnumber, e); |
763 | } |
764 | |
765 | double ecorr = e; |
766 | if( DO_NONLINEAR_CORRECTION ) ecorr = getCorrectedEnergy( e, idmax ); |
767 | |
768 | if( SHOWER_DEBUG ) { |
769 | cout << "\n\nShower energy before correction: " << e << " GeV" << endl; |
770 | cout << "Shower energy after correction: " << ecorr << " GeV\n\n" << endl; |
771 | } |
772 | |
773 | |
774 | |
775 | |
776 | |
777 | |
778 | double se = sqrt( 0.9*0.9*ecorr*ecorr + 2.5*2.5*ecorr + 1.0 ); |
779 | |
780 | se /= 100.; |
781 | |
782 | if( (type%10)==1 ) |
783 | se *= 1.5; |
784 | else if( (type%10)==2 ) |
785 | se *= 1.25; |
786 | |
787 | |
788 | |
789 | |
790 | |
791 | locCluster.type = type; |
792 | locCluster.nhits = dime; |
793 | locCluster.id = id; |
794 | locCluster.idmax = idmax; |
795 | |
796 | locCluster.E = ecorr; |
797 | locCluster.Esum = e1; |
798 | locCluster.x = x; |
799 | locCluster.y = y; |
800 | locCluster.x1 = x1; |
801 | locCluster.y1 = y1; |
802 | locCluster.chi2 = chi2; |
803 | locCluster.time = showerTime; |
804 | locCluster.emax = ecellmax; |
805 | locCluster.sigma_E = se; |
806 | |
807 | clusterStorage.push_back( locClusterStorage ); |
808 | ccalClusters.push_back( locCluster ); |
809 | |
810 | } |
811 | |
812 | |
813 | return; |
814 | } |
815 | |
816 | |
817 | |
818 | |
819 | |
820 | |
821 | |
822 | |
823 | |
824 | |
825 | double DCCALShower_factory::getEnergyWeightedTime( cluster_t clusterStorage, int nHits ) |
826 | { |
827 | |
828 | double weightedtime = 0.; |
829 | double totEn = 0; |
830 | for( int j = 0; j < (nHits > MAX_CC60 ? MAX_CC60 : nHits); j++ ) { |
831 | weightedtime += clusterStorage.t[j]*clusterStorage.E[j]; |
832 | totEn += clusterStorage.E[j]; |
833 | } |
834 | weightedtime /= totEn; |
835 | |
836 | return weightedtime; |
837 | |
838 | } |
839 | |
840 | |
841 | |
842 | |
843 | |
844 | |
845 | |
846 | |
847 | |
848 | |
849 | double DCCALShower_factory::getCorrectedTime( double time, double energy ) |
850 | { |
851 | |
852 | |
853 | int iPar; |
854 | if( energy < 1.0 ) iPar = 0; |
855 | else iPar = 1; |
856 | |
857 | double dt = timewalk_p0[iPar]*exp( timewalk_p1[iPar] + timewalk_p2[iPar]*energy ) + timewalk_p3[iPar]; |
858 | double t_cor = time - dt; |
859 | |
860 | return t_cor; |
861 | |
862 | } |
863 | |
864 | |
865 | |
866 | |
867 | |
868 | |
869 | |
870 | |
871 | |
872 | |
873 | double DCCALShower_factory::getShowerDepth( double energy ) |
874 | { |
875 | |
876 | double z0 = CCAL_RADIATION_LENGTH, e0 = CCAL_CRITICAL_ENERGY; |
877 | double depth = (energy > 0.) ? z0*log(1.+energy/e0) : 0.; |
878 | return depth; |
879 | |
880 | } |
881 | |
882 | |
883 | |
884 | |
885 | |
886 | |
887 | |
888 | |
889 | |
890 | |
891 | double DCCALShower_factory::getCorrectedEnergy( double energy, int id ) |
892 | { |
893 | |
894 | if( id < 0 ) return energy; |
895 | |
896 | if( Nonlin_p1[id] == 0. && Nonlin_p2[id] == 0. && Nonlin_p3[id] == 0.) return energy; |
897 | if( Nonlin_p0[id] == 0. ) return energy; |
898 | |
899 | |
900 | |
901 | double emin = 0., emax = 12.; |
902 | double e0 = (emin+emax)/2.; |
903 | |
904 | double de1 = energy - emin*nonlin_func( emin, id ); |
905 | double de2 = energy - emax*nonlin_func( emax, id ); |
906 | double de = energy - e0*nonlin_func( e0, id ); |
907 | |
908 | while( fabs(emin-emax) > 1.e-5 ) { |
909 | if( de1*de > 0. && de2*de < 0.) { |
910 | emin = e0; |
911 | de1 = energy - emin*nonlin_func( emin, id ); |
912 | } else { |
913 | emax = e0; |
914 | de2 = energy - emax*nonlin_func( emax, id ); |
915 | } |
916 | e0 = (emin+emax)/2.; |
917 | de = energy - e0*nonlin_func( e0, id ); |
918 | } |
919 | |
920 | return e0; |
921 | |
922 | } |
923 | |
924 | |
925 | |
926 | |
927 | |
928 | |
929 | |
930 | |
931 | |
932 | |
933 | double DCCALShower_factory::nonlin_func( double e, int id ) |
934 | { |
935 | |
936 | return pow( (e/Nonlin_p0[id]), Nonlin_p1[id] + Nonlin_p2[id]*e + Nonlin_p3[id]*e*e ); |
937 | |
938 | } |
939 | |
940 | |
941 | |
942 | |
943 | |
944 | |
945 | |
946 | |
947 | |
948 | |
949 | |
950 | |
951 | |
952 | |
953 | |
954 | |
955 | |
956 | |
957 | |
958 | |
959 | |
960 | |
961 | |
962 | |
963 | |
964 | |
965 | |
966 | |
967 | |
968 | |
969 | |
970 | void DCCALShower_factory::main_island( vector<int> &ia, vector<int> &id, vector<gamma_t> &gammas ) |
971 | { |
972 | |
973 | |
974 | |
975 | |
976 | int nhits; |
977 | int ncl; |
978 | vector<int> lencl; |
979 | |
980 | |
981 | nhits = static_cast<int>(ia.size()); |
982 | if( nhits == 0 ) return; |
983 | |
984 | |
985 | |
986 | |
987 | |
988 | |
989 | |
990 | |
991 | |
992 | |
993 | |
994 | |
995 | |
996 | ncl = clus_hyc( nhits, ia, id, lencl ); |
997 | |
998 | |
999 | |
1000 | |
1001 | |
1002 | |
1003 | |
1004 | |
1005 | |
1006 | |
1007 | |
1008 | |
1009 | |
1010 | |
1011 | int nadcgam = 0; |
1012 | |
1013 | if( ncl <= 0 ) return; |
1014 | |
1015 | int ipncl = 0; |
1016 | for( int icl = 0; icl < ncl; icl++ ) { |
1017 | |
1018 | int ecl = 0; |
1019 | for( int ii = 0; ii < lencl[icl]; ii++ ) |
1020 | ecl += id[ipncl+ii]; |
1021 | if( ecl > MIN_ENERGY0. ) { |
1022 | |
1023 | vector< int > icl_a; |
1024 | vector< int > icl_d; |
1025 | |
1026 | icl_a.insert( icl_a.begin(), ia.begin()+ipncl, ia.begin()+ipncl+lencl[icl] ); |
1027 | icl_d.insert( icl_d.begin(), id.begin()+ipncl, id.begin()+ipncl+lencl[icl] ); |
1028 | |
1029 | |
1030 | if( SHOWER_DEBUG ) { |
1031 | |
1032 | cout << "\n\n======================" << endl; |
1033 | cout << "Processing Cluster " << icl << ":" << endl; |
1034 | for( unsigned int ih = 0; ih < icl_a.size(); ih++ ) { |
1035 | cout << icl_a[ih] << " " << icl_d[ih] << endl; |
1036 | } |
1037 | |
1038 | } |
1039 | |
1040 | int before = nadcgam; |
1041 | gams_hyc( lencl[icl], icl_a, icl_d, nadcgam, gammas ); |
1042 | int after = nadcgam; |
1043 | |
1044 | |
1045 | if( SHOWER_DEBUG ) { |
1046 | cout << "Reconstructed " << after-before << " gammas. \n\n" << endl; |
1047 | } |
1048 | |
1049 | |
1050 | if( nadcgam > MADCGAM50 ) { |
1051 | nadcgam = MADCGAM50; |
1052 | break; |
1053 | } |
1054 | } |
1055 | ipncl = ipncl + lencl[icl]; |
1056 | } |
1057 | |
1058 | |
1059 | |
1060 | |
1061 | |
1062 | |
1063 | |
1064 | for( int ig = 0; ig < nadcgam; ig++ ) { |
1065 | gammas[ig].energy /= 10000.; |
1066 | gammas[ig].x = (static_cast<double>(MCOL12+1)-2.*gammas[ig].x)*xsize2.09/2.; |
1067 | gammas[ig].y = (static_cast<double>(MROW12+1)-2.*gammas[ig].y)*ysize2.09/2.; |
1068 | gammas[ig].xc = -gammas[ig].xc*xsize2.09; |
1069 | gammas[ig].yc = -gammas[ig].yc*ysize2.09; |
1070 | } |
1071 | |
1072 | |
1073 | if( nadcgam < 1 ) return; |
1074 | |
1075 | |
1076 | |
1077 | |
1078 | for( int ig = 1; ig < nadcgam; ig++ ) |
1079 | { |
1080 | if( gammas[ig].energy > gammas[ig-1].energy ) { |
1081 | gamma_t ref_gam = gammas[ig]; |
1082 | |
1083 | for( int ii = ig-1; ii >= -1; ii-- ) { |
1084 | if( ii >= 0 ) { |
1085 | if( ref_gam.energy > gammas[ii].energy ) { |
1086 | gammas[ii+1] = gammas[ii]; |
1087 | } else { |
1088 | gammas[ii+1] = ref_gam; |
1089 | break; |
1090 | } |
1091 | } else { |
1092 | gammas[0] = ref_gam; |
1093 | } |
1094 | } |
1095 | } |
1096 | } |
1097 | |
1098 | |
1099 | |
1100 | return; |
1101 | } |
1102 | |
1103 | |
1104 | |
1105 | |
1106 | |
1107 | |
1108 | |
1109 | |
1110 | |
1111 | |
1112 | |
1113 | int DCCALShower_factory::clus_hyc( int nw, vector<int> &ia, vector<int> &id, vector<int> &lencl ) |
1114 | { |
1115 | |
1116 | |
1117 | |
1118 | int maxcl = 200; |
1119 | int ncl; |
1120 | int next, iak; |
1121 | int ib, ie; |
1122 | int ias, iaf; |
1123 | int last, lastcl; |
1124 | int leng; |
1125 | |
1126 | int loclencl[200]; |
1127 | |
1128 | |
1129 | |
1130 | |
1131 | ncl = 0; |
1132 | if( nw < 1 ) return ncl; |
1133 | if( nw < 2 ) { |
1134 | ncl = 1; |
1135 | lencl.push_back(1); |
1136 | return ncl; |
1137 | } |
1138 | |
1139 | order_hyc( nw, ia, id ); |
1140 | |
1141 | ncl = 0; |
1142 | next = 0; |
1143 | |
1144 | for( int k = 1; k < (nw+1); k++ ) { |
1145 | |
1146 | if( k < nw ) iak = ia[k]; |
1147 | if( (iak-ia[k-1] <= 1) && (k < nw) ) continue; |
1148 | |
1149 | ib = next; |
1150 | ie = k-1; |
1151 | next = k; |
1152 | |
1153 | if( ncl >= maxcl ) return ncl; |
1154 | ncl++; |
1155 | |
1156 | loclencl[ncl-1] = next-ib; |
1157 | if(ncl == 1) continue; |
1158 | |
1159 | ias = ia[ib]; |
1160 | iaf = ia[ie]; |
1161 | last = ib-1; |
1162 | lastcl = ncl-1; |
1163 | |
1164 | for( int icl = lastcl; icl > 0; icl-- ) { |
1165 | |
1166 | leng = loclencl[icl-1]; |
1167 | if( (ias-ia[last]) > 100 ) break; |
1168 | |
1169 | for( int ii = last; ii >= last-leng+1; ii-- ) { |
1170 | if( ( ias-ia[ii] ) > 100 ) break; |
1171 | if( ( iaf-ia[ii] ) >= 100 ) { |
1172 | |
1173 | if( (icl < (ncl-1)) && (leng <= 10800) ) { |
1174 | |
1175 | vector< int > iawork; |
1176 | ucopy1( ia, iawork, last-leng+1, leng ); |
1177 | ucopy2( ia, last+1, last+1-leng, ib-last-1 ); |
1178 | ucopy3( iawork, ia, ib-leng, leng ); |
1179 | |
1180 | vector< int > idwork; |
1181 | ucopy1( id, idwork, last-leng+1, leng ); |
1182 | ucopy2( id, last+1, last+1-leng, ib-last-1 ); |
1183 | ucopy3( idwork, id, ib-leng, leng ); |
1184 | |
1185 | for( int jj = icl; jj < ncl-1; jj++ ) { |
1186 | loclencl[jj-1] = loclencl[jj]; |
1187 | } |
1188 | |
1189 | } |
1190 | |
1191 | ib = ib-leng; |
1192 | loclencl[ncl-2] = loclencl[ncl-1]+leng; |
1193 | ncl = ncl-1; |
1194 | break; |
1195 | } |
1196 | } |
1197 | |
1198 | last = last-leng; |
1199 | |
1200 | } |
1201 | } |
1202 | |
1203 | for( int icl = 0; icl < ncl; icl++ ) { |
1204 | lencl.push_back( loclencl[icl] ); |
1205 | } |
1206 | |
1207 | |
1208 | return ncl; |
1209 | } |
1210 | |
1211 | |
1212 | |
1213 | |
1214 | |
1215 | |
1216 | |
1217 | |
1218 | |
1219 | |
1220 | |
1221 | void DCCALShower_factory::order_hyc( int nw, vector<int> &ia, vector<int> &id ) |
1222 | { |
1223 | |
1224 | |
1225 | if( nw < 2 ) return; |
1226 | |
1227 | for( int k = 1; k < nw; k++ ) { |
1228 | |
1229 | if( ia[k] <= ia[k-1] ) { |
1230 | int iak = ia[k]; |
1231 | int idk = id[k]; |
1232 | |
1233 | for( int ii = k-1; ii >= -1; ii = ii-1 ) { |
1234 | |
1235 | if( ii >= 0 ) { |
1236 | if( iak < ia[ii] ) { |
1237 | ia[ii+1] = ia[ii]; |
1238 | id[ii+1] = id[ii]; |
1239 | } else { |
1240 | ia[ii+1] = iak; |
1241 | id[ii+1] = idk; |
1242 | break; |
1243 | } |
1244 | } else { |
1245 | ia[0] = iak; |
1246 | id[0] = idk; |
1247 | } |
1248 | |
1249 | } |
1250 | } |
1251 | } |
1252 | |
1253 | |
1254 | return; |
1255 | } |
1256 | |
1257 | |
1258 | |
1259 | |
1260 | |
1261 | |
1262 | |
1263 | |
1264 | |
1265 | |
1266 | |
1267 | void DCCALShower_factory::gams_hyc( int nadc, vector<int> &ia, vector<int> &id, |
1268 | int &nadcgam, vector<gamma_t> &gammas ) |
1269 | { |
1270 | |
1271 | |
1272 | |
1273 | int ngam0; |
1274 | int niter; |
1275 | int npk; |
1276 | int ipnpk[10]; |
1277 | int igmpk[2][10]; |
1278 | int minpk; |
1279 | int idelta; |
1280 | int itype; |
1281 | int leng; |
1282 | int ixypk, ixpk, iypk; |
1283 | int ic, idc, in; |
1284 | int ixy, ixymax, ixymin; |
1285 | int ix, iy, iyc; |
1286 | int iwk; |
1287 | int iia; |
1288 | int ide; |
1289 | int idecorr; |
1290 | |
1291 | double fw[13]; |
1292 | double chisq; |
1293 | double chisq1; |
1294 | double chisq2; |
1295 | double ratio; |
1296 | double eg; |
1297 | double epk[10]; |
1298 | double xpk[10]; |
1299 | double ypk[10]; |
1300 | double a, dx, dy; |
1301 | double e1, x1, y1; |
1302 | double e2, x2, y2; |
1303 | double fe, fia; |
1304 | |
1305 | int idsum; |
1306 | |
1307 | idelta = 0; |
1308 | chisq1 = 90.; |
1309 | chisq2 = 50.; |
1310 | niter = 6; |
1311 | |
1312 | ngam0 = nadcgam; |
1313 | |
1314 | |
1315 | |
1316 | vector<int> iwrk[13]; |
1317 | vector<int> idp[13]; |
1318 | vector<double> fwrk[13]; |
1319 | |
1320 | |
1321 | |
1322 | for( int ii=0; ii<13; ii++ ) { |
1323 | iwrk[ii].reserve(nadc); |
1324 | fwrk[ii].reserve(nadc); |
1325 | idp[ii].reserve(nadc); |
1326 | for( int ih=0; ih<nadc; ih++ ) { |
1327 | iwrk[ii].push_back(0); |
1328 | fwrk[ii].push_back(0.); |
1329 | idp[ii].push_back(0); |
1330 | } |
1331 | } |
1332 | |
1333 | |
1334 | |
1335 | |
1336 | |
1337 | |
1338 | |
1339 | order_hyc( nadc, ia, id ); |
1340 | |
1341 | idsum = 0; |
1342 | for( int ic = 0; ic < nadc; ic++ ) |
1343 | idsum += id[ic]; |
1344 | |
1345 | if( nadc < 3 ) |
1346 | minpk = 1; |
1347 | else { |
1348 | int trial = 7.*log(1. + 0.0001*static_cast<double>(idsum)) + 0.5; |
1349 | if( trial > 1 ) minpk = trial; |
1350 | else minpk = 1; |
1351 | } |
1352 | minpk *= 100; |
1353 | |
1354 | |
1355 | npk = 0; |
1356 | |
1357 | for( ic = 0; ic < nadc; ic++ ) { |
1358 | idc = id[ic]; |
1359 | if( idc < minpk ) continue; |
1360 | ixy = ia[ic]; |
1361 | ixymax = ixy + 100 + 1; |
1362 | ixymin = ixy - 100 - 1; |
1363 | iyc = ixy - (ixy/100)*100; |
1364 | |
1365 | int peakVal = 1; |
1366 | |
1367 | in = ic+1; |
1368 | while( in < nadc && ia[in] <= ixymax ) { |
1369 | iy = ia[in] - (ia[in]/100)*100; |
1370 | if( abs(iy-iyc) <= 1 && id[in] >= idc ) peakVal = 0; |
1371 | in++; |
1372 | } |
1373 | |
1374 | in = ic-1; |
1375 | while( in >= 0 && ia[in] >= ixymin ) { |
1376 | iy = ia[in] - (ia[in]/100)*100; |
1377 | if( abs(iy-iyc) <= 1 && id[in] > idc ) peakVal = 0; |
1378 | in -= 1; |
1379 | } |
1380 | |
1381 | if( !peakVal ) continue; |
1382 | |
1383 | npk += 1; |
1384 | ipnpk[npk-1] = ic; |
1385 | if( npk == 10 || npk >= 10000/nadc-3 ) break; |
1386 | |
1387 | } |
1388 | |
1389 | if( npk <= 0 ) return; |
1390 | |
1391 | if( SHOWER_DEBUG ) cout << "Found " << npk << " peaks. Now processing..." << endl; |
1392 | |
1393 | |
1394 | |
1395 | |
1396 | if( npk == 1 ) { |
1397 | |
1398 | if( nadcgam >= MADCGAM50-1 ) return; |
1399 | nadcgam = nadcgam+1; |
1400 | chisq = chisq2; |
1401 | |
1402 | ic = ipnpk[0]; |
1403 | ix = ia[ic]/100; |
1404 | iy = ia[ic] - ix*100; |
1405 | |
1406 | itype = peak_type( ix, iy ); |
1407 | |
1408 | e2 = 0.; |
1409 | gamma_hyc( nadc, ia, id, chisq, |
1410 | e1, x1, y1, |
1411 | e2, x2, y2 ); |
1412 | |
1413 | gamma_t gam1; |
1414 | gamma_t gam2; |
1415 | |
1416 | gam1.type = itype; |
1417 | gam1.dime = nadc; |
1418 | gam1.id = 0; |
1419 | |
1420 | gam1.chi2 = chisq; |
1421 | gam1.energy = e1; |
1422 | gam1.x = x1; |
1423 | gam1.y = y1; |
1424 | gam1.xc = 0.; |
1425 | gam1.yc = 0.; |
1426 | |
1427 | if( e2 > 0. && nadcgam <= MADCGAM50-1 ) { |
1428 | nadcgam = nadcgam+1; |
1429 | |
1430 | gam2.type = itype+10; |
1431 | gam2.dime = nadc; |
1432 | gam2.id = 2; |
1433 | |
1434 | gam2.chi2 = chisq; |
1435 | gam2.energy = e2; |
1436 | gam2.x = x2; |
1437 | gam2.y = y2; |
1438 | gam2.xc = 0.5*(x2-x1); |
1439 | gam2.yc = 0.5*(y2-y1); |
1440 | |
1441 | gam1.type = itype+10; |
1442 | gam1.id = 1; |
1443 | gam1.xc = 0.5*(x1-x2); |
1444 | gam1.yc = 0.5*(y1-y2); |
1445 | |
1446 | for( int jj = 0; jj < nadc; jj++ ) { |
1447 | if( jj < MAX_CC60 ) { |
1448 | gam1.icl_in[jj] = ia[jj]; |
1449 | gam2.icl_in[jj] = ia[jj]; |
1450 | gam1.icl_en[jj] = static_cast<int>(static_cast<double>(id[jj])*e1/(e1+e2) + 0.5); |
1451 | gam2.icl_en[jj] = static_cast<int>(static_cast<double>(id[jj])*e2/(e1+e2) + 0.5); |
1452 | } |
1453 | } |
1454 | |
1455 | gammas.push_back( gam1 ); |
1456 | gammas.push_back( gam2 ); |
1457 | |
1458 | } else { |
1459 | for( int jj = 0; jj < nadc; jj++ ) { |
1460 | if( jj < MAX_CC60 ) { |
1461 | gam1.icl_in[jj] = ia[jj]; |
1462 | gam1.icl_en[jj] = id[jj]; |
1463 | } |
1464 | } |
1465 | |
1466 | gammas.push_back( gam1 ); |
1467 | |
1468 | } |
1469 | |
1470 | |
1471 | } else { |
1472 | |
1473 | |
1474 | |
1475 | |
1476 | |
1477 | |
1478 | |
1479 | |
1480 | |
1481 | |
1482 | |
1483 | |
1484 | if( nadcgam >= MADCGAM50-1 ) return; |
1485 | |
1486 | ratio = 1.; |
1487 | for( int iter = 0; iter < niter; iter++ ) { |
1488 | for( int ii = 0; ii < nadc; ii++ ) { |
1489 | iwrk[0].push_back(0); |
1490 | fwrk[0].push_back(0.); |
1491 | } |
1492 | |
1493 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1494 | |
1495 | ic = ipnpk[ipk]; |
1496 | if( iter > 0 ) ratio = fwrk[ipk+1][ic]/fwrk[npk+1][ic]; |
1497 | eg = ratio*static_cast<double>(id[ic]); |
1498 | ixypk = ia[ic]; |
1499 | ixpk = ixypk/100; |
1500 | iypk = ixypk - ixpk*100; |
1501 | epk[ipk] = eg; |
1502 | xpk[ipk] = eg*static_cast<double>(ixpk); |
1503 | ypk[ipk] = eg*static_cast<double>(iypk); |
1504 | |
1505 | if( ic < nadc-1 ) { |
1506 | for( int ii = ic+1; ii < nadc; ii++ ) { |
1507 | ixy = ia[ii]; |
1508 | ix = ixy/100; |
1509 | iy = ixy - ix*100; |
1510 | if( ixy-ixypk > 100+1 ) break; |
1511 | if( abs(iy-iypk) <= 1 ) { |
1512 | if( iter != 0 ) ratio = fwrk[ipk+1][ii]/fwrk[npk+1][ii]; |
1513 | eg = ratio*static_cast<double>(id[ii]); |
1514 | epk[ipk] = epk[ipk] + eg; |
1515 | xpk[ipk] = xpk[ipk] + eg*static_cast<double>(ix); |
1516 | ypk[ipk] = ypk[ipk] + eg*static_cast<double>(iy); |
1517 | } |
1518 | } |
1519 | } |
1520 | |
1521 | if( ic > 0 ) { |
1522 | for( int ii = ic-1; ii >= 0; ii--) { |
1523 | ixy = ia[ii]; |
1524 | ix = ixy/100; |
1525 | iy = ixy - ix*100; |
1526 | if( ixypk-ixy > 100+1 ) break; |
1527 | if( abs(iy-iypk) <= 1 ) { |
1528 | if( iter != 0 ) ratio = fwrk[ipk+1][ii]/fwrk[npk+1][ii]; |
1529 | eg = ratio*static_cast<double>(id[ii]); |
1530 | epk[ipk] = epk[ipk] + eg; |
1531 | xpk[ipk] = xpk[ipk] + eg*static_cast<double>(ix); |
1532 | ypk[ipk] = ypk[ipk] + eg*static_cast<double>(iy); |
1533 | } |
1534 | } |
1535 | } |
1536 | |
1537 | if( epk[ipk] > 0. ) { |
1538 | xpk[ipk] = xpk[ipk]/epk[ipk]; |
1539 | ypk[ipk] = ypk[ipk]/epk[ipk]; |
1540 | } |
1541 | |
1542 | for( int ii = 0; ii < nadc; ii++ ) { |
1543 | ixy = ia[ii]; |
1544 | ix = ixy/100; |
1545 | iy = ixy - ix*100; |
1546 | dx = fabs( static_cast<double>(ix) - xpk[ipk] ); |
1547 | dy = fabs( static_cast<double>(iy) - ypk[ipk] ); |
1548 | |
1549 | a = epk[ipk]*cell_hyc( dx, dy ); |
1550 | |
1551 | fwrk[ipk+1][ii] = a; |
1552 | fwrk[0][ii] = fwrk[0][ii] + fwrk[ipk+1][ii]; |
1553 | iwrk[ipk+1][ii] = static_cast<int>(a + 0.5); |
1554 | iwrk[0][ii] = iwrk[0][ii] + iwrk[ipk+1][ii]; |
1555 | } |
1556 | |
1557 | |
1558 | } |
1559 | |
1560 | for( int ii = 0; ii < nadc; ii++ ) { |
1561 | iwk = iwrk[0][ii]; |
1562 | if( iwk < 1 ) iwk = 1; |
1563 | iwrk[npk+1][ii] = iwk; |
1564 | |
1565 | if( fwrk[0][ii] > 1.e-2 ) |
1566 | fwrk[npk+1][ii] = fwrk[0][ii]; |
1567 | else |
1568 | fwrk[npk+1][ii] = 1.e-2; |
1569 | |
1570 | } |
1571 | } |
1572 | |
1573 | |
1574 | |
1575 | if( SHOWER_DEBUG ) { |
1576 | |
1577 | |
1578 | cout << "\n\n\nAfter 6 iterations: " << endl; |
1579 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1580 | cout << "peak " << ipk+1 << ": " << endl; |
1581 | for( int jj = 0; jj < nadc; jj++ ) { |
1582 | cout << ia[jj] <<"; "<< id[jj] <<"; "<< fwrk[ipk+1][jj] << endl; |
1583 | } |
1584 | } |
1585 | |
1586 | } |
1587 | |
1588 | |
1589 | |
1590 | |
1591 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1592 | |
1593 | vector<int> iwrk_a; |
1594 | vector<int> iwrk_d; |
1595 | |
1596 | leng = 0; |
1597 | for( int ii = 0; ii < nadc; ii++ ) { |
1598 | if( fwrk[0][ii] > 1.e-2 ) { |
1599 | ixy = ia[ii]; |
1600 | fe = static_cast<double>(id[ii])*fwrk[ipk+1][ii]/fwrk[0][ii]; |
1601 | |
1602 | if( fe > idelta ) { |
1603 | leng = leng+1; |
1604 | iwrk_a.push_back( ixy ); |
1605 | iwrk_d.push_back( static_cast<int>(fe+0.5) ); |
1606 | } |
1607 | |
1608 | } |
1609 | } |
1610 | |
1611 | if( nadcgam >= MADCGAM50-1 ) return; |
1612 | |
1613 | igmpk[1][ipk] = 0; |
1614 | if( leng == 0 ) continue; |
1615 | nadcgam = nadcgam + 1; |
1616 | chisq = chisq1; |
1617 | |
1618 | ic = ipnpk[ipk]; |
1619 | ix = ia[ic]/100; |
1620 | iy = ia[ic] - ix*100; |
1621 | |
1622 | itype = peak_type( ix, iy ); |
1623 | |
1624 | e2 = 0.; |
1625 | gamma_hyc( leng, iwrk_a, iwrk_d, chisq, |
1626 | e1, x1, y1, e2, x2, y2 ); |
1627 | |
1628 | gamma_t gam1; |
1629 | gamma_t gam2; |
1630 | |
1631 | gam1.chi2 = chisq; |
1632 | gam1.type = itype+40; |
1633 | gam1.energy = e1; |
1634 | gam1.x = x1; |
1635 | gam1.y = y1; |
1636 | gam1.dime = leng; |
1637 | |
1638 | gam1.id = 90; |
1639 | |
1640 | igmpk[0][ipk] = nadcgam; |
1641 | igmpk[1][ipk] = nadcgam; |
1642 | |
1643 | if( e2 > 0. && nadcgam <= MADCGAM50-1 ) { |
1644 | nadcgam = nadcgam+1; |
1645 | |
1646 | gam2.chi2 = chisq; |
1647 | gam2.type = itype+50; |
1648 | gam2.energy = e2; |
1649 | gam2.x = x2; |
1650 | gam2.y = y2; |
1651 | gam2.xc = 0.5*(x2-x1); |
1652 | gam2.yc = 0.5*(y2-y1); |
1653 | gam1.xc = 0.5*(x1-x2); |
1654 | gam1.yc = 0.5*(y1-y2); |
1655 | |
1656 | gam2.id = 92; |
1657 | gam1.id = 91; |
1658 | |
1659 | gam2.dime = leng; |
1660 | igmpk[1][ipk] = nadcgam; |
1661 | |
1662 | gammas.push_back( gam1 ); |
1663 | gammas.push_back( gam2 ); |
1664 | |
1665 | } else { gammas.push_back( gam1 ); } |
1666 | |
1667 | } |
1668 | |
1669 | |
1670 | |
1671 | |
1672 | |
1673 | |
1674 | |
1675 | |
1676 | for( int ii = 0; ii < nadc; ii++ ) { |
1677 | iwrk[0][ii] = 0; |
1678 | fwrk[0][ii] = 0.; |
1679 | idp[0][ii] = 0; |
1680 | } |
1681 | |
1682 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1683 | for( int ii = 0; ii < nadc; ii++ ) { |
1684 | iwrk[ipk+1][ii] = 0; |
1685 | fwrk[ipk+1][ii] = 0.; |
1686 | idp[ipk+1][ii] = 0; |
1687 | |
1688 | if( igmpk[1][ipk] == 0 ) continue; |
1689 | |
1690 | for( int ig = igmpk[0][ipk]; ig <= igmpk[1][ipk]; ig++ ) { |
1691 | ixy = ia[ii]; |
1692 | ix = ixy/100; |
1693 | iy = ixy-(ix*100); |
1694 | dx = static_cast<double>(ix) - gammas[ig-1].x; |
1695 | dy = static_cast<double>(iy) - gammas[ig-1].y; |
1696 | |
1697 | fia = gammas[ig-1].energy*cell_hyc( dx, dy ); |
1698 | iia = static_cast<int>(fia+0.5); |
1699 | |
1700 | |
1701 | iwrk[ipk+1][ii] += iia; |
1702 | fwrk[ipk+1][ii] += fia; |
1703 | idp[ipk+1][ii] += iia; |
1704 | |
1705 | iwrk[0][ii] += iia; |
1706 | fwrk[0][ii] += fia; |
1707 | |
1708 | } |
1709 | |
1710 | } |
1711 | } |
1712 | |
1713 | |
1714 | |
1715 | |
1716 | |
1717 | |
1718 | for( int ii = 0; ii < nadc; ii++ ) { |
1719 | |
1720 | idp[0][ii] = 0; |
1721 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1722 | idp[0][ii] += idp[ipk+1][ii]; |
1723 | } |
1724 | |
1725 | ide = id[ii] - idp[0][ii]; |
1726 | if( ide == 0 ) continue; |
1727 | if( fwrk[0][ii] == 0. ) continue; |
1728 | |
1729 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1730 | fw[ipk+1] = fwrk[ipk+1][ii]/fwrk[0][ii]; |
1731 | } |
1732 | |
1733 | idecorr = 0; |
1734 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1735 | fia = ide*fw[ipk+1]; |
1736 | if( (fwrk[ipk+1][ii] + fia) > 0. ) { |
1737 | fwrk[ipk+1][ii] += fia; |
1738 | fwrk[0][ii] += fia; |
1739 | } |
1740 | iia = static_cast<int>(fia+0.5); |
1741 | if( (iwrk[ipk+1][ii] + iia) > 0 ) { |
1742 | iwrk[ipk+1][ii] += iia; |
1743 | iwrk[0][ii] += iia; |
1744 | idecorr += iia; |
1745 | } else if( (iwrk[ipk+1][ii] + iia) < 0 ) { |
1746 | |
1747 | } |
1748 | } |
1749 | |
1750 | } |
1751 | |
1752 | |
1753 | |
1754 | |
1755 | |
1756 | nadcgam = ngam0; |
1757 | gammas.erase( gammas.begin()+nadcgam, gammas.end() ); |
1758 | |
1759 | |
1760 | |
1761 | |
1762 | |
1763 | |
1764 | for( int ipk = 0; ipk < npk; ipk++ ) { |
1765 | leng = 0; |
1766 | |
1767 | vector<int> iwrk_a; |
1768 | vector<int> iwrk_d; |
1769 | |
1770 | for( int ii = 0; ii < nadc; ii++ ){ |
1771 | if( iwrk[0][ii] > 0 ) { |
1772 | fe = id[ii]*fwrk[ipk+1][ii]/fwrk[0][ii]; |
1773 | if( fe > idelta ) { |
1774 | leng++; |
1775 | iwrk_a.push_back( ia[ii] ); |
1776 | iwrk_d.push_back( static_cast<int>(fe+0.5) ); |
1777 | } |
1778 | } |
1779 | } |
1780 | |
1781 | if( nadcgam >= MADCGAM50-1 ) return; |
1782 | |
1783 | |
1784 | if( leng == 0 ) continue; |
1785 | |
1786 | nadcgam++; |
1787 | |
1788 | chisq = chisq2; |
1789 | |
1790 | ic = ipnpk[ipk]; |
1791 | ix = ia[ic]/100; |
1792 | iy = ia[ic] - ix*100; |
1793 | |
1794 | itype = peak_type( ix, iy ); |
1795 | |
1796 | e2 = 0.; |
1797 | gamma_hyc( leng, iwrk_a, iwrk_d, chisq, |
1798 | e1, x1, y1, e2, x2, y2 ); |
1799 | |
1800 | gamma_t gam1; |
1801 | gamma_t gam2; |
1802 | |
1803 | gam1.type = itype+20; |
1804 | gam1.dime = leng; |
1805 | gam1.id = 10; |
1806 | |
1807 | gam1.chi2 = chisq; |
1808 | gam1.energy = e1; |
1809 | gam1.x = x1; |
1810 | gam1.y = y1; |
1811 | gam1.xc = 0.; |
1812 | gam1.yc = 0.; |
1813 | |
1814 | if( e2 > 0. && nadcgam <= MADCGAM50-1 ) { |
1815 | nadcgam = nadcgam+1; |
1816 | |
1817 | gam2.type = itype+30; |
1818 | gam2.dime = leng; |
1819 | gam2.id = 12; |
1820 | |
1821 | gam2.chi2 = chisq; |
1822 | gam2.energy = e2; |
1823 | gam2.x = x2; |
1824 | gam2.y = y2; |
1825 | gam2.xc = 0.5*(x2-x1); |
1826 | gam2.yc = 0.5*(y2-y1); |
1827 | |
1828 | gam1.type = itype+30; |
1829 | gam1.id = 11; |
1830 | gam1.xc = 0.5*(x1-x2); |
1831 | gam1.yc = 0.5*(x1-x2); |
1832 | |
1833 | for( int jj = 0; jj < leng; jj++ ) { |
1834 | if( jj < MAX_CC60 ) { |
1835 | gam1.icl_in[jj] = iwrk_a[jj]; |
1836 | gam2.icl_in[jj] = iwrk_a[jj]; |
1837 | gam1.icl_en[jj] = static_cast<int>( |
1838 | static_cast<double>(iwrk_d[jj])*e1/(e1+e2) + 0.5 ); |
1839 | gam2.icl_en[jj] = static_cast<int>( |
1840 | static_cast<double>(iwrk_d[jj])*e2/(e1+e2) + 0.5 ); |
1841 | } |
1842 | } |
1843 | |
1844 | gammas.push_back( gam1 ); |
1845 | gammas.push_back( gam2 ); |
1846 | |
1847 | } else { |
1848 | for( int jj = 0; jj < leng; jj++ ) { |
1849 | if( jj < MAX_CC60 ) { |
1850 | gam1.icl_in[jj] = iwrk_a[jj]; |
1851 | gam1.icl_en[jj] = iwrk_d[jj]; |
1852 | } |
1853 | } |
1854 | |
1855 | gammas.push_back( gam1 ); |
1856 | |
1857 | } |
1858 | |
1859 | |
1860 | } |
1861 | |
1862 | |
1863 | |
1864 | } |
1865 | |
1866 | |
1867 | |
1868 | |
1869 | |
1870 | return; |
1871 | } |
1872 | |
1873 | |
1874 | |
1875 | |
1876 | |
1877 | |
1878 | |
1879 | |
1880 | |
1881 | |
1882 | |
1883 | int DCCALShower_factory::peak_type( int ix, int iy ) |
1884 | { |
1885 | |
1886 | |
1887 | |
1888 | |
1889 | |
1890 | |
1891 | int itype = 0; |
1892 | if( (ix == MCOL12/2-1 || ix == MCOL12/2+2) && (iy >= MROW12/2-1 && iy <= MROW12/2+2) ) itype = 1; |
1893 | if( (iy == MROW12/2-1 || iy == MROW12/2+2) && (ix >= MCOL12/2-1 && ix <= MCOL12/2+2) ) itype = 1; |
1894 | if( ix == 1 || ix == MCOL12 || iy == 1 || iy == MROW12 ) itype = 2; |
1895 | |
1896 | return itype; |
1897 | } |
1898 | |
1899 | |
1900 | |
1901 | |
1902 | |
1903 | |
1904 | |
1905 | |
1906 | |
1907 | |
1908 | |
1909 | void DCCALShower_factory::gamma_hyc( int nadc, vector<int> ia, vector<int> id, double &chisq, |
1910 | double &e1, double &x1, double &y1, |
1911 | double &e2, double &x2, double &y2 ) |
1912 | { |
1913 | |
1914 | |
1915 | int dof; |
1916 | |
1917 | double dxy; |
1918 | double stepmin; |
1919 | double stepx, stepy; |
1920 | double parx, pary; |
1921 | double chimem, chisq0, chi0; |
1922 | double chiold; |
1923 | double chi00; |
1924 | double x0, y0; |
1925 | double ee, xx, yy; |
1926 | double d2, xm2, xm2cut; |
1927 | double chir, chil, chiu, chid; |
1928 | |
1929 | dxy = 0.05; |
1930 | stepmin = 0.002; |
1931 | xm2cut = 1.7; |
1932 | |
1933 | int nzero; |
1934 | vector<int> iaz; |
1935 | |
1936 | |
1937 | |
1938 | e2 = 0.; |
1939 | x2 = 0.; |
1940 | y2 = 0.; |
1941 | |
1942 | |
1943 | fill_zeros( nadc, ia, nzero, iaz ); |
1944 | mom1_pht( nadc, ia, id, nzero, iaz, e1, x1, y1 ); |
1945 | |
1946 | if( nadc <= 0 ) return; |
1947 | |
1948 | chimem = chisq; |
1949 | chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x1, y1, chi0 ); |
1950 | |
1951 | |
1952 | chisq0 = chi0; |
1953 | dof = nzero + nadc - 2; |
1954 | if( dof < 1 ) dof = 1; |
1955 | chisq = chi0/dof; |
1956 | x0 = x1; |
1957 | y0 = y1; |
1958 | |
1959 | |
1960 | |
1961 | |
1962 | int rounds = 0; |
1963 | while( 1 ) { |
1964 | |
1965 | chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0+dxy, y0, chir ); |
1966 | chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0-dxy, y0, chil ); |
1967 | chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0, y0+dxy, chiu ); |
1968 | chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0, y0-dxy, chid ); |
1969 | |
1970 | if( chi0 > chir || chi0 > chil ) { |
1971 | stepx = dxy; |
1972 | if( chir > chil ) stepx = -stepx; |
1973 | } else { |
1974 | stepx = 0.; |
1975 | parx = chir + chil - 2.*chi0; |
1976 | if( parx > 0. ) stepx = -dxy*(chir-chil)/(2.*parx); |
1977 | } |
1978 | if( chi0 > chiu || chi0 > chid ) { |
1979 | stepy = dxy; |
1980 | if( chiu > chid ) stepy = -stepy; |
1981 | } else { |
1982 | stepy = 0.; |
1983 | pary = chiu + chid - 2.*chi0; |
1984 | if( pary > 0. ) stepy = -dxy*(chiu-chid)/(2.*pary); |
1985 | } |
1986 | |
1987 | |
1988 | |
1989 | |
1990 | if( fabs(stepx) < stepmin && fabs(stepy) < stepmin ) break; |
1991 | |
1992 | chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0+stepx, y0+stepy, chi00 ); |
1993 | |
1994 | |
1995 | |
1996 | if( chi00 >= chi0 ) break; |
1997 | |
1998 | chi0 = chi00; |
1999 | x0 = x0+stepx; |
2000 | y0 = y0+stepy; |
2001 | |
2002 | rounds++; |
2003 | if( rounds > 10000 ) { cout << "max rounds" << endl; break; } |
2004 | |
2005 | } |
2006 | |
2007 | if( chi0 < chisq0 ) { |
2008 | x1 = x0; |
2009 | y1 = y0; |
2010 | chisq = chi0/dof; |
2011 | } |
2012 | |
2013 | |
2014 | |
2015 | |
2016 | if( chisq <= chimem ) return; |
2017 | |
2018 | chiold = chisq; |
2019 | tgamma_hyc( nadc, ia, id, nzero, iaz, chisq, ee, xx, yy, e2, x2, y2 ); |
2020 | |
2021 | if( e2 > 0. ) { |
2022 | |
2023 | |
2024 | |
2025 | d2 = pow((xx-x2)*xsize2.09, 2.0) + pow((yy-y2)*ysize2.09, 2.0); |
2026 | xm2 = ee*e2*d2; |
2027 | |
2028 | if( xm2 > 0. ) xm2 = sqrt(xm2)/1270.*0.1; |
2029 | |
2030 | if( xm2 > xm2cut*xsize2.09 ) { |
2031 | e1 = ee; |
2032 | x1 = xx; |
2033 | y1 = yy; |
2034 | } else { |
2035 | e2 = 0.; |
2036 | chisq = chiold; |
2037 | } |
2038 | } |
2039 | |
2040 | return; |
2041 | } |
2042 | |
2043 | |
2044 | |
2045 | |
2046 | |
2047 | |
2048 | |
2049 | |
2050 | |
2051 | |
2052 | |
2053 | void DCCALShower_factory::fill_zeros( int nadc, vector<int> ia, int &nneib, vector<int> &iaz ) |
2054 | { |
2055 | |
2056 | int ix, iy, nneibnew; |
2057 | vector<int> ian; |
2058 | |
2059 | nneib = 0; |
2060 | for( int ii = 0; ii < nadc; ii++ ) { |
2061 | ix = ia[ii]/100; |
2062 | iy = ia[ii] - ix*100; |
2063 | |
2064 | if( ix > 1 ) { |
2065 | nneib = nneib+1; |
2066 | ian.push_back(iy + (ix-1)*100); |
2067 | |
2068 | if( iy > 1 ) { |
2069 | nneib = nneib+1; |
2070 | ian.push_back(iy-1 + (ix-1)*100); |
2071 | } |
2072 | if( iy < MROW12 ) { |
2073 | nneib = nneib+1; |
2074 | ian.push_back(iy+1 + (ix-1)*100); |
2075 | } |
2076 | } |
2077 | if( ix < MCOL12 ) { |
2078 | nneib = nneib+1; |
2079 | ian.push_back(iy + (ix+1)*100); |
2080 | |
2081 | if( iy > 1 ) { |
2082 | nneib = nneib+1; |
2083 | ian.push_back(iy-1 + (ix+1)*100); |
2084 | } |
2085 | if( iy < MROW12 ) { |
2086 | nneib = nneib+1; |
2087 | ian.push_back(iy+1 + (ix+1)*100); |
2088 | } |
2089 | } |
2090 | |
2091 | if( iy > 1 ) { |
2092 | nneib = nneib+1; |
2093 | ian.push_back(iy-1 + ix*100); |
2094 | } |
2095 | if( iy < MROW12 ) { |
2096 | nneib = nneib+1; |
2097 | ian.push_back(iy+1 + ix*100); |
2098 | } |
2099 | } |
2100 | |
2101 | |
2102 | for( int ii = 0; ii < nneib; ii++ ) { |
2103 | for( int jj = 0; jj < nadc; jj++ ) { |
2104 | if( ia[jj] == ian[ii] ) ian[ii] = -1; |
2105 | } |
2106 | } |
2107 | |
2108 | for( int ii = 0; ii < nneib; ii++ ) { |
2109 | if( ian[ii] == -1 ) continue; |
2110 | for( int jj = ii+1; jj < nneib; jj++ ) { |
2111 | if( ian[jj] == ian[ii] ) ian[jj] = -1; |
2112 | } |
2113 | } |
2114 | |
2115 | nneibnew = 0; |
2116 | for( int ii = 0; ii < nneib; ii++ ) { |
2117 | ix = ian[ii]/100; |
2118 | iy = ian[ii] - ix*100; |
2119 | if( ian[ii] != -1 ) { |
2120 | if( stat_ch[iy-1][ix-1] == 0 ) { |
2121 | nneibnew = nneibnew+1; |
2122 | iaz.push_back( ian[ii] ); |
2123 | } |
2124 | } |
2125 | } |
2126 | nneib = nneibnew; |
2127 | |
2128 | |
2129 | return; |
2130 | } |
2131 | |
2132 | |
2133 | |
2134 | |
2135 | |
2136 | |
2137 | |
2138 | |
2139 | |
2140 | |
2141 | |
2142 | void DCCALShower_factory::mom1_pht( int nadc, vector<int> ia, vector<int> id, |
2143 | int nzero, vector<int> iaz, double &e1, double &x1, double &y1 ) |
2144 | { |
2145 | |
2146 | |
2147 | int ix, iy; |
2148 | double dx, dy; |
2149 | |
2150 | double a; |
2151 | double corr; |
2152 | |
2153 | |
2154 | |
2155 | e1 = 0.; |
2156 | x1 = 0.; |
2157 | y1 = 0.; |
2158 | |
2159 | if( nadc <= 0 ) return; |
2160 | for( int ii = 0; ii < nadc; ii++ ) { |
2161 | a = static_cast<double>(id[ii]); |
2162 | ix = ia[ii]/100; |
2163 | iy = ia[ii] - ix*100; |
2164 | e1 = e1 + a; |
2165 | x1 = x1 + a*static_cast<double>(ix); |
2166 | y1 = y1 + a*static_cast<double>(iy); |
2167 | } |
2168 | if( e1 <= 0. ) return; |
2169 | x1 = x1/e1; |
2170 | y1 = y1/e1; |
2171 | corr = 0.; |
2172 | for( int ii = 0; ii < nadc; ii++ ) { |
2173 | ix = ia[ii]/100; |
2174 | iy = ia[ii] - ix*100; |
2175 | dx = static_cast<double>(ix) - x1; |
2176 | dy = static_cast<double>(iy) - y1; |
2177 | corr = corr + cell_hyc( dx, dy ); |
2178 | } |
2179 | |
2180 | if( COUNT_ZERO1 ) { |
2181 | for( int ii = 0; ii < nzero; ii++ ) { |
2182 | ix = iaz[ii]/100; |
2183 | iy = iaz[ii] - ix*100; |
2184 | dx = static_cast<double>(ix) - x1; |
2185 | dy = static_cast<double>(iy) - y1; |
2186 | corr = corr + cell_hyc( dx, dy ); |
2187 | } |
2188 | } |
2189 | |
2190 | corr = corr/1.006; |
2191 | if( SHOWER_DEBUG ) cout << "e, corr = " << e1 << ", " << corr << endl; |
2192 | |
2193 | |
2194 | if( corr < 0.8 ) { |
2195 | if( SHOWER_DEBUG ) { |
2196 | cout << "corr = " << corr << ", " << e1 << ", " << x1 << ", " << y1; |
2197 | cout << "! Too many around central hole!" << endl; |
2198 | } |
2199 | corr = 0.8; |
2200 | } else if( corr > 1.0 ) { |
2201 | corr = 1.; |
2202 | } |
2203 | e1 = e1/corr; |
2204 | |
2205 | |
2206 | return; |
2207 | } |
2208 | |
2209 | |
2210 | |
2211 | |
2212 | |
2213 | |
2214 | |
2215 | |
2216 | |
2217 | |
2218 | |
2219 | void DCCALShower_factory::chisq1_hyc( int nadc, vector<int> ia, vector<int> id, |
2220 | int nneib, vector<int> iaz, double e1, double x1, double y1, double &chisq ) |
2221 | { |
2222 | |
2223 | |
2224 | int ix, iy; |
2225 | double dx, dy; |
2226 | double fa, fcell; |
2227 | |
2228 | |
2229 | |
2230 | |
2231 | chisq = 0.; |
2232 | |
2233 | for( int ii = 0; ii < nadc; ii++ ) { |
2234 | fa = static_cast<double>(id[ii]); |
2235 | ix = ia[ii]/100; |
2236 | iy = ia[ii] - ix*100; |
2237 | dx = x1 - static_cast<double>(ix); |
2238 | dy = y1 - static_cast<double>(iy); |
2239 | if( e1 != 0. ) { |
2240 | if( fabs(dx) <= 6.0 && fabs(dy) <= 6.0 ) { |
2241 | fcell = cell_hyc( dx, dy ); |
2242 | chisq = chisq + e1*pow((fcell-(fa/e1)), 2.)/sigma2(dx, dy, fcell, e1); |
2243 | } |
2244 | } else { |
2245 | chisq = chisq + fa*fa/9.; |
2246 | |
2247 | } |
2248 | } |
2249 | |
2250 | for(int ii = 0; ii < nneib; ii++ ) { |
2251 | ix = iaz[ii]/100; |
2252 | iy = iaz[ii] - ix*100; |
2253 | dx = x1 - static_cast<double>(ix); |
2254 | dy = y1 - static_cast<double>(iy); |
2255 | if( e1 != 0. ) { |
2256 | if( fabs(dx) < 6.0 && fabs(dy) < 6.0 ) { |
2257 | fcell = cell_hyc( dx, dy ); |
2258 | chisq = chisq + e1*fcell*fcell/sigma2(dx, dy, fcell, e1); |
2259 | } |
2260 | } else { |
2261 | chisq = chisq + id[ii]*id[ii]/9.; |
2262 | |
2263 | } |
2264 | } |
2265 | |
2266 | |
2267 | return; |
2268 | } |
2269 | |
2270 | |
2271 | |
2272 | |
2273 | |
2274 | |
2275 | |
2276 | |
2277 | |
2278 | |
2279 | |
2280 | double DCCALShower_factory::sigma2( double dx, double dy, double fc, double e ) |
2281 | { |
2282 | double sig2; |
2283 | double alp = 0.816; |
2284 | double bet1 = 32.1; |
2285 | double bet2 = 1.72; |
2286 | |
2287 | double r = dx*dx + dy*dy; |
2288 | if( r > 25. ) { |
2289 | sig2 = 100.; |
2290 | return sig2; |
2291 | } |
2292 | |
2293 | sig2 = alp*fc + (bet1 + bet2*sqrt(e/100.))*d2c( dx, dy ) + 0.2/(e/100.); |
2294 | if( TEST_P1 ) sig2 = sig2/pow(0.0001*e, 0.166); |
2295 | sig2 *= 100.; |
2296 | |
2297 | return sig2; |
2298 | } |
2299 | |
2300 | |
2301 | double DCCALShower_factory::d2c( double dx, double dy ) |
2302 | { |
2303 | int i, j; |
2304 | double ax, ay, wx, wy; |
2305 | double d2c; |
2306 | |
2307 | ax = fabs(dx*100.); |
2308 | ay = fabs(dy*100.); |
2309 | i = int(ax); |
2310 | j = int(ay); |
2311 | |
2312 | if( (i < 499) && (j < 499) && (i >= 0) && (j >= 0) ) { |
2313 | |
2314 | |
2315 | wx = ax-static_cast<double>(i); |
2316 | wy = ay-static_cast<double>(j); |
2317 | |
2318 | d2c = ad2c[i][j] * (1.-wx) * (1.-wy) + |
2319 | ad2c[i+1][j] * wx * (1.-wy) + |
2320 | ad2c[i][j+1] * (1.-wx) * wy + |
2321 | ad2c[i+1][j+1] * wx * wy; |
2322 | |
2323 | } else d2c = 1.; |
2324 | |
2325 | return d2c; |
2326 | } |
2327 | |
2328 | |
2329 | double DCCALShower_factory::cell_hyc( double dx, double dy ) |
2330 | { |
2331 | int i, j; |
2332 | double ax, ay, wx, wy; |
2333 | double cell_hyc; |
2334 | |
2335 | ax = fabs(dx*100.); |
2336 | ay = fabs(dy*100.); |
2337 | i = static_cast<int>(ax); |
2338 | j = static_cast<int>(ay); |
2339 | |
2340 | if( i < 499 && j < 499 && i >= 0 && j >= 0 ) { |
2341 | |
2342 | wx = ax-static_cast<double>(i); |
2343 | wy = ay-static_cast<double>(j); |
2344 | |
2345 | cell_hyc = acell[i][j] * (1.-wx) * (1.-wy) + |
2346 | acell[i+1][j] * wx * (1.-wy) + |
2347 | acell[i][j+1] * (1.-wx) * wy + |
2348 | acell[i+1][j+1] * wx * wy; |
2349 | |
2350 | } else cell_hyc = 0.; |
2351 | |
2352 | return cell_hyc; |
2353 | } |
2354 | |
2355 | |
2356 | |
2357 | |
2358 | |
2359 | |
2360 | |
2361 | |
2362 | |
2363 | |
2364 | |
2365 | |
2366 | void DCCALShower_factory::tgamma_hyc( int nadc, vector<int> ia, vector<int> id, |
2367 | int nzero, vector<int> iaz, double &chisq, double &e1, double &x1, double &y1, |
2368 | double &e2, double &x2, double &y2 ) |
2369 | { |
2370 | |
2371 | |
2372 | int ix, iy; |
2373 | int dof; |
2374 | |
2375 | double dx, dy; |
2376 | double dxy; |
2377 | double dxc, dyc; |
2378 | double dx0, dy0; |
2379 | double dx1, dy1; |
2380 | double dx2, dy2; |
2381 | double u, r, rsq, rsq2; |
2382 | double epsc, eps0, eps1, eps2; |
2383 | double stepmin, epsmax; |
2384 | double delch; |
2385 | double step; |
2386 | double cosi, scal; |
2387 | double chisq2, chisqc; |
2388 | double dchi, dchida, dchi0; |
2389 | double a1, a2; |
2390 | double ex; |
2391 | double e1c, x1c, y1c; |
2392 | double e2c, x2c, y2c; |
2393 | double gr, gre, grx, gry; |
2394 | double grc; |
2395 | double grec, grxc, gryc; |
2396 | double gx1, gy1; |
2397 | double gx2, gy2; |
2398 | double e0, x0, y0; |
2399 | double xx, yy, yx; |
2400 | |
2401 | double f1c, f2c, f1x, f2x, f1y, f2y; |
2402 | double chisqt, chisqt0, chisqt1, chisqt2; |
2403 | double chisqtx1, chisqtx2, chisqty1, chisqty2; |
2404 | double dchidax1, dchidax2, dchiday1, dchiday2; |
2405 | |
2406 | stepmin = 0.5; |
2407 | epsmax = 0.9999; |
2408 | delch = 10.; |
2409 | |
2410 | |
2411 | |
2412 | yx = 0.; |
2413 | grx = 0.; gry = 0.; |
2414 | gre = 0.; gr = 0.; |
2415 | dx0 = 0.; dy0 = 0.; |
2416 | eps0 = 0.; |
2417 | |
2418 | e2 = 0.; |
2419 | x2 = 0.; |
2420 | y2 = 0.; |
2421 | |
2422 | if( nadc <= 0 ) return; |
2423 | |
2424 | mom2_pht( nadc, ia, id, nzero, iaz, e0, x0, y0, xx, yy, yx ); |
2425 | |
2426 | if( e0 <= 0 ) return; |
2427 | |
2428 | |
2429 | |
2430 | dxy = xx - yy; |
2431 | rsq2 = dxy*dxy + 4.*yx*yx; |
2432 | if( rsq2 < 1.e-20 ) rsq2 = 1.e-20; |
2433 | rsq = sqrt(rsq2); |
2434 | dxc = -sqrt((rsq+dxy)*2.); |
2435 | dyc = sqrt((rsq-dxy)*2.); |
2436 | if( yx >= 0. ) dyc = -dyc; |
2437 | r = sqrt(dxc*dxc + dyc*dyc); |
2438 | epsc = 0.; |
2439 | for( int ii = 0; ii < nadc; ii++ ) { |
2440 | ix = ia[ii]/100; |
2441 | iy = ia[ii] - ix*100; |
2442 | dx = static_cast<double>(ix) - x0; |
2443 | dy = static_cast<double>(iy) - y0; |
2444 | u = dx*dxc/r + dy*dyc/r; |
2445 | epsc = epsc - 0.01*id[ii]*u*fabs(u); |
2446 | } |
2447 | if (rsq != 0) |
2448 | epsc = epsc/(0.01*e0*rsq); |
2449 | if( epsc > 0.8 ) epsc = 0.8; |
2450 | if( epsc < -0.8 ) epsc = -0.8; |
2451 | dxc = dxc/sqrt(1.-(epsc*epsc)); |
2452 | dyc = dyc/sqrt(1.-(epsc*epsc)); |
2453 | |
2454 | |
2455 | |
2456 | |
2457 | |
2458 | step = 0.1; |
2459 | cosi = 0.0; |
| Value stored to 'cosi' is never read |
2460 | chisq2 = 1.e35; |
2461 | |
2462 | for( int iter = 0; iter < 100; iter++ ) { |
2463 | |
2464 | c3to5_pht( e0, x0, y0, epsc, dxc, dyc, e1c, x1c, y1c, e2c, x2c, y2c ); |
2465 | |
2466 | eps1 = (1. + epsc)/2.; |
2467 | eps2 = (1. - epsc)/2.; |
2468 | chisqc = 0.; |
2469 | |
2470 | for( int ii = 0; ii < nadc+nzero; ii++ ) { |
2471 | |
2472 | if( ii < nadc ) { |
2473 | ex = static_cast<double>(id[ii]); |
2474 | ix = ia[ii]/100; |
2475 | iy = ia[ii] - ix*100; |
2476 | } else { |
2477 | ex = 0.; |
2478 | ix = iaz[ii-nadc]/100; |
2479 | iy = iaz[ii-nadc] - ix*100; |
2480 | } |
2481 | |
2482 | dx1 = x1c - static_cast<double>(ix); |
2483 | dy1 = y1c - static_cast<double>(iy); |
2484 | dx2 = x2c - static_cast<double>(ix); |
2485 | dy2 = y2c - static_cast<double>(iy); |
2486 | |
2487 | f1c = cell_hyc( dx1, dy1 ); |
2488 | f2c = cell_hyc( dx2, dy2 ); |
2489 | |
2490 | chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt ); |
2491 | chisqc += chisqt; |
2492 | |
2493 | } |
2494 | |
2495 | if( chisqc >= chisq2 ) { |
2496 | |
2497 | if( iter > 0 ) { |
2498 | dchi = chisqc - chisq2; |
2499 | dchi0 = gr*step; |
2500 | step = 0.5*step/sqrt(1. + dchi/dchi0); |
2501 | } |
2502 | step = 0.5*step; |
2503 | |
2504 | } else { |
2505 | |
2506 | grec = 0.; |
2507 | grxc = 0.; |
2508 | gryc = 0.; |
2509 | |
2510 | for( int ii = 0; ii < nadc+nzero; ii++ ) { |
2511 | |
2512 | if( ii < nadc ) { |
2513 | ex = static_cast<double>(id[ii]); |
2514 | ix = ia[ii]/100; |
2515 | iy = ia[ii] - ix*100; |
2516 | } else { |
2517 | ex = 0.; |
2518 | ix = iaz[ii-nadc]/100; |
2519 | iy = iaz[ii-nadc] - ix*100; |
2520 | } |
2521 | |
2522 | dx1 = x1c - static_cast<double>(ix); |
2523 | dy1 = y1c - static_cast<double>(iy); |
2524 | dx2 = x2c - static_cast<double>(ix); |
2525 | dy2 = y2c - static_cast<double>(iy); |
2526 | |
2527 | f1c = cell_hyc( dx1, dy1 ); |
2528 | f2c = cell_hyc( dx2, dy2 ); |
2529 | |
2530 | a1 = e1c*f1c; |
2531 | a2 = e2c*f2c; |
2532 | |
2533 | |
2534 | chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt0 ); |
2535 | chisq2t_hyc( ex, e1c+1., dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt1 ); |
2536 | chisq2t_hyc( ex, e1c, dx1, dy1, e2c+1., dx2, dy2, f1c, f2c, chisqt2 ); |
2537 | |
2538 | f1x = cell_hyc( dx1+0.05, dy1 ); |
2539 | f2x = cell_hyc( dx2+0.05, dy2 ); |
2540 | f1y = cell_hyc( dx1, dy1+0.05 ); |
2541 | f2y = cell_hyc( dx2, dy2+0.05 ); |
2542 | |
2543 | chisq2t_hyc( ex, e1c, dx1+0.05, dy1, e2c, dx2, dy2, f1x, f2c, chisqtx1 ); |
2544 | chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2+0.05, dy2, f1c, f2x, chisqtx2 ); |
2545 | chisq2t_hyc( ex, e1c, dx1, dy1+0.05, e2c, dx2, dy2, f1y, f2c, chisqty1 ); |
2546 | chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2+0.05, f1c, f2y, chisqty2 ); |
2547 | |
2548 | dchidax1 = 20.*(chisqtx1 - chisqt0); |
2549 | dchidax2 = 20.*(chisqtx2 - chisqt0); |
2550 | dchiday1 = 20.*(chisqty1 - chisqt0); |
2551 | dchiday2 = 20.*(chisqty2 - chisqt0); |
2552 | dchida = 0.5*(chisqt1 + chisqt2 - chisqt0); |
2553 | |
2554 | gx1 = (e1c*f1x - a1)*dchidax1; |
2555 | gx2 = (e2c*f2x - a2)*dchidax2; |
2556 | gy1 = (e1c*f1y - a1)*dchiday1; |
2557 | gy2 = (e2c*f2y - a2)*dchiday2; |
2558 | |
2559 | grec += dchida*(f1c-f2c)*e0 - ( (gx1+gx2)*dxc + (gy1+gy2)*dyc ); |
2560 | grxc += gx1*eps2 - gx2*eps1; |
2561 | gryc += gy1*eps2 - gy2*eps1; |
2562 | |
2563 | } |
2564 | |
2565 | grc = sqrt( grec*grec + grxc*grxc + gryc*gryc ); |
2566 | if( grc < 1.e-6 ) grc = 1.e-6; |
2567 | if( iter > 0 ) { |
2568 | cosi = ( gre*grec + grx*grxc + gry*gryc )/(gr*grc); |
2569 | scal = fabs((gr/grc) - cosi); |
2570 | if( scal < 0.1 ) scal = 0.1; |
2571 | step = step/scal; |
2572 | } |
2573 | chisq2 = chisqc; |
2574 | eps0 = epsc; |
2575 | dx0 = dxc; |
2576 | dy0 = dyc; |
2577 | gre = grec; |
2578 | grx = grxc; |
2579 | gry = gryc; |
2580 | gr = grc; |
2581 | } |
2582 | epsc = eps0 - step*gre/gr; |
2583 | while( fabs(epsc) > epsmax ) { |
2584 | step = step/2.; |
2585 | epsc = eps0 - step*gre/gr; |
2586 | } |
2587 | dxc = dx0 - step*grx/gr; |
2588 | dyc = dy0 - step*gry/gr; |
2589 | if( step*gr < stepmin ) break; |
2590 | } |
2591 | |
2592 | if( (chisq*(nadc+nzero-2) - chisq2) < delch ) return; |
2593 | dof = nzero+nadc-5; |
2594 | if( dof < 1 ) dof = 1; |
2595 | chisq = chisq2/dof; |
2596 | |
2597 | c3to5_pht( e0, x0, y0, eps0, dx0, dy0, e1, x1, y1, e2, x2, y2 ); |
2598 | |
2599 | return; |
2600 | } |
2601 | |
2602 | |
2603 | |
2604 | |
2605 | |
2606 | |
2607 | |
2608 | |
2609 | |
2610 | |
2611 | |
2612 | |
2613 | void DCCALShower_factory::mom2_pht( int nadc, vector<int> ia, vector<int> id, |
2614 | int nzero, vector<int> iaz, double &a0, double &x0, double &y0, |
2615 | double &xx, double &yy, double &yx) |
2616 | { |
2617 | |
2618 | |
2619 | int ix, iy; |
2620 | double dx, dy; |
2621 | double a; |
2622 | double corr; |
2623 | |
2624 | |
2625 | |
2626 | a0 = 0.; |
2627 | x0 = 0.; |
2628 | y0 = 0.; |
2629 | xx = 0.; |
2630 | yy = 0.; |
2631 | yx = 0.; |
2632 | |
2633 | if( nadc <= 0 ) return; |
2634 | for( int ii = 0; ii < nadc; ii++ ) { |
2635 | a = static_cast<double>(id[ii]); |
2636 | ix = ia[ii]/100; |
2637 | iy = ia[ii] - ix*100; |
2638 | a0 = a0 + a; |
2639 | x0 = x0 + a*static_cast<double>(ix); |
2640 | y0 = y0 + a*static_cast<double>(iy); |
2641 | } |
2642 | if( a0 <= 0. ) return; |
2643 | x0 = x0/a0; |
2644 | y0 = y0/a0; |
2645 | |
2646 | for( int ii = 0; ii < nadc; ii++ ) { |
2647 | a = static_cast<double>(id[ii])/a0; |
2648 | ix = ia[ii]/100; |
2649 | iy = ia[ii] - ix*100; |
2650 | dx = static_cast<double>(ix) - x0; |
2651 | dy = static_cast<double>(iy) - y0; |
2652 | xx = xx + a*dx*dx; |
2653 | yy = yy + a*dy*dy; |
2654 | yx = yx + a*dx*dy; |
2655 | } |
2656 | |
2657 | corr = 0.; |
2658 | for( int ii = 0; ii < nadc; ii++ ) { |
2659 | ix = ia[ii]/100; |
2660 | iy = ia[ii] - ix*100; |
2661 | dx = static_cast<double>(ix) - x0; |
2662 | dy = static_cast<double>(iy) - y0; |
2663 | corr = corr + cell_hyc( dx, dy ); |
2664 | } |
2665 | if( COUNT_ZERO1 ) { |
2666 | for( int ii = 0; ii < nzero; ii++ ) { |
2667 | ix = iaz[ii]/100; |
2668 | iy = iaz[ii] - ix*100; |
2669 | dx = static_cast<double>(ix) - x0; |
2670 | dy = static_cast<double>(iy) - y0; |
2671 | corr = corr + cell_hyc( dx, dy ); |
2672 | } |
2673 | } |
2674 | |
2675 | corr = corr/1.006; |
2676 | if( corr < 0.8 ) { |
2677 | corr = 0.8; |
2678 | } else if( corr > 1. ) { |
2679 | corr = 1.0; |
2680 | } |
2681 | a0 = a0/corr; |
2682 | |
2683 | |
2684 | return; |
2685 | } |
2686 | |
2687 | |
2688 | |
2689 | |
2690 | |
2691 | |
2692 | |
2693 | |
2694 | |
2695 | |
2696 | |
2697 | |
2698 | void DCCALShower_factory::c3to5_pht( double e0, double x0, double y0, double eps, |
2699 | double dx, double dy, double &e1, double &x1, double &y1, double &e2, |
2700 | double &x2, double &y2 ) |
2701 | { |
2702 | e1 = e0*(1.+eps)/2.; |
2703 | e2 = e0 - e1; |
2704 | x1 = x0 + 0.5*dx*(1.-eps); |
2705 | y1 = y0 + 0.5*dy*(1.-eps); |
2706 | x2 = x0 - 0.5*dx*(1.+eps); |
2707 | y2 = y0 - 0.5*dy*(1.+eps); |
2708 | |
2709 | return; |
2710 | } |
2711 | |
2712 | |
2713 | |
2714 | |
2715 | |
2716 | |
2717 | |
2718 | |
2719 | |
2720 | |
2721 | |
2722 | |
2723 | void DCCALShower_factory::chisq2t_hyc( double ecell, double e1, double dx1, double dy1, |
2724 | double e2, double dx2, double dy2, double f1, double f2, double &chisqt ) |
2725 | { |
2726 | double s; |
2727 | double p1, p2; |
2728 | |
2729 | if( TEST_P1 ) { |
2730 | p1 = pow(0.0001*e1, 0.166); |
2731 | p2 = pow(0.0001*e2, 0.166); |
2732 | } else { |
2733 | p1 = 1.; |
2734 | p2 = 1.; |
2735 | } |
2736 | |
2737 | if( e1 != 0. && e2 != 0. ) |
2738 | s = e1*sigma2(dx1, dy1, f1, e1)/p1 + e2*sigma2(dx2, dy2, f2, e2)/p2; |
2739 | else if( e1 == 0. && e2 == 0. ) |
2740 | s = 90000.; |
2741 | else if( e1 == 0. ) |
2742 | s = e2*sigma2(dx2, dy2, f2, e2)/p2; |
2743 | else |
2744 | s = e1*sigma2(dx1, dy1, f1, e1)/p1; |
2745 | |
2746 | chisqt = pow((e1*f1 + e2*f2 - ecell), 2.)/s; |
2747 | |
2748 | return; |
2749 | } |
2750 | |
2751 | |
2752 | |
2753 | |
2754 | |
2755 | |
2756 | |
2757 | |
2758 | |
2759 | |
2760 | |
2761 | |
2762 | void DCCALShower_factory::ucopy1( vector<int> &ia, vector<int> &iwork, int start, int length ) |
2763 | { |
2764 | |
2765 | for( int ii = 0; ii < length; ii++ ) { |
2766 | iwork.push_back( ia[start+ii] ); |
2767 | } |
2768 | |
2769 | |
2770 | return; |
2771 | } |
2772 | |
2773 | |
2774 | |
2775 | |
2776 | |
2777 | |
2778 | |
2779 | |
2780 | |
2781 | void DCCALShower_factory::ucopy2( vector<int> &ia, int start1, int start2, int length ) |
2782 | { |
2783 | |
2784 | vector<int> work; |
2785 | for( int ii = 0; ii < length; ii++ ) { |
2786 | work.push_back( ia[start1+ii] ); |
2787 | } |
2788 | |
2789 | for( int ii = 0; ii < length; ii++ ) { |
2790 | ia[start2+ii] = work[ii]; |
2791 | } |
2792 | |
2793 | |
2794 | return; |
2795 | } |
2796 | |
2797 | |
2798 | |
2799 | |
2800 | |
2801 | |
2802 | |
2803 | |
2804 | |
2805 | |
2806 | void DCCALShower_factory::ucopy3( vector<int> &iwork, vector<int> &ia, int start, int length ) |
2807 | { |
2808 | |
2809 | for( int ii = 0; ii < length; ii++ ) { |
2810 | ia[start+ii] = iwork[ii]; |
2811 | } |
2812 | |
2813 | |
2814 | return; |
2815 | } |
2816 | |
2817 | |
2818 | |
2819 | |
2820 | |
2821 | |
2822 | |
2823 | |
2824 | |
2825 | |
2826 | |
2827 | |
2828 | |
2829 | |
2830 | |
2831 | |
2832 | |
2833 | |
2834 | |
2835 | |
2836 | |
2837 | |