Bug Summary

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