Bug Summary

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