Bug Summary

File:libraries/CCAL/DCCALShower_factory.cc
Location:line 2457, 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 cout << "case 0 ch" << endl;
2247 //if( SHOWER_DEBUG ) cout << "case 0 ch" << endl;
2248 }
2249 }
2250
2251 for(int ii = 0; ii < nneib; ii++ ) {
2252 ix = iaz[ii]/100;
2253 iy = iaz[ii] - ix*100;
2254 dx = x1 - static_cast<double>(ix);
2255 dy = y1 - static_cast<double>(iy);
2256 if( e1 != 0. ) {
2257 if( fabs(dx) < 6.0 && fabs(dy) < 6.0 ) {
2258 fcell = cell_hyc( dx, dy );
2259 chisq = chisq + e1*fcell*fcell/sigma2(dx, dy, fcell, e1);
2260 }
2261 } else {
2262 chisq = chisq + id[ii]*id[ii]/9.;
2263 //if( SHOWER_DEBUG ) cout << "case 0 ch" << endl;
2264 }
2265 }
2266
2267
2268 return;
2269}
2270
2271
2272
2273
2274
2275//==========================================================
2276//
2277// sigma2, d2c, & cell_hyc
2278//
2279//==========================================================
2280
2281double DCCALShower_factory::sigma2( double dx, double dy, double fc, double e )
2282{
2283 double sig2;
2284 double alp = 0.816;
2285 double bet1 = 32.1;
2286 double bet2 = 1.72;
2287
2288 double r = dx*dx + dy*dy;
2289 if( r > 25. ) {
2290 sig2 = 100.;
2291 return sig2;
2292 }
2293
2294 sig2 = alp*fc + (bet1 + bet2*sqrt(e/100.))*d2c( dx, dy ) + 0.2/(e/100.);
2295 if( TEST_P1 ) sig2 = sig2/pow(0.0001*e, 0.166);
2296 sig2 *= 100.;
2297
2298 return sig2;
2299}
2300
2301
2302double DCCALShower_factory::d2c( double dx, double dy )
2303{
2304 int i, j;
2305 double ax, ay, wx, wy;
2306 double d2c;
2307
2308 ax = fabs(dx*100.);
2309 ay = fabs(dy*100.);
2310 i = int(ax);
2311 j = int(ay);
2312
2313 if( i < 499 && j < 499 && i >= 0 && j >= 0 ) {
2314 //if( i < 500. && j < 500. ) {
2315
2316 wx = ax-static_cast<double>(i);
2317 wy = ay-static_cast<double>(j);
2318
2319 d2c = ad2c[i][j] * (1.-wx) * (1.-wy) +
2320 ad2c[i+1][j] * wx * (1.-wy) +
2321 ad2c[i][j+1] * (1.-wx) * wy +
2322 ad2c[i+1][j+1] * wx * wy;
2323
2324 } else d2c = 1.;
2325
2326 return d2c;
2327}
2328
2329
2330double DCCALShower_factory::cell_hyc( double dx, double dy )
2331{
2332 int i, j;
2333 double ax, ay, wx, wy;
2334 double cell_hyc;
2335
2336 ax = fabs(dx*100.);
2337 ay = fabs(dy*100.);
2338 i = static_cast<int>(ax);
2339 j = static_cast<int>(ay);
2340
2341 if( i < 499 && j < 499 && i >= 0 && j >= 0 ) {
2342
2343 wx = ax-static_cast<double>(i);
2344 wy = ay-static_cast<double>(j);
2345 //std::cout << "i " << i << " j " << j << " wx " << wx << " wy " << wy << std::endl;
2346 cell_hyc = acell[i][j] * (1.-wx) * (1.-wy) +
2347 acell[i+1][j] * wx * (1.-wy) +
2348 acell[i][j+1] * (1.-wx) * wy +
2349 acell[i+1][j+1] * wx * wy;
2350
2351 } else cell_hyc = 0.;
2352
2353 return cell_hyc;
2354}
2355
2356
2357
2358
2359
2360//==========================================================
2361//
2362// tgamma_hyc( int, vector<int>, vector<int>, int, vector<int>, double,
2363// double, double, double, double, double, double )
2364//
2365//==========================================================
2366
2367void DCCALShower_factory::tgamma_hyc( int nadc, vector<int> ia, vector<int> id,
2368 int nzero, vector<int> iaz, double &chisq, double &e1, double &x1, double &y1,
2369 double &e2, double &x2, double &y2 )
2370{
2371 //-------------- Local Declarations -------------//
2372
2373 int ix, iy;
2374 int dof;
2375
2376 double dx, dy;
2377 double dxy;
2378 double dxc, dyc;
2379 double dx0, dy0;
2380 double dx1, dy1;
2381 double dx2, dy2;
2382 double u, r, rsq, rsq2;
2383 double epsc, eps0, eps1, eps2;
2384 double stepmin, epsmax;
2385 double delch;
2386 double step;
2387 double cosi, scal;
2388 double chisq2, chisqc;
2389 double dchi, dchida, dchi0;
2390 double a1, a2;
2391 double ex;
2392 double e1c, x1c, y1c;
2393 double e2c, x2c, y2c;
2394 double gr, gre, grx, gry;
2395 double grc;
2396 double grec, grxc, gryc;
2397 double gx1, gy1;
2398 double gx2, gy2;
2399 double e0, x0, y0;
2400 double xx, yy, yx;
2401
2402 double f1c, f2c, f1x, f2x, f1y, f2y;
2403 double chisqt, chisqt0, chisqt1, chisqt2;
2404 double chisqtx1, chisqtx2, chisqty1, chisqty2;
2405 double dchidax1, dchidax2, dchiday1, dchiday2;
2406
2407 stepmin = 0.5;
2408 epsmax = 0.9999;
2409 delch = 10.;
2410
2411 //-------------- Event Analysis Code -------------//
2412
2413 yx = 0.;
2414 grx = 0.; gry = 0.;
2415 gre = 0.; gr = 0.;
2416 dx0 = 0.; dy0 = 0.;
2417 eps0 = 0.;
2418
2419 mom2_pht( nadc, ia, id, nzero, iaz, e0, x0, y0, xx, yy, yx );
2420
2421 e2 = 0.;
2422 x2 = 0.;
2423 y2 = 0.;
2424
2425 if( nadc <= 0 ) return;
2426
2427 // choosing of the starting point
2428
2429 dxy = xx - yy;
2430 rsq2 = dxy*dxy + 4.*yx*yx;
2431 if( rsq2 < 1.e-20 ) rsq2 = 1.e-20;
2432 rsq = sqrt(rsq2);
2433 dxc = -sqrt((rsq+dxy)*2.);
2434 dyc = sqrt((rsq-dxy)*2.);
2435 if( yx >= 0. ) dyc = -dyc;
2436 r = sqrt(dxc*dxc + dyc*dyc);
2437 epsc = 0.;
2438 for( int ii = 0; ii < nadc; ii++ ) {
2439 ix = ia[ii]/100;
2440 iy = ia[ii] - ix*100;
2441 dx = static_cast<double>(ix) - x0;
2442 dy = static_cast<double>(iy) - y0;
2443 u = dx*dxc/r + dy*dyc/r;
2444 epsc = epsc - 0.01*id[ii]*u*fabs(u);
2445 }
2446 epsc = epsc/(0.01*e0*rsq);
2447 if( epsc > 0.8 ) epsc = 0.8;
2448 if( epsc < -0.8 ) epsc = -0.8;
2449 dxc = dxc/sqrt(1.-(epsc*epsc));
2450 dyc = dyc/sqrt(1.-(epsc*epsc));
2451
2452
2453
2454 // start of iterations:
2455
2456 step = 0.1;
2457 cosi = 0.0;
Value stored to 'cosi' is never read
2458 chisq2 = 1.e35;
2459
2460 for( int iter = 0; iter < 100; iter++ ) {
2461
2462 c3to5_pht( e0, x0, y0, epsc, dxc, dyc, e1c, x1c, y1c, e2c, x2c, y2c );
2463
2464 eps1 = (1. + epsc)/2.;
2465 eps2 = (1. - epsc)/2.;
2466 chisqc = 0.;
2467
2468 for( int ii = 0; ii < nadc+nzero; ii++ ) { // chi2 calculation
2469
2470 if( ii < nadc ) {
2471 ex = static_cast<double>(id[ii]);
2472 ix = ia[ii]/100;
2473 iy = ia[ii] - ix*100;
2474 } else {
2475 ex = 0.;
2476 ix = iaz[ii-nadc]/100;
2477 iy = iaz[ii-nadc] - ix*100;
2478 }
2479
2480 dx1 = x1c - static_cast<double>(ix);
2481 dy1 = y1c - static_cast<double>(iy);
2482 dx2 = x2c - static_cast<double>(ix);
2483 dy2 = y2c - static_cast<double>(iy);
2484
2485 f1c = cell_hyc( dx1, dy1 );
2486 f2c = cell_hyc( dx2, dy2 );
2487
2488 chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt );
2489 chisqc += chisqt;
2490
2491 }
2492
2493 if( chisqc >= chisq2 ) { // new step if no improvement
2494
2495 if( iter > 0 ) {
2496 dchi = chisqc - chisq2;
2497 dchi0 = gr*step;
2498 step = 0.5*step/sqrt(1. + dchi/dchi0);
2499 }
2500 step = 0.5*step;
2501
2502 } else { // calculate gradient
2503
2504 grec = 0.;
2505 grxc = 0.;
2506 gryc = 0.;
2507
2508 for( int ii = 0; ii < nadc+nzero; ii++ ) {
2509
2510 if( ii < nadc ) {
2511 ex = static_cast<double>(id[ii]);
2512 ix = ia[ii]/100;
2513 iy = ia[ii] - ix*100;
2514 } else {
2515 ex = 0.;
2516 ix = iaz[ii-nadc]/100;
2517 iy = iaz[ii-nadc] - ix*100;
2518 }
2519
2520 dx1 = x1c - static_cast<double>(ix);
2521 dy1 = y1c - static_cast<double>(iy);
2522 dx2 = x2c - static_cast<double>(ix);
2523 dy2 = y2c - static_cast<double>(iy);
2524
2525 f1c = cell_hyc( dx1, dy1 );
2526 f2c = cell_hyc( dx2, dy2 );
2527
2528 a1 = e1c*f1c;
2529 a2 = e2c*f2c;
2530 //a0 = a1 + a2;
2531
2532 chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt0 );
2533 chisq2t_hyc( ex, e1c+1., dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt1 );
2534 chisq2t_hyc( ex, e1c, dx1, dy1, e2c+1., dx2, dy2, f1c, f2c, chisqt2 );
2535
2536 f1x = cell_hyc( dx1+0.05, dy1 );
2537 f2x = cell_hyc( dx2+0.05, dy2 );
2538 f1y = cell_hyc( dx1, dy1+0.05 );
2539 f2y = cell_hyc( dx2, dy2+0.05 );
2540
2541 chisq2t_hyc( ex, e1c, dx1+0.05, dy1, e2c, dx2, dy2, f1x, f2c, chisqtx1 );
2542 chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2+0.05, dy2, f1c, f2x, chisqtx2 );
2543 chisq2t_hyc( ex, e1c, dx1, dy1+0.05, e2c, dx2, dy2, f1y, f2c, chisqty1 );
2544 chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2+0.05, f1c, f2y, chisqty2 );
2545
2546 dchidax1 = 20.*(chisqtx1 - chisqt0);
2547 dchidax2 = 20.*(chisqtx2 - chisqt0);
2548 dchiday1 = 20.*(chisqty1 - chisqt0);
2549 dchiday2 = 20.*(chisqty2 - chisqt0);
2550 dchida = 0.5*(chisqt1 + chisqt2 - chisqt0);
2551
2552 gx1 = (e1c*f1x - a1)*dchidax1;
2553 gx2 = (e2c*f2x - a2)*dchidax2;
2554 gy1 = (e1c*f1y - a1)*dchiday1;
2555 gy2 = (e2c*f2y - a2)*dchiday2;
2556
2557 grec += dchida*(f1c-f2c)*e0 - ( (gx1+gx2)*dxc + (gy1+gy2)*dyc );
2558 grxc += gx1*eps2 - gx2*eps1;
2559 gryc += gy1*eps2 - gy2*eps1;
2560
2561 }
2562
2563 grc = sqrt( grec*grec + grxc*grxc + gryc*gryc );
2564 if( grc < 1.e-6 ) grc = 1.e-6;
2565 if( iter > 0 ) {
2566 cosi = ( gre*grec + grx*grxc + gry*gryc )/(gr*grc);
2567 scal = fabs((gr/grc) - cosi);
2568 if( scal < 0.1 ) scal = 0.1;
2569 step = step/scal;
2570 }
2571 chisq2 = chisqc;
2572 eps0 = epsc;
2573 dx0 = dxc;
2574 dy0 = dyc;
2575 gre = grec;
2576 grx = grxc;
2577 gry = gryc;
2578 gr = grc;
2579 }
2580 epsc = eps0 - step*gre/gr;
2581 while( fabs(epsc) > epsmax ) {
2582 step = step/2.;
2583 epsc = eps0 - step*gre/gr;
2584 }
2585 dxc = dx0 - step*grx/gr;
2586 dyc = dy0 - step*gry/gr;
2587 if( step*gr < stepmin ) break;
2588 }
2589
2590 if( (chisq*(nadc+nzero-2) - chisq2) < delch ) return;
2591 dof = nzero+nadc-5;
2592 if( dof < 1 ) dof = 1;
2593 chisq = chisq2/dof;
2594
2595 c3to5_pht( e0, x0, y0, eps0, dx0, dy0, e1, x1, y1, e2, x2, y2 );
2596
2597 return;
2598}
2599
2600
2601
2602
2603
2604//==========================================================
2605//
2606// mom2_pht( int, vector<int>, vector<int>, int, vector<int>,
2607// double, double, double, double, double, double )
2608//
2609//==========================================================
2610
2611void DCCALShower_factory::mom2_pht( int nadc, vector<int> ia, vector<int> id,
2612 int nzero, vector<int> iaz, double &a0, double &x0, double &y0,
2613 double &xx, double &yy, double &yx)
2614{
2615 //-------------- Local Declarations --------------//
2616
2617 int ix, iy;
2618 double dx, dy;
2619 double a;
2620 double corr;
2621
2622 //-------------- Event Analysis Code -------------//
2623
2624 a0 = 0.;
2625 x0 = 0.;
2626 y0 = 0.;
2627 xx = 0.;
2628 yy = 0.;
2629 yx = 0.;
2630
2631 if( nadc <= 0 ) return;
2632 for( int ii = 0; ii < nadc; ii++ ) {
2633 a = static_cast<double>(id[ii]);
2634 ix = ia[ii]/100;
2635 iy = ia[ii] - ix*100;
2636 a0 = a0 + a;
2637 x0 = x0 + a*static_cast<double>(ix);
2638 y0 = y0 + a*static_cast<double>(iy);
2639 }
2640 if( a0 <= 0. ) return;
2641 x0 = x0/a0;
2642 y0 = y0/a0;
2643
2644 for( int ii = 0; ii < nadc; ii++ ) {
2645 a = static_cast<double>(id[ii])/a0;
2646 ix = ia[ii]/100;
2647 iy = ia[ii] - ix*100;
2648 dx = static_cast<double>(ix) - x0;
2649 dy = static_cast<double>(iy) - y0;
2650 xx = xx + a*dx*dx;
2651 yy = yy + a*dy*dy;
2652 yx = yx + a*dx*dy;
2653 }
2654
2655 corr = 0.;
2656 for( int ii = 0; ii < nadc; ii++ ) {
2657 ix = ia[ii]/100;
2658 iy = ia[ii] - ix*100;
2659 dx = static_cast<double>(ix) - x0;
2660 dy = static_cast<double>(iy) - y0;
2661 corr = corr + cell_hyc( dx, dy );
2662 }
2663 if( COUNT_ZERO1 ) {
2664 for( int ii = 0; ii < nzero; ii++ ) {
2665 ix = iaz[ii]/100;
2666 iy = iaz[ii] - ix*100;
2667 dx = static_cast<double>(ix) - x0;
2668 dy = static_cast<double>(iy) - y0;
2669 corr = corr + cell_hyc( dx, dy );
2670 }
2671 }
2672
2673 corr = corr/1.006;
2674 if( corr < 0.8 ) {
2675 corr = 0.8;
2676 } else if( corr > 1. ) {
2677 corr = 1.0;
2678 }
2679 a0 = a0/corr;
2680
2681
2682 return;
2683}
2684
2685
2686
2687
2688
2689//==========================================================
2690//
2691// c3to5_pht( double, double, double, double, double, double,
2692// double, double, double, double, double, double )
2693//
2694//==========================================================
2695
2696void DCCALShower_factory::c3to5_pht( double e0, double x0, double y0, double eps,
2697 double dx, double dy, double &e1, double &x1, double &y1, double &e2,
2698 double &x2, double &y2 )
2699{
2700 e1 = e0*(1.+eps)/2.;
2701 e2 = e0 - e1;
2702 x1 = x0 + 0.5*dx*(1.-eps);
2703 y1 = y0 + 0.5*dy*(1.-eps);
2704 x2 = x0 - 0.5*dx*(1.+eps);
2705 y2 = y0 - 0.5*dy*(1.+eps);
2706
2707 return;
2708}
2709
2710
2711
2712
2713
2714//==========================================================
2715//
2716// chisq2t_hyc( double, double, double, double, double, double,
2717// double, double, double, double )
2718//
2719//==========================================================
2720
2721void DCCALShower_factory::chisq2t_hyc( double ecell, double e1, double dx1, double dy1,
2722 double e2, double dx2, double dy2, double f1, double f2, double &chisqt )
2723{
2724 double s;
2725 double p1, p2;
2726
2727 if( TEST_P1 ) {
2728 p1 = pow(0.0001*e1, 0.166);
2729 p2 = pow(0.0001*e2, 0.166);
2730 } else {
2731 p1 = 1.;
2732 p2 = 1.;
2733 }
2734
2735 if( e1 != 0. && e2 != 0. )
2736 s = e1*sigma2(dx1, dy1, f1, e1)/p1 + e2*sigma2(dx2, dy2, f2, e2)/p2;
2737 else if( e1 == 0. && e2 == 0. )
2738 s = 90000.;
2739 else if( e1 == 0. )
2740 s = e2*sigma2(dx2, dy2, f2, e2)/p2;
2741 else
2742 s = e1*sigma2(dx1, dy1, f1, e1)/p1;
2743
2744 chisqt = pow((e1*f1 + e2*f2 - ecell), 2.)/s;
2745
2746 return;
2747}
2748
2749
2750
2751
2752
2753
2754//==========================================================
2755//
2756// ucopy1( vector<int>, vector<int>, int, int )
2757//
2758//==========================================================
2759
2760void DCCALShower_factory::ucopy1( vector<int> &ia, vector<int> &iwork, int start, int length )
2761{
2762
2763 for( int ii = 0; ii < length; ii++ ) {
2764 iwork.push_back( ia[start+ii] );
2765 }
2766
2767
2768 return;
2769}
2770
2771
2772
2773//==========================================================
2774//
2775// ucopy2( vector<int>, int, int, int )
2776//
2777//==========================================================
2778
2779void DCCALShower_factory::ucopy2( vector<int> &ia, int start1, int start2, int length )
2780{
2781
2782 vector<int> work;
2783 for( int ii = 0; ii < length; ii++ ) {
2784 work.push_back( ia[start1+ii] );
2785 }
2786
2787 for( int ii = 0; ii < length; ii++ ) {
2788 ia[start2+ii] = work[ii];
2789 }
2790
2791
2792 return;
2793}
2794
2795
2796
2797
2798//==========================================================
2799//
2800// ucopy3( vector<int>, vector<int>, int, int )
2801//
2802//==========================================================
2803
2804void DCCALShower_factory::ucopy3( vector<int> &iwork, vector<int> &ia, int start, int length )
2805{
2806
2807 for( int ii = 0; ii < length; ii++ ) {
2808 ia[start+ii] = iwork[ii];
2809 }
2810
2811
2812 return;
2813}
2814
2815
2816
2817
2818
2819
2820
2821
2822
2823
2824
2825
2826
2827
2828
2829
2830
2831
2832
2833
2834
2835