Bug Summary

File:/volatile/halld/gluex/nightly/2024-04-14/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/CCAL/DCCALShower_factory.cc
Location:line 2459, column 2
Description:Value stored to 'cosi' is never read

Annotated Source Code

1
2/*
3 * File: DCCALHit_factory.cc
4 *
5 * Created on 11/25/18 by A.S.
6 * use structure similar to FCAL
7 */
8
9#include "CCAL/DCCALShower_factory.h"
10
11
12static mutex CCAL_MUTEX;
13//static bool CCAL_PROFILE_LOADED = false;
14
15
16//==========================================================
17//
18// Constructor
19//
20//==========================================================
21
22DCCALShower_factory::DCCALShower_factory()
23{
24 // Set defaults:
25
26 MIN_CLUSTER_BLOCK_COUNT = 2;
27 MAX_HITS_FOR_CLUSTERING = 80;
28 MIN_CLUSTER_SEED_ENERGY = 0.035; /* [GeV] */
29 MIN_CLUSTER_ENERGY = 0.05; /* [GeV] */
30 MAX_CLUSTER_ENERGY = 15.9; /* [GeV] */
31 TIME_CUT = 15.0; /* [ns] */
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; ///< >0 once off info ; >2 event by event ; >3 everything
62 gPARMS->SetDefaultParameter("DFCALShower:VERBOSE", VERBOSE, "Verbosity level for DFCALShower objects and factories");
63}
64
65
66
67
68//==========================================================
69//
70// brun
71//
72//==========================================================
73
74jerror_t DCCALShower_factory::brun(JEventLoop *locEventLoop, int32_t runnumber)
75{
76 // Only print messages for one thread whenever run number change
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 //----------- Read in shower profile data -----------//
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 // follow similar procedure as other resources (DMagneticFieldMapFineMesh)
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 // check to see if we actually have a file
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 //---------- Initialize channel status array ----------//
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 //---------- Read shower timewalk parameters ----------//
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 //--------- Read in non-linearity parameters ----------//
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// evnt
260//
261//==========================================================
262
263jerror_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 //------- Get the CCALHits and organize hit pattern -------//
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; // will hold all clusters
293 vector< cluster_t > clusterStorage; // will hold the constituents of every cluster
294
295
296
297
298 //------- Call island clusterizer on each pattern -------//
299
300
301 for( int ipat = 0; ipat < n_patterns; ipat++ )
302 {
303
304
305 /*----- prepare data in dimensionless format -----*/
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 /*------------- call to island -------------*/
330
331 vector< gamma_t > gammas; gammas.clear();// Output of main_island (holds reconstructed photons)
332 main_island( ia, id, gammas );
333
334
335 /*------------ post-island processing ------------*/
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 } // end looping over hit patterns
345
346
347
348
349 //-------------- Fill DCCALShower Object --------------//
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 //shower->ExyztCovariance(i, j)
381 float xerr = 0.1016;// / sqrt(ccalClusters[k].E) + 0.2219; // HARD CODED VALUE!!!!
382 float yerr = 0.1016;// / sqrt(ccalClusters[k].E) + 0.2219; // HARD CODED VALUE!!!!
383 float zerr = 2.5; // HARD CODED VALUE!!!!
384 float terr = 0.2; // HARD CODED VALUE!!!!
385 float eerr = 0.5;//pow(ccalClusters[k].E, 2) * (0.01586 / ccalClusters[k].E + 0.0002342 / pow(ccalClusters[k].E, 2) + 1.696e-6); // HARD CODED VALUE!!!!
386 //copy xyz errors into covariance matrix
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// getHitPatterns
436//
437//==========================================================
438
439void DCCALShower_factory::getHitPatterns( vector< const DCCALHit* > hitarray,
440 vector< vector< const DCCALHit* > > &hitPatterns )
441{
442
443 /*------------------------------
444
445 Method for sorting hit patterns:
446
447 1. Find hit with largest energy.
448 2. Sort hit vector in order of increasing time-difference from this hit
449 3. Push all hits within 15 ns of max hit to a new vector (first few elements of sorted vec)
450 4. Erase the elements that were moved from the previous vector
451 - maybe steps 3&4 can be done simultaneously?
452 5. Push the new vector back to the hitPatterns vector
453 6. Repeat steps 1-5 until no hits remain in original hitarray
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// sortByTime
522//
523//==========================================================
524
525void 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; // nothing to sort
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 } // end if statement
562
563 } // end loop over hits
564
565 return;
566}
567
568
569
570
571//==========================================================
572//
573// cleanHitPattern
574//
575//==========================================================
576
577void 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// processShowers
611//
612//==========================================================
613
614void 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 //------------- Do some post-island processing -------------//
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; // stores cluster parameters
630 cluster_t locClusterStorage; // stores hit information of cluster cells
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 // check that shower is not just a single module and that energy is reasonable:
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 //------------ Find cell with max energy ------------//
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 //----------- Loop over constituent hits ------------//
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 // The shower position is calculated using logarithmic weighting:
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++ ) // zero the rest
729 locClusterStorage.id[j] = -1;
730
731 double weightedTime = getEnergyWeightedTime( locClusterStorage, dime );
732 double showerTime = getCorrectedTime( weightedTime, e );
733
734
735
736 //------- Get position at surface of Calorimeter -------//
737
738 DVector3 vertex(0.0, 0.0, m_zTarget); // for now, use center of target as vertex
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 //-------- Calculate nonlinear-corrected energy --------//
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 //-------- Get energy resolution (needs updating) --------//
777
778 double se = sqrt( 0.9*0.9*ecorr*ecorr + 2.5*2.5*ecorr + 1.0 );
779 // from HYCAL reconstruction, need to tune
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 //----------------- Fill cluster bank -----------------//
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// getEnergyWeightedTime
822//
823//==========================================================
824
825double 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// getCorrectedTime
846//
847//==========================================================
848
849double DCCALShower_factory::getCorrectedTime( double time, double energy )
850{
851 // timewalk correction:
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// getShowerDepth
870//
871//==========================================================
872
873double 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// getCorrectedEnergy
888//
889//==========================================================
890
891double 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 //if( energy < 0.5 || energy > 12. ) return energy;
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// nonlin_func
930//
931//==========================================================
932
933double 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// main_island()
967//
968//==========================================================
969
970void DCCALShower_factory::main_island( vector<int> &ia, vector<int> &id, vector<gamma_t> &gammas )
971{
972
973 //----------------------------------------------
974 // Initialize some useful variables:
975
976 int nhits; // number of hits in detector
977 int ncl; // number of clusters
978 vector<int> lencl; // length of a cluster
979
980
981 nhits = static_cast<int>(ia.size());
982 if( nhits == 0 ) return;
983
984 /*
985 The vector 'ia' will hold the addresses of modules that were hit.
986 They're defined as 100*(i+1)+(j+1) where i is column and j is row of the hit module.
987 Row 0, column 0 is bottom right corner of CCAL (looking upstream).
988
989 The 'id' vector holds the energies of the hit modules (in units of 0.1 MeVs).
990 */
991
992
993 //------------------ Search for clusters: ------------------//
994
995
996 ncl = clus_hyc( nhits, ia, id, lencl );
997
998 /*
999 At this point, the vectors ia and id are sorted such that they are grouped
1000 together by simply-connected clusters. For example, the first lencl[0] elements of ia give
1001 the addresses of the elements in the first cluster. The next lencl[1] elements
1002 after that give the addresses of the elements in the second cluster.
1003
1004 The clusters here are just created by grouping everything that is physically
1005 next to each other into one cluster - all simply-connected cells are joined in a cluster.
1006 */
1007
1008
1009 //-------------------- Process Cluster: --------------------//
1010
1011 int nadcgam = 0; // number of found clusters
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; // addresses of current cluster
1024 vector< int > icl_d; // energies of current cluster
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 //-------------------- Prepare gammas for final processing --------------------//
1060
1061
1062 // convert position to units of cm:
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 // sort gammas in order of decreasing energy:
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// clus_hyc
1110//
1111//==========================================================
1112
1113int DCCALShower_factory::clus_hyc( int nw, vector<int> &ia, vector<int> &id, vector<int> &lencl )
1114{
1115
1116 //---------------- Local Declarations ---------------//
1117
1118 int maxcl = 200; // maximumum number of clusters allowed
1119 int ncl; // number of clusters
1120 int next, iak; //
1121 int ib, ie; //
1122 int ias, iaf; //
1123 int last, lastcl; //
1124 int leng; //
1125 //
1126 int loclencl[200]; // stores the lengths of clusters locally
1127
1128
1129 //---------------- Event Analysis Code --------------//
1130
1131 ncl = 0;
1132 if( nw < 1 ) return ncl;
1133 if( nw < 2 ) { // if only one hit
1134 ncl = 1;
1135 lencl.push_back(1);
1136 return ncl;
1137 }
1138
1139 order_hyc( nw, ia, id ); // sort the addresses (ia) in increasing order
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; // no subclusters to be glued
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 } // end loop over previous subclusters
1201 } // end loop over all hits
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// order_hyc
1218//
1219//==========================================================
1220
1221void DCCALShower_factory::order_hyc( int nw, vector<int> &ia, vector<int> &id )
1222{
1223 // sort ia and id in order of increasing address
1224
1225 if( nw < 2 ) return; // only one hit
1226
1227 for( int k = 1; k < nw; k++ ) { // loop over hits
1228
1229 if( ia[k] <= ia[k-1] ) { // check if address is less than previous entry
1230 int iak = ia[k];
1231 int idk = id[k];
1232
1233 for( int ii = k-1; ii >= -1; ii = ii-1 ) { // loop over the previous entries
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 } // end loop over previous entries
1250 } // endif
1251 } // end loop over hits
1252
1253
1254 return;
1255}
1256
1257
1258
1259
1260
1261//==========================================================
1262//
1263// gams_hyc
1264//
1265//==========================================================
1266
1267void DCCALShower_factory::gams_hyc( int nadc, vector<int> &ia, vector<int> &id,
1268 int &nadcgam, vector<gamma_t> &gammas )
1269{
1270
1271 //-------------Local Declarations------------//
1272
1273 int ngam0; // number of gammas found before
1274 int niter; // max number of iterations (6)
1275 int npk; // number of peaks in the cluster
1276 int ipnpk[10]; // counter number with max energy in a peak
1277 int igmpk[2][10]; //
1278 int minpk; // min energy of a counter in a cluster
1279 int idelta; //
1280 int itype; // type of peak
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; // current value of chi2
1293 double chisq1; // value of chi2 for preliminary gammas separation
1294 double chisq2; // value of chi2 for final gammas separation
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]; // working array for resolved peaks
1317 vector<int> idp[13]; // energy of each cell of the island belonging to each peak
1318 vector<double> fwrk[13]; // working array for resolved peaks
1319
1320 // allocate memory for each vector:
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 // peaks search:
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 // gammas search for one peak:
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 { // cluster with more than one peak
1472
1473 //------------------------------------------
1474 /*
1475 First step - 1 gamma in each peak.
1476 Do a preliminary estimation of (E,x,y) of each peak, and split each peak into two hits
1477 only if it is badly needed (chi2 improvement is too high).
1478 If this split occurs, it is only for better (E,x,y) estimation, as it will be
1479 rejoined and reanalyzed in the second step.
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 } // end loop over peaks
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 } // end of iterations to separate peaks in a cluster
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 } // end loop over peaks
1668
1669
1670
1671 /*
1672 This is the second step: ( 1 or 2 gamma in each peak )
1673 (e,x,y) of hits were preliminarily estimated in the first step.
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 // part of gamma 'ig' energy belonging to cell 'i' from peak 'ipk':
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 } // end loop over hits
1711 } // end loop over peaks
1712
1713
1714
1715
1716 // Recover working array and renormalize total sum to the original cell energy:
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 //cout << "WARNING NEGATIVE CORR: ia = " << ia[ii] << "; id = " << id[ii] << endl;
1747 }
1748 } // end loop over peaks
1749
1750 } // end loop over hits
1751
1752
1753
1754 // erase the gammas found in the previous step:
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 } // end loop over peaks
1861
1862
1863
1864 } // end looping over multi-peak cluster
1865
1866
1867
1868
1869
1870 return;
1871}
1872
1873
1874
1875
1876
1877//==========================================================
1878//
1879// peak_type
1880//
1881//==========================================================
1882
1883int DCCALShower_factory::peak_type( int ix, int iy )
1884{
1885 /*
1886 itype = 2 : if the peak is in the most outer layer
1887 itype = 1 : if the peak is in the most inner layer
1888 itype = 0 : if the peak is anywhere else
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// gamma_hyc
1906//
1907//==========================================================
1908
1909void 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 //-------------Local Declarations------------//
1914
1915 int dof;
1916
1917 double dxy; // initial step for iteration
1918 double stepmin; // minimum step for iteration
1919 double stepx, stepy; // current steps
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 //-------------Event Analysis Code------------//
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 ); // calculate initial values of (E,x,y)
1945
1946 if( nadc <= 0 ) return;
1947
1948 chimem = chisq;
1949 chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x1, y1, chi0 ); // initial value of chi2
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 // start of iterations
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 // if steps at minimum, end iterations
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 // if chi2 at minimum, end iterations
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 ) { // if chi2 improved, then fix the improved values
2008 x1 = x0;
2009 y1 = y0;
2010 chisq = chi0/dof;
2011 }
2012
2013 // if chi2 is less than maximum allowed for one gamma in a peak, return.
2014 // otherwise, try separating the peak into two gammas:
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. ) { // if chi2 improved, decide if the separation
2022 // has physical meaning by calculating the
2023 // effective mass of the two gammas
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; // mass in MeV; 1270 = zccal
2029
2030 if( xm2 > xm2cut*xsize2.09 ) { // if the separation has physical meaning
2031 e1 = ee; // fix the parameters of the first gamma
2032 x1 = xx;
2033 y1 = yy;
2034 } else { // if no physical meaning e2=0
2035 e2 = 0.; // (second gamma is killed)
2036 chisq = chiold;
2037 }
2038 }
2039
2040 return;
2041}
2042
2043
2044
2045
2046
2047//==========================================================
2048//
2049// fill_zeros
2050//
2051//==========================================================
2052
2053void 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 ) { // fill left neib
2065 nneib = nneib+1;
2066 ian.push_back(iy + (ix-1)*100);
2067
2068 if( iy > 1 ) { // fill bottom left neib
2069 nneib = nneib+1;
2070 ian.push_back(iy-1 + (ix-1)*100);
2071 }
2072 if( iy < MROW12 ) { // fill top left neib
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 ) { // fill bottom right neib
2082 nneib = nneib+1;
2083 ian.push_back(iy-1 + (ix+1)*100);
2084 }
2085 if( iy < MROW12 ) { // fill top right neib
2086 nneib = nneib+1;
2087 ian.push_back(iy+1 + (ix+1)*100);
2088 }
2089 }
2090
2091 if( iy > 1 ) { // fill bottom neib
2092 nneib = nneib+1;
2093 ian.push_back(iy-1 + ix*100);
2094 }
2095 if( iy < MROW12 ) { // fill top neib
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// mom1_pht
2139//
2140//==========================================================
2141
2142void 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 //-------------Local Declarations------------//
2146
2147 int ix, iy;
2148 double dx, dy;
2149
2150 double a; //
2151 double corr; // correction to energy
2152
2153 //--------------Event Analysis Code-------------//
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// chisq1_hyc
2216//
2217//==========================================================
2218
2219void 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 //-------------Local Declarations------------//
2223
2224 int ix, iy;
2225 double dx, dy;
2226 double fa, fcell;
2227
2228
2229 //--------------Event Analysis Code-------------//
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 //if( SHOWER_DEBUG ) cout << "case 0 ch" << endl;
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 //if( SHOWER_DEBUG ) cout << "case 0 ch" << endl;
2263 }
2264 }
2265
2266
2267 return;
2268}
2269
2270
2271
2272
2273
2274//==========================================================
2275//
2276// sigma2, d2c, & cell_hyc
2277//
2278//==========================================================
2279
2280double 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
2301double 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 //if( i < 500. && j < 500. ) {
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
2329double 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 //std::cout << "i " << i << " j " << j << " wx " << wx << " wy " << wy << std::endl;
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// tgamma_hyc( int, vector<int>, vector<int>, int, vector<int>, double,
2362// double, double, double, double, double, double )
2363//
2364//==========================================================
2365
2366void 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 //-------------- Local Declarations -------------//
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 //-------------- Event Analysis Code -------------//
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 // choosing of the starting point
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) // e0 has been checked above
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 // start of iterations:
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++ ) { // chi2 calculation
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 ) { // new step if no improvement
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 { // calculate gradient
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 //a0 = a1 + a2;
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// mom2_pht( int, vector<int>, vector<int>, int, vector<int>,
2609// double, double, double, double, double, double )
2610//
2611//==========================================================
2612
2613void 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 //-------------- Local Declarations --------------//
2618
2619 int ix, iy;
2620 double dx, dy;
2621 double a;
2622 double corr;
2623
2624 //-------------- Event Analysis Code -------------//
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// c3to5_pht( double, double, double, double, double, double,
2694// double, double, double, double, double, double )
2695//
2696//==========================================================
2697
2698void 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// chisq2t_hyc( double, double, double, double, double, double,
2719// double, double, double, double )
2720//
2721//==========================================================
2722
2723void 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// ucopy1( vector<int>, vector<int>, int, int )
2759//
2760//==========================================================
2761
2762void 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// ucopy2( vector<int>, int, int, int )
2778//
2779//==========================================================
2780
2781void 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// ucopy3( vector<int>, vector<int>, int, int )
2803//
2804//==========================================================
2805
2806void 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