AMSDST
myEcalShowerPlus.cxx
Go to the documentation of this file.
1 // Authors: M.Duranti - INFN di Perugia
2 #include "myEcalShowerPlus.h"
3 #include "debug.h"
4 #include "TClass.h"
5 #include "myEvent.h"
6 
7 using namespace std;
8 
9 //--------------------------------------------------------------------
10 ClassImp(myEcalShowerPlus);
11 //--------------------------------------------------------------------
12 
13 myEcalShowerPlus::myEcalShowerPlus(){
14 #ifdef PDEBUG
15  printf("In myEcalShowerPlus::myEcalShowerPlus\n");
16 #endif
17  PRINTDEBUG;
18  init();
19  PRINTDEBUG;
20 }
21 
22 myEcalShowerPlus::~myEcalShowerPlus(){
23 #ifdef PDEBUG
24  printf("In myEcalShowerPlus::~myEcalShowerPlus\n");
25 #endif
26  PRINTDEBUG;
27  Clear();
28  PRINTDEBUG;
29 }
30 
31 float myEcalShowerPlus::EdepUpToLayer(int ilay){
32 
33  if (ilay>17) printf("EdepUpToLayer(%d): ECAL has only 18 layers (0-17)...\n", ilay);
34 
35  float _EdepUpToLayer=0;
36 
37  for (int ii=0; ii<=ilay; ii++) {
38  _EdepUpToLayer+=Elayer[ii];
39  }
40 
41  return _EdepUpToLayer;
42 }
43 
44 float myEcalShowerPlus::EdepUpToSuperLayer(int islay){
45  return EdepUpToLayer((islay*2)+1);
46 }
47 
48 void myEcalShowerPlus::Clear(Option_t* option){
49 #ifdef PDEBUG
50  printf("In myEcalShowerPlus::Clear\n");
51 #endif
52  PRINTDEBUG;
53  _ecalhitmap.clear();
54  PRINTDEBUG;
55  fill_n(Elayer, 18, 0.);
56  fill_n(Elayer_corr, 18, 0.);
57  PRINTDEBUG;
58  fill_n(Energy3C, 3, 0.);
59  PRINTDEBUG;
60  //memset(EcalAxis_chi2_2, 0, 18*sizeof(EcalAxis_chi2_2[0]));
61  //memset(EcalAxis_chi2_8, 0, 18*sizeof(EcalAxis_chi2_2[0]));
62  fill_n(EcalAxis_chi2_2, 18, -99.);
63  fill_n(EcalAxis_chi2_8, 18, -99.);
64  PRINTDEBUG;
65  BDT = -9999999;
66  BDTs = -9999999;
67  BDTE = -9999999;
68  BDTEs = -9999999;
69  BDT_v4 = -9999999;
70  BDTChi2 = -9999999;
71  BDTChi2s = -9999999;
72  BDTChi2E = -9999999;
73  BDTChi2Es = -9999999;
74  BDTChi2_v2 = -9999999;
75  PRINTDEBUG;
76  ESEV2 = -9999999;
77  ESEV3 = -9999999;
78  ESEv3Efficiency = -9999999;
79  ESEv3Rejection = -9999999;
80  PRINTDEBUG;
81  return;
82 }
83 
84 void myEcalShowerPlus::init(){
85 #ifdef PDEBUG
86  printf("In myEcalShowerPlus::init\n");
87 #endif
88  PRINTDEBUG;
89  Clear();
90  PRINTDEBUG;
91 }
92 
93 float myEcalShowerPlus::GetEcalBDT(unsigned int iBDTVERSION, int TMVAClassifier, int EnergyFlag){
94 
95  float _bdt=-9999999;
96 
97  if (iBDTVERSION==5) {
98  if (TMVAClassifier==0) {
99  if (EnergyFlag==0) return BDT;
100  else if (EnergyFlag==1) return BDTE;
101  }
102  else if (TMVAClassifier==1) {
103  if (EnergyFlag==0) return BDTs;
104  else if (EnergyFlag==1) return BDTEs;
105  }
106  }
107  else if (iBDTVERSION==4) {
108  if (TMVAClassifier==0 && EnergyFlag==0) return BDT_v4;
109  }
110 
111  printf("BDT(%u, %d, %d) not stored, sorry...\n", iBDTVERSION, TMVAClassifier, EnergyFlag);
112  return _bdt;
113 }
114 
115 float myEcalShowerPlus::GetEcalBDTCHI2(unsigned int iBDTVERSION, int TMVAClassifier, int EnergyFlag){
116 
117  float _bdt=-9999999;
118 
119  if (iBDTVERSION==3) {
120  if (TMVAClassifier==0) {
121  if (EnergyFlag==0) return BDTChi2;
122  else if (EnergyFlag==1) return BDTChi2E;
123  }
124  else if (TMVAClassifier==1) {
125  if (EnergyFlag==0) return BDTChi2s;
126  else if (EnergyFlag==1) return BDTChi2Es;
127  }
128  }
129  else if (iBDTVERSION==2) {
130  if (TMVAClassifier==0 && EnergyFlag==0) return BDTChi2_v2;
131  }
132 
133  printf("BDTCHI2(%u, %d, %d) not stored, sorry...\n", iBDTVERSION, TMVAClassifier, EnergyFlag);
134  return _bdt;
135 }
136 
137 float myEcalShowerPlus::GetEcalHitADC(short int Cell, short int Plane, int component){
138  static std::map<short int, ecalhit>::iterator it;
139  short int ittp = ComputeIndex(Cell, Plane);
140  it = _ecalhitmap.find(ittp);
141  if (it!=_ecalhitmap.end()) return (float)(it->second.ADC[component]);
142  return -1;
143 }
144 
145 float myEcalShowerPlus::GetEcalHitEdep(short int Cell, short int Plane){
146  static std::map<short int, ecalhit>::iterator it;
147  short int ittp = ComputeIndex(Cell, Plane);
148  it = _ecalhitmap.find(ittp);
149  if (it!=_ecalhitmap.end()) return (float)(it->second.Edep);
150  return -1e-9;
151 }
152 
153 unsigned int myEcalShowerPlus::GetEcalHitStatus(short int Cell, short int Plane){
154  static std::map<short int, ecalhit>::iterator it;
155  short int ittp = ComputeIndex(Cell, Plane);
156  it = _ecalhitmap.find(ittp);
157  if (it!=_ecalhitmap.end()) return (unsigned int)(it->second.Status);
158  return 0;
159 }
160 
161 #if 0
162 void myEcalShowerPlus::sync()
163 {
164  exportCell( edep, EDEP );
165  exportCell( adc_highgain, ADC_HIGHGAIN );
166  exportCell( adc_lowgain, ADC_LOWGAIN );
167  exportCell( adc_dynodeC, ADC_DYNODE );
168  exportPMTADCDynode( adc_dynode );
169  exportCellStatus( cell_status );
170  /*
171  for( int i=0; i<18; ++i )
172  for( int j=0; j<72; ++j )
173  {
174  edep[i][j] = GetEcalHitEdep(j, i);
175  adc_highgain[i][j] = GetEcalHitADC(j, i, 0);
176  adc_lowgain[i][j] = GetEcalHitADC(j, i, 1);
177  adc_dynodeC[i][j] = GetEcalHitADC(j, i, 2);
178  adc_dynode[i/2][j/2] = adc_dynodeC[i][j];
179  }
180  */
181 }
182 #endif
183 
184 void myEcalShowerPlus::exportCell(float cell_arr[18][72], CellType t)
185 {
186  switch(t){
187  case EDEP: for(int i=0; i<18; ++i) for(int j=0; j<72; ++j) cell_arr[i][j] = GetEcalHitEdep(j, i); break;
188  case ADC_HIGHGAIN: for(int i=0; i<18; ++i) for(int j=0; j<72; ++j) cell_arr[i][j] = GetEcalHitADC(j, i, 0);break;
189  case ADC_LOWGAIN: for(int i=0; i<18; ++i) for(int j=0; j<72; ++j) cell_arr[i][j] = GetEcalHitADC(j, i, 1);break;
190  case ADC_DYNODE: for(int i=0; i<18; ++i) for(int j=0; j<72; ++j) cell_arr[i][j] = GetEcalHitADC(j, i, 2);break;
191  }
192 }
193 void myEcalShowerPlus::exportPMTADCDynode(float pmt_dy[9][36])
194 {
195  for(int i=0; i<9; ++i) for(int j=0; j<36; ++j) pmt_dy[i][j] = GetEcalHitADC(j*2, i*2, 2);
196 }
197 void myEcalShowerPlus::exportCellStatus(UInt_t cell_st[18][72])
198 {
199  for(int i=0; i<18; ++i) for(int j=0; j<72; ++j) cell_st[i][j] = GetEcalHitStatus(j, i);
200 }
201 
202 float myEcalShowerPlus::RecomputeBDT() {
203 
204  // La lista e' la seguente.
205  // 1. per ogni layer: LayerEneFrac[18], LayerS1S3[18], LayerS3Frac[18], LayerSigma[18], LayerChi2[18]
206  // 2. variabili generali: ShowerMean, F2SLEneDep, L2LFrac, R3cmFrac, S3totx, S3toty, NecalHits
207 
208  const unsigned int nLAYERs = 18;
209  const unsigned int nCELLs = 72;
210  const Float_t ecalZEntry = -142.792;
211  const Float_t ecalZExit = -158.457;
212  const Float_t EneDepThreshold = 2.;//threshold on single cell in MeV (1MeV~2ADC)
213 
214  float MapEneDep[nLAYERs][nCELLs]; // Energy deposit in every cell of ECAL [GeV]
215  float ClusterEnergy[nLAYERs]; // Lateral leak corrected energy deposit [GeV]
216  float LayerEneDep[nLAYERs]; // Energy deposit [GeV] in each layer (sum of every cell of each layer)
217 
218  float LayerEnergy = 0.; // sum_i ClusterEnergy[i]
219 
220  float EneDep = 0.;
221  float EneDepXY[2]; // 0 = x, 1 = y
222 
223  float ClusterEnergyXY[2]; // 0 = x, 1 = y
224 
225  float LayerMean[nLAYERs]; // Mean [cell] for each layer: (sum_j j*MapEneDep[i][j])/(sum_ij MapEneDep[i][j])
226  float LayerSigma[nLAYERs];
227  float LayerEneFrac[nLAYERs];
228  float LayerS1S3[nLAYERs];
229  float LayerS3Frac[nLAYERs];
230 
231  float ShowerMean = 0.;
232  float ShowerSigma = 0.;
233  float L2LFrac = 0.;
234  float S3totx = 0.;
235  float S3toty = 0.;
236  float R3cmFrac = Energy3C[0];
237  float ShowerFootprintX = 0.;
238  float ShowerFootprintY = 0.;
239  float NEcalHits = 0.;
240 
241  float ShowerMeanXY[2]; // 0 = x, 1 = y
242 
243  float sigmaXY[2]; // 0 = x, 1 = y
244  float sigmaXYZ[2]; // 0 = x, 1 = y
245  float sigmaZ[2]; // 0 = x, 1 = y
246 
247  float F2SLEneDep = 0.;
248 
249  // Initialize array variables
250  for (unsigned int iproj = 0; iproj < 2; ++iproj) {
251  EneDepXY[iproj] = 0.;
252  ClusterEnergyXY[iproj] = 0.;
253  ShowerMeanXY[iproj] = 0.;
254  sigmaXY[iproj] = 0.;
255  sigmaXYZ[iproj] = 0.;
256  sigmaZ[iproj] = 0.;
257  }
258 
259  for (unsigned int ilayer = 0; ilayer < nLAYERs; ++ilayer) {
260  ClusterEnergy[ilayer] = 0.;
261  LayerEneDep[ilayer] = 0.;
262  LayerEneFrac[ilayer] = 0.;
263 
264  LayerMean[ilayer] = 0.;
265  LayerSigma[ilayer] = 0.;
266 
267  LayerS1S3[ilayer] = 1.;
268  LayerS3Frac[ilayer] = 1.;
269 
270  for (unsigned int icell = 0; icell < nCELLs; ++icell)
271  MapEneDep[ilayer][icell] = 0.;
272  }
273 
274  //---------------------------------------------------------------
275 
276  std::map<short int, ecalhit>::iterator it;
277  for (it=_ecalhitmap.begin(); it!=_ecalhitmap.end(); ++it) {
278  ecalhit hit = it->second;
279  short int index = it->first;
280  ClusterEnergy[GetPlane(index)] += hit.Edep;
281  MapEneDep[GetPlane(index)][GetCell(index)] = hit.Edep;
282  }
283 
284  // Compute mean, s1, s3 and energies
285  for (unsigned int ilayer = 0; ilayer < nLAYERs; ++ilayer)
286  {
287  if (ClusterEnergy[ilayer] == 0.) continue;
288 
289  UChar_t proj = !((ilayer/2) % 2);
290 
291  ClusterEnergyXY[proj] += ClusterEnergy[ilayer];
292  LayerEnergy += ClusterEnergy[ilayer];
293 
294  ShowerMean += (ilayer + 1)*ClusterEnergy[ilayer];
295 
296  int imaxcell = 0;
297  float maxcellene = 0.;
298 
299  for (unsigned int icell = 0; icell < nCELLs; ++icell)
300  {
301  if (MapEneDep[ilayer][icell] <= EneDepThreshold ) continue;
302 
303  LayerEneDep[ilayer] += MapEneDep[ilayer][icell];
304 
305  ++NEcalHits;
306 
307  if (MapEneDep[ilayer][icell] >= maxcellene)
308  {
309  maxcellene = MapEneDep[ilayer][icell];
310  imaxcell = icell;
311  }
312 
313  LayerMean[ilayer] += (icell + 1)*MapEneDep[ilayer][icell];
314  ShowerMeanXY[proj] += (icell + 1)*MapEneDep[ilayer][icell];
315  }
316 
317  //
318  // Now apply corrections for lateral leakage -- ADDED 27may2013 MI
319  //
320  float LatLeak[10]={0.};
321  //case 1: cell at left border -> simmetrize shower
322  if (imaxcell==0)
323  {
324  int iLatLeak=1;
325  while(MapEneDep[ilayer][imaxcell+iLatLeak]>0.&&iLatLeak<10)
326  {
327  LatLeak[iLatLeak-1] = MapEneDep[ilayer][imaxcell+iLatLeak];
328  iLatLeak++;
329  }
330  }
331  //case 2: cell close to left border within 10 cells
332  //-> simmetrize shower using ratio of adjacent cells
333  else if (imaxcell<9)
334  {
335  //find ratio
336  float LatRatio = MapEneDep[ilayer][imaxcell-1]/MapEneDep[ilayer][imaxcell+1] ;
337  int iLatLeak=1;
338  while(MapEneDep[ilayer][imaxcell+iLatLeak]>0.&&iLatLeak<10)
339  {
340  if ((int)imaxcell-iLatLeak < 0) LatLeak[iLatLeak-imaxcell-1] = LatRatio*MapEneDep[ilayer][imaxcell+iLatLeak];
341  iLatLeak++;
342  }
343  }
344  //case 3: cell at right border -> simmetrize shower
345  if (imaxcell==71)
346  {
347  int iLatLeak=1;
348  while(MapEneDep[ilayer][imaxcell-iLatLeak]>0.&&iLatLeak<10)
349  {
350  LatLeak[iLatLeak-1] = MapEneDep[ilayer][imaxcell-iLatLeak];
351  iLatLeak++;
352  }
353  }
354  //case 4: cell close to left border within 10 cells
355  //-> simmetrize shower using ratio of adjacent cells
356  else if (imaxcell>62)
357  {
358  //find ratio
359  float LatRatio = MapEneDep[ilayer][imaxcell+1]/MapEneDep[ilayer][imaxcell-1] ;
360  int iLatLeak=1;
361  while(MapEneDep[ilayer][imaxcell-iLatLeak]>0.&&iLatLeak<10)
362  {
363  if ((int)imaxcell-iLatLeak < 0) LatLeak[iLatLeak-imaxcell-1] = LatRatio*MapEneDep[ilayer][imaxcell-iLatLeak];
364  iLatLeak++;
365  }
366  }
367  //
368  for(int i=0;i<10;i++) LayerEneDep[ilayer] += LatLeak[i];
369  //End of correction for lateral leakage
370  //
371 
372  EneDep += LayerEneDep[ilayer];
373  EneDepXY[proj] += LayerEneDep[ilayer];
374 
375  float S1 = MapEneDep[ilayer][imaxcell];
376  float S3 = S1;
377  if (imaxcell==0 || imaxcell==71) S3 += LatLeak[0];
378  if (imaxcell > 0) S3 += MapEneDep[ilayer][imaxcell - 1];
379  if (imaxcell < (int)(nCELLs - 1)) S3 += MapEneDep[ilayer][imaxcell + 1];
380 
381  if (!proj) S3totx += S3;
382  else S3toty += S3;
383 
384  if (S1*S3 > 0.) LayerS1S3[ilayer] = S1/S3;
385 
386  if (LayerEneDep[ilayer] > 0.)
387  {
388  LayerS3Frac[ilayer] = S3/LayerEneDep[ilayer];
389  LayerMean[ilayer] = LayerMean[ilayer]/LayerEneDep[ilayer];
390  }
391  else
392  {
393  LayerMean[ilayer] = -1.;
394  }
395  } //end of loop on layers
396 
397  if (EneDep <= 0. || LayerEnergy <= 0. || EneDepXY[0] <= 0. || EneDepXY[1] <= 0.) return -0.999;
398 
399  ShowerMeanXY[0] /= EneDepXY[0];
400  ShowerMeanXY[1] /= EneDepXY[1];
401  ShowerMean /= LayerEnergy;
402 
403  F2SLEneDep = (LayerEneDep[0] + LayerEneDep[1] + LayerEneDep[2] + LayerEneDep[3])/1000.;
404  L2LFrac = (LayerEneDep[16] + LayerEneDep[17])/EneDep;
405 
406  S3totx /= EneDepXY[0];
407  S3toty /= EneDepXY[1];
408 
409  for (unsigned int ilayer = 0; ilayer < nLAYERs; ++ilayer)
410  {
411  if (ClusterEnergy[ilayer] == 0.) continue;
412 
413  UChar_t proj = !((ilayer/2) % 2);
414 
415  LayerEneFrac[ilayer] = LayerEneDep[ilayer]/EneDep;
416 
417  ShowerSigma += TMath::Power(ilayer + 1 - ShowerMean, 2)*ClusterEnergy[ilayer];
418 
419  for (unsigned int icell = 0; icell < nCELLs; ++icell)
420  {
421  if (MapEneDep[ilayer][icell] <= EneDepThreshold) continue;
422 
423  LayerSigma[ilayer] += TMath::Power(icell + 1 - LayerMean[ilayer], 2)*MapEneDep[ilayer][icell];
424 
425  sigmaXY[proj] += TMath::Power(icell + 1 - ShowerMeanXY[proj], 2)*MapEneDep[ilayer][icell];
426  sigmaXYZ[proj] += (icell + 1 - ShowerMeanXY[proj])*(ilayer + 1 - ShowerMean)*MapEneDep[ilayer][icell];
427  sigmaZ[proj] += TMath::Power(ilayer + 1 - ShowerMean, 2)*MapEneDep[ilayer][icell];
428  }
429 
430  if (LayerEneDep[ilayer] > 0.)
431  LayerSigma[ilayer] = TMath::Sqrt(LayerSigma[ilayer]/LayerEneDep[ilayer]);
432  else LayerSigma[ilayer] = -1.;
433  }
434 
435  ShowerSigma = TMath::Sqrt(ShowerSigma/LayerEnergy); // ShowerSigma from 0 to 17
436  ShowerMean = ShowerMean - 1.;
437 
438  for (unsigned int iproj = 0; iproj < 2; ++iproj)
439  {
440  sigmaXY[iproj] /= EneDepXY[iproj];
441  sigmaXYZ[iproj] /= EneDepXY[iproj];
442  sigmaZ[iproj] /= EneDepXY[iproj];
443  }
444 
445  ShowerFootprintX = TMath::Sqrt(TMath::Abs(sigmaXY[0]*sigmaZ[0] - TMath::Power(sigmaXYZ[0], 2)));
446  ShowerFootprintY = TMath::Sqrt(TMath::Abs(sigmaXY[1]*sigmaZ[1] - TMath::Power(sigmaXYZ[1], 2)));
447 
448  return -999;
449 }
450 
451 //---------------------------------------------------------------------
452 
453 #ifdef _WITHGBATCH_
454 
455 //--------------------------------------------------------------------
456 ClassImp(myEcalShowerPlusFiller);
457 //--------------------------------------------------------------------
458 
459 myEcalShowerPlusFiller::myEcalShowerPlusFiller(){
460 #ifdef PDEBUG
461  printf("In myEcalShowerPlusFiller::myEcalShowerPlusFiller\n");
462 #endif
463  PRINTDEBUG;
464  init();
465  PRINTDEBUG;
466 }
467 
468 myEcalShowerPlusFiller::~myEcalShowerPlusFiller(){
469 #ifdef PDEBUG
470  printf("In myEcalShowerPlusFiller::~myEcalShowerPlusFiller\n");
471 #endif
472  PRINTDEBUG;
473  Clear();
474  PRINTDEBUG;
475 }
476 
477 void myEcalShowerPlusFiller::Clear(Option_t* option){
478 #ifdef PDEBUG
479  printf("In myEcalShowerPlusFiller::Clear\n");
480 #endif
481  PRINTDEBUG;
482  shower=0;
483  PRINTDEBUG;
484  return;
485 }
486 
487 void myEcalShowerPlusFiller::init(){
488 #ifdef PDEBUG
489  printf("In myEcalShowerPlusFiller::init\n");
490 #endif
491  PRINTDEBUG;
492  Clear();
493  PRINTDEBUG;
494 }
495 
496 void myEcalShowerPlusFiller::Fill(EcalShowerR* _shower, bool kShort){
497 #ifdef PDEBUG
498  printf("In myEcalShowerPlusFiller::Fill\n");
499 #endif
500  PRINTDEBUG;
501 
502  shower = _shower;
503 
504  AMSEventR *evt= AMSEventR::Head();
505 
506  if (shower) {
507 
508  PRINTDEBUG;
509 
510  copy_n(shower->Elayer_corr, 18, Elayer_corr);
511 
512  if (!kShort) {
513 
514  // this is a workaround to remove some NaN present in the pass4
515  for(int i=0;i<(int)(evt->NEcalHit());i++){
516  if(evt->pEcalHit(i)->Edep!=evt->pEcalHit(i)->Edep)
517  evt->pEcalHit(i)->Edep=0;
518  }
519 
520  //------------------------------------------------
521 
522 
523  copy_n(shower->Energy3C, 3, Energy3C);
524 
525 
526  PRINTDEBUG;
527 
528 #ifdef ECALBDT
529  BDT = shower->GetEcalBDT(5, 0, 0);
530  BDTs = shower->GetEcalBDT(5, 1, 0);
531  BDTE = shower->GetEcalBDT(5, 0, 1);
532  BDTEs = shower->GetEcalBDT(5, 1, 1);
533 
534  BDT_v4 = shower->GetEcalBDT(4);
535 #endif
536 
537  PRINTDEBUG;
538 
539 #ifdef ECALAXIS
540  EcalAxis::Version=3;
541  static EcalAxis ecalaxis;
542 #endif
543 
544  PRINTDEBUG;
545 
546 #ifdef ECALAXIS
547  // Performing before GetEcalBDTCHI2(>=3) then the EcalAxis::process will be NOT redone.
548  // If EcalAxis::process is called with another algorithm then the 'cached' value is destroyed
549  ecalaxis.process(shower,8);//use Cell method for finding COG
550  for(int ilayer=0; ilayer<18; ilayer++) EcalAxis_chi2_8[ilayer] = ecalaxis.ecalchi2->get_chi2(ilayer,0);
551 #endif
552 
553  PRINTDEBUG;
554 
555 #ifdef ECALBDT
556  BDTChi2 = shower->GetEcalBDTCHI2(3, 0, 0);
557  BDTChi2s = shower->GetEcalBDTCHI2(3, 1, 0);
558  BDTChi2E = shower->GetEcalBDTCHI2(3, 0, 1);
559  BDTChi2Es = shower->GetEcalBDTCHI2(3, 1, 1);
560 #endif
561 
562  PRINTDEBUG;
563 
564 #ifdef ECALAXIS
565  // Performing before GetEcalBDTCHI2(<3) then the EcalAxis::process will be NOT redone.
566  // If EcalAxis::process is called with another algorithm then the 'cached' value is destroyed
567  ecalaxis.process(shower,2);//use Lateral Fit algorithm for finding COG in the layer -> time consuming
568  for(int ilayer=0; ilayer<18; ilayer++) EcalAxis_chi2_2[ilayer] = ecalaxis.ecalchi2->get_chi2(ilayer,0);
569 #endif
570 
571  PRINTDEBUG;
572 
573 #ifdef ECALBDT
574  BDTChi2_v2 = shower->GetEcalBDTCHI2(2);
575 #endif
576 
577  PRINTDEBUG;
578 
579 #ifdef ECALESE
580  ESEV2 = shower->EcalStandaloneEstimatorV2();
581  ESEV3 = shower->EcalStandaloneEstimatorV3();
582  ESEv3Efficiency = shower->GetESEv3Efficiency();
583  ESEv3Rejection = shower->GetESEv3Rejection();
584 #endif
585 
586  PRINTDEBUG;
587 
588  // for (int i2Dclu=0; i2Dclu<shower->NEcal2DCluster(); i2Dclu++){
589  // Ecal2DClusterR* _2Dclu = shower->pEcal2DCluster(i2Dclu);
590  // if (_2Dclu){
591  // for (int iclu=0; iclu<_2Dclu->NEcalCluster(); iclu++){
592  // EcalClusterR* clu = _2Dclu->pEcalCluster(iclu);
593  // if (clu){
594  // for (int ihit=0; ihit<clu->NEcalHit(); ihit++){
595  // EcalHitR* hit = clu->pEcalHit(ihit);
596  // if (hit){
597  // for (int iL=0; iL<18; iL++) {
598  // if (hit->Plane==iL) {
599  // EdepLayer[iL]+=hit->Edep;
600  // for (int ii=iL; ii<18; ii++) {
601  // EdepUpToLayer[ii]+=hit->Edep;
602  // }
603  // }
604  // }
605  // }
606  // }
607  // }
608  // }
609  // }
610  // }
611 
612  //for (int i2Dclu=0; i2Dclu<shower->NEcal2DCluster(); i2Dclu++){
613  // Ecal2DClusterR* _2Dclu = shower->pEcal2DCluster(i2Dclu);
614  // if (_2Dclu){
615  // for (int iclu=0; iclu<_2Dclu->NEcalCluster(); iclu++){
616  // EcalClusterR* clu = _2Dclu->pEcalCluster(iclu);
617  // if (clu){
618  // for (int iL=0; iL<18; iL++) {
619  // if (clu->Plane==iL) {
620  // EdepLayer[iL]+=clu->Edep;
621  // }
622  // }
623  // }
624  // }
625  // }
626  //}
627 
628  //Build Hit Map
629  for(int i2dcl=0; i2dcl<shower->NEcal2DCluster(); i2dcl++) {
630  Ecal2DClusterR *ecal2dcl = shower->pEcal2DCluster(i2dcl);
631  if (ecal2dcl) {
632  for(int i1dcl=0; i1dcl<ecal2dcl->NEcalCluster(); i1dcl++) {
633  EcalClusterR *ecal1dcl = ecal2dcl->pEcalCluster(i1dcl);
634  if (ecal1dcl) {
635  if( ecal1dcl->Plane >= 0 and ecal1dcl->Plane < 18 )
636  Elayer[ecal1dcl->Plane] += ecal1dcl->Edep;
637  for(int ihit=0; ihit<ecal1dcl->NEcalHit(); ihit++){
638  EcalHitR *hit = ecal1dcl->pEcalHit(ihit);
639  if (hit) {
640  // Status bit 5 (6th) means "used".
641  // Threshold to 4 ADC to high gain is essentially the same as in the zero suppression:
642  // any residual hit, after zero suppression, is however negligible and useless.
643  // Suggestions from M. Incagli, on 08/01/2014
644  if ((hit->Status & 32) && hit->ADC[0] > 4) {
645  AddResults(hit);
646  }
647  }
648  else {
649  Warning(__func__, "Run=%d, Event=%d) pEcalHit(%d) does not exist...\n", myEvent::gethead()->Run, myEvent::gethead()->Event, ihit);
650  }
651  }
652  }
653  else {
654  Warning(__func__, "Run=%d, Event=%d) pEcalCluster(%d) does not exist...\n", myEvent::gethead()->Run, myEvent::gethead()->Event, i1dcl);
655  }
656  }
657  }
658  else {
659  Warning(__func__, "Run=%d, Event=%d) pEcal2DCluster(%d) does not exist...\n", myEvent::gethead()->Run, myEvent::gethead()->Event, i2dcl);
660  }
661  }
662 
663  PRINTDEBUG;
664 
665  }
666  else { //kShort
667  PRINTDEBUG;
668  BDT = shower->GetEcalBDT(5, 0, 0);
669  PRINTDEBUG;
670  }
671 
672  PRINTDEBUG;
673  }
674 
675  return;
676 }
677 
678 void myEcalShowerPlusFiller::AddResults(EcalHitR* hit) {
679 
680  if (hit) {
681  ecalhit adct;
682  adct.Edep = hit->Edep;
683  adct.Status = hit->Status;
684  //memcpy(adct.ADC, hit->ADC, 3*sizeof(hit->ADC[0]));
685  std::copy(hit->ADC, hit->ADC+3, adct.ADC);
686  _ecalhitmap.insert(std::pair<short int, ecalhit>(ComputeIndex(hit->Cell, hit->Plane), adct));
687  }
688 
689  return;
690 }
691 
692 #endif //#ifdef _WITHGBATCH_
693 
Generated on Tue Dec 29 2015 23:02:50 for AMSDST by   doxygen 1.8.3