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