/[PAMELA software]/calo/flight/CaloEnergy/src/CaloEnergy.cpp
ViewVC logotype

Contents of /calo/flight/CaloEnergy/src/CaloEnergy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.13 - (show annotations) (download)
Thu Aug 20 09:02:13 2009 UTC (15 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.12: +378 -234 lines
Cleanup and new features

1 #include <CaloEnergy.h>
2 #include <PamLevel2.h>
3
4 //--------------------------------------
5 /**
6 * Default constructor
7 */
8 CaloEnergy::CaloEnergy(){
9 Clear();
10 }
11
12 CaloEnergy::CaloEnergy(PamLevel2 *l2p){
13 //
14 Clear();
15 //
16 L2 = l2p;
17 //
18 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
19 //
20 fSimu = false;
21 this->Set();
22 //
23 }
24
25 CaloEnergy::CaloEnergy(PamLevel2 *l2p, Bool_t simulation){
26 //
27 Clear();
28 //
29 L2 = l2p;
30 //
31 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
32 //
33 fSimu = simulation;
34 this->Set();
35 //
36 }
37
38 void CaloEnergy::Delete(){
39 Clear();
40 delete this;
41 }
42
43 void CaloEnergy::UseLevel2(){
44 if ( clong ){
45 delete clong;
46 clong = new CaloLong(L2);
47 clong->SplitInto(0,22);
48 };
49 if ( cp ) delete cp;
50 cp = NULL;
51 }
52
53 void CaloEnergy::UseCaloPreSampler(){
54 //
55 // use the presampler setting forcefitmode to 1000 means to force the DV routine to find the track inside the calorimeter using the "shower" approach developed for electrons
56 //
57 if ( !cp ) cp = new CaloPreSampler(L2);
58 cp->SplitInto(0,22);
59 cp->SetForceFitMode(1000);
60 // cp->UseTracker(false);
61 // cp->ForceCaloFit();
62 // cp->SetDebug(true);
63 // cp->Process();
64 if ( clong ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer());
65 }
66
67
68 void CaloEnergy::UseLongFit(){
69 fPl = 0;
70 fLong = true;
71 if ( !clong ){
72 clong = new CaloLong(L2);
73 if ( cp ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer());
74 clong->SplitInto(0,22);
75 };
76 //
77 }
78
79 void CaloEnergy::Set(){
80 //
81 // set default values, NB default conversion factor for energy is just very approximated!
82 //
83 OBT = 0;
84 PKT = 0;
85 atime = 0;
86 sntr = "start";
87 //
88 AOBT = 0;
89 APKT = 0;
90 aatime = 0;
91 asntr = "start";
92 //
93 debug = false;
94 //
95 indep = false;
96 //
97 fAllpl = true;
98 fLong = false;
99 fPl = 1;
100 fRad = -1;
101 cp = NULL;
102 clong = NULL;
103 x0max = -1.;
104 //
105 multicol = false;
106 //
107 this->DefineGeometry();
108 fXosel =true;
109 fXesel = true;
110 fYosel = true;
111 fYesel = true;
112 fConv_rxe = 44.4;
113 fConv_rxo = 44.4;
114 fConv_ryo = 44.4;
115 fConv_rye = 44.4;
116 fXomin = 1000;
117 fXemin = 1000;
118 fYomin = 1000;
119 fYemin = 1000;
120 //
121 }
122 void CaloEnergy::DefineGeometry(){
123 //
124 // Use CaloStrip to determine once the position of border strips for each section
125 //
126 // fM = 2. + 0.096; // real position from cbar BUG the 0.096 is already taken into account in the border calculation made by Giovanna
127 fM = 2. ; // real position from cbar
128 // fM1 = 2. - 0.122 - 0.096; // due to calculation of xe1 etc. BUG! this way we count from the silicon border not from the silicon sensitive area
129 fM1 = 2. - 0.122 - 0.096 + 0.096; // due to calculation of xe1 etc.
130 if ( fM1 < 0. ) fM1 = 0.;
131 //
132 CaloStrip *cs = new CaloStrip(fSimu);
133 //
134 // view y plane 0 strip 0
135 cs->Set(1,0,0);
136 xe1= cs->GetY();
137 // view y plane 0 strip 31
138 cs->Set(1,0,31);
139 xe2= cs->GetY();
140 // view y plane 0 strip 32
141 cs->Set(1,0,32);
142 xe3= cs->GetY();
143 // view y plane 0 strip 63
144 cs->Set(1,0,63);
145 xe4= cs->GetY();
146 // view y plane 0 strip 64
147 cs->Set(1,0,64);
148 xe5= cs->GetY();
149 // view y plane 0 strip 95
150 cs->Set(1,0,95);
151 xe6= cs->GetY();
152 // view x plane 0 strip 0
153 cs->Set(0,0,0);
154 yo1= cs->GetX();
155 // view x plane 0 strip 31
156 cs->Set(0,0,31);
157 yo2= cs->GetX();
158 // view x plane 0 strip 32
159 cs->Set(0,0,32);
160 yo3= cs->GetX();
161 // view x plane 0 strip 63
162 cs->Set(0,0,63);
163 yo4= cs->GetX();
164 // view x plane 0 strip 64
165 cs->Set(0,0,64);
166 yo5= cs->GetX();
167 // view x plane 0 strip 95
168 cs->Set(0,0,95);
169 yo6= cs->GetX();
170 // view y plane 1 strip 0
171 cs->Set(1,1,0);
172 xo1= cs->GetY();
173 // view y plane 1 strip 31
174 cs->Set(1,1,31);
175 xo2= cs->GetY();
176 // view y plane 1 strip 32
177 cs->Set(1,1,32);
178 xo3= cs->GetY();
179 // view y plane 1 strip 63
180 cs->Set(1,1,63);
181 xo4= cs->GetY();
182 // view y plane 1 strip 64
183 cs->Set(1,1,64);
184 xo5= cs->GetY();
185 // view y plane 1 strip 95
186 cs->Set(1,1,95);
187 xo6= cs->GetY();
188 // view x plane 1 strip 0
189 cs->Set(0,1,0);
190 ye1= cs->GetX();
191 // view x plane 1 strip 31
192 cs->Set(0,1,31);
193 ye2= cs->GetX();
194 // view x plane 1 strip 32
195 cs->Set(0,1,32);
196 ye3= cs->GetX();
197 // view x plane 1 strip 63
198 cs->Set(0,1,63);
199 ye4= cs->GetX();
200 // view x plane 1 strip 64
201 cs->Set(0,1,64);
202 ye5= cs->GetX();
203 // view x plane 1 strip 95
204 cs->Set(0,1,95);
205 ye6= cs->GetX();
206 //
207 for (Int_t p = 0; p<22; p ++){
208 for (Int_t v = 0; v<2; v++ ){
209 cs->Set(v,p,0);
210 trk_z[p][v]= cs->GetZ(); // Z coord for each plane
211 };
212 };
213 //
214 delete cs;
215 //
216 }
217
218 void CaloEnergy::Clear(){
219 //
220 // clear variables
221 //
222 fPartsel = false;
223 fSel = false;
224 fXosel = false;
225 fXesel = false;
226 fYosel = false;
227 fYesel = false;
228 fCount = 0.;
229 fEnergy = 0.;
230 fEnergyxe = 0.;
231 fEnergyxo = 0.;
232 fEnergyye = 0.;
233 fEnergyyo = 0.;
234 fMax_plane = 0;
235 fMax_planexo = 0;
236 fMax_planexe = 0;
237 fMax_planeyo = 0;
238 fMax_planeye = 0;
239 xomax_en= 0.;
240 xemax_en= 0.;
241 yomax_en= 0.;
242 yemax_en= 0.;
243 //
244 memset(enstrip,0,2*22*96*(sizeof(Float_t)));
245 en = 0.;
246 view = 0;
247 plane = 0;
248 strip = 0;
249 energyxe = 0.;
250 energyyo = 0.;
251 energyxo = 0.;
252 energyye = 0.;
253 fYoout = 0;
254 fYeout = 0;
255 fXoout = 0;
256 fXeout = 0;
257 fXEen_maxplane = 0.;
258 fXOen_maxplane = 0.;
259 fYEen_maxplane = 0.;
260 fYOen_maxplane = 0.;
261 memset(en_xep,0,11*sizeof(Float_t));
262 memset(en_yep,0,11*sizeof(Float_t));
263 memset(en_xop,0,11*sizeof(Float_t));
264 memset(en_yop,0,11*sizeof(Float_t));
265 //
266 fColumn = -1;
267 fColXE = -1;
268 fColXO = -1;
269 fColYE = -1;
270 fColYO = -1;
271 memset(encol,0,2*3*sizeof(Float_t));
272 entot[0] = 0.;
273 entot[1] = 0.;
274 //
275 }
276
277 void CaloEnergy::Print(){
278 //
279 printf("========================================================================\n");
280 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
281 printf(" fEnergy :.............. %f \n",fEnergy);
282 printf(" fMax_plane :........... %f \n",fMax_plane);
283 printf(" fMax_planexo :......... %i \n",fMax_planexo);
284 printf(" fMax_planexe :......... %i \n",fMax_planexe);
285 printf(" fMax_planeyo :......... %i \n",fMax_planeyo);
286 printf(" fMax_planeye :......... %i \n",fMax_planeye);
287 printf(" fCount :.............. %f \n",fCount);
288 printf(" fSel :.............. %i \n",fSel);
289 printf(" fPartsel:.............. %i \n",fPartsel);
290 printf(" fXesel :.............. %i \n",fXesel);
291 printf(" fXosel :.............. %i \n",fXosel);
292 printf(" fYesel :.............. %i \n",fYesel);
293 printf(" fYosel :.............. %i \n",fYosel);
294 printf(" fXemin :.............. %i \n",fXemin);
295 printf(" fXomin :.............. %i \n",fXomin);
296 printf(" fYemin :.............. %i \n",fYemin);
297 printf(" fYomin :.............. %i \n",fYomin);
298 printf(" fXeout :.............. %i \n",fXeout);
299 printf(" fXoout :.............. %i \n",fXoout);
300 printf(" fYeout :.............. %i \n",fYeout);
301 printf(" fYoout :.............. %i \n",fYoout);
302 printf(" fSimu :.............. %i \n",fSimu);
303 printf(" fM :.............. %f \n",fM);
304 printf(" fM1 :.............. %f \n",fM1);
305 printf(" fRad :.............. %i \n",fRad);
306 printf(" fPl :.............. %i \n",fPl);
307 printf(" fColumn :.............. %i \n",fColumn);
308 printf(" multicol:.............. %i \n",multicol);
309 printf(" encol x :.............. %f \n",this->GetEncol(0));
310 printf(" encol y :.............. %f \n",this->GetEncol(1));
311 printf(" entot x :.............. %f \n",this->GetEntot(0));
312 printf(" entot y :.............. %f \n",this->GetEntot(1));
313 printf(" fColXE :.............. %i \n",fColXE);
314 printf(" fColXO :.............. %i \n",fColXO);
315 printf(" fColYE :.............. %i \n",fColYE);
316 printf(" fColYO :.............. %i \n",fColYO);
317 printf(" fConv_rxe ............. %f \n",fConv_rxe);
318 printf(" fConv_rxo ............. %f \n",fConv_rxo);
319 printf(" fConv_ryo ............. %f \n",fConv_ryo);
320 printf(" fConv_rye ............. %f \n",fConv_rye);
321 printf(" fLong :.............. %i \n",fLong);
322 printf(" energyxe:.............. %f \n",energyxe);
323 printf(" energyxo:.............. %f \n",energyxo);
324 printf(" energyye:.............. %f \n",energyye);
325 printf(" energyyo:.............. %f \n",energyyo);
326 printf(" fXEen_maxplane:........ %f \n",fXEen_maxplane);
327 printf(" fXOen_maxplane:........ %f \n",fXOen_maxplane);
328 printf(" fYEen_maxplane:........ %f \n",fYEen_maxplane);
329 printf(" fYOen_maxplane:........ %f \n",fYOen_maxplane);
330 printf(" x0max :.............. %f \n",x0max);
331 printf(" debug :.............. %i \n",debug);
332
333 printf("========================================================================\n");
334 //
335 }
336
337 void CaloEnergy::SetMinimumContainment(TString section, Int_t plane){
338 section.ToUpper();
339 if ( section.Contains("XO") ) fXomin = plane;
340 if ( section.Contains("XE") ) fXemin = plane;
341 if ( section.Contains("YO") ) fYomin = plane;
342 if ( section.Contains("YE") ) fYemin = plane;
343 }
344
345 void CaloEnergy::SetMinimumContainment(Int_t plane){
346 this->SetMinimumContainment("XEXOYEYO",plane);
347 }
348
349 void CaloEnergy::SetConversionFactor(TString section, Float_t conv){
350 section.ToUpper();
351 if ( section.Contains("XO") ) fConv_rxo = conv;
352 if ( section.Contains("XE") ) fConv_rxe = conv;
353 if ( section.Contains("YO") ) fConv_ryo = conv;
354 if ( section.Contains("YE") ) fConv_rye = conv;
355 }
356
357 void CaloEnergy::SetConversionFactor(Float_t conv){
358 this->SetConversionFactor("XEXOYEYO",conv);
359 }
360
361 Int_t CaloEnergy::GetMinimumContainment(TString section){
362 section.ToUpper();
363 if ( section.Contains("XO") ) return(fXomin);
364 if ( section.Contains("XE") ) return(fXemin);
365 if ( section.Contains("YE") ) return(fYemin);
366 if ( section.Contains("YO") ) return(fYomin);
367 printf(" ERROR: section not recognized \n");
368 return(-1000);
369 }
370
371 Float_t CaloEnergy::GetConversionFactor(TString section){
372 section.ToUpper();
373 if ( section.Contains("XO") ) return(fConv_rxo);
374 if ( section.Contains("XE") ) return(fConv_rxe);
375 if ( section.Contains("YO") ) return(fConv_ryo);
376 if ( section.Contains("YE") ) return(fConv_rye);
377 printf(" ERROR: section not recognized \n");
378 return(-1000.);
379 }
380
381 Int_t CaloEnergy::GetMaxplane(TString section){
382 section.ToUpper();
383 if ( section.Contains("XO") ) return fMax_planexo;
384 if ( section.Contains("XE") ) return fMax_planexe;
385 if ( section.Contains("YO") ) return fMax_planeyo;
386 if ( section.Contains("YE") ) return fMax_planeye;
387 return(-1);
388 }
389
390 Int_t CaloEnergy::GetColumn(TString section){
391 section.ToUpper();
392 if ( section.Contains("XO") ) return fColXO;
393 if ( section.Contains("XE") ) return fColXE;
394 if ( section.Contains("YO") ) return fColYO;
395 if ( section.Contains("YE") ) return fColYE;
396 return(-1);
397 }
398
399 Float_t CaloEnergy::GetMipEnergyAtMaxplane(TString section){
400 printf(" WARNING: OBSOLETE METHOD, use GetMipEnergyAtMaxplane(TString) instead! \n");
401 return (this->GetEnergyAtMaxplane(section));
402 }
403
404 Float_t CaloEnergy::GetEnergyAtMaxplane(TString section){
405 section.ToUpper();
406 if ( section.Contains("XO") ) return xomax_en;
407 if ( section.Contains("XE") ) return xemax_en;
408 if ( section.Contains("YO") ) return yomax_en;
409 if ( section.Contains("YE") ) return yemax_en;
410 return(-1);
411 }
412
413 Float_t CaloEnergy::GetMaxEnergy(TString section){
414 printf(" WARNING: OBSOLETE METHOD, use GetMipEnergy(TString) instead! \n");
415 return (this->GetMipEnergy(section));
416 }
417
418 Float_t CaloEnergy::GetMipEnergy(TString section){
419 section.ToUpper();
420 if ( fLong ){
421 this->Process(section);
422 return fXOen_maxplane;
423 } else {
424 if ( section.Contains("XO") ) return fXOen_maxplane;
425 if ( section.Contains("XE") ) return fXEen_maxplane;
426 if ( section.Contains("YO") ) return fYOen_maxplane;
427 if ( section.Contains("YE") ) return fYEen_maxplane;
428 };
429 return(-1);
430 }
431
432 Float_t CaloEnergy::GetEncol(Int_t i){
433 if ( fColumn > -1 && (((fXesel || fXosel)&&i==1) || ((fYesel || fYosel)&&i==0)) ){
434 Int_t t = -1;
435 if ( i == 0 ){
436 if ( fColumn == 0 || fColumn == 3 || fColumn == 6 ) t = 0;
437 if ( fColumn == 1 || fColumn == 4 || fColumn == 7 ) t = 1;
438 if ( fColumn == 2 || fColumn == 5 || fColumn == 8 ) t = 2;
439 } else {
440 if ( fColumn == 0 || fColumn == 1 || fColumn == 2 ) t = 0;
441 if ( fColumn == 3 || fColumn == 4 || fColumn == 5 ) t = 1;
442 if ( fColumn == 6 || fColumn == 7 || fColumn == 8 ) t = 2;
443 };
444 if ( debug ) printf(" encol: i %i t %i encol %f \n",i,t,encol[i][t]);
445 return encol[i][t];
446 };
447 return(-1.);
448 }
449
450 Float_t CaloEnergy::GetMaxEnergy(){
451 printf(" WARNING: OBSOLETE METHOD, use GetMipEnergy() instead! \n");
452 return (this->GetMipEnergy());
453 }
454
455 Float_t CaloEnergy::GetMipEnergy(){
456 if ( fLong ){
457 if ( debug ) printf(" oh! call process! with asntr %s and sntr %s \n",asntr.Data(),sntr.Data());
458 this->Process(asntr);
459 };
460 return((fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane));
461 }
462
463
464 Bool_t CaloEnergy::IsInsideAcceptance(TString section){
465 //
466 // check if the event is inside the acceptance of the given section(s)
467 //
468 TString ntr = section;
469 if ( !L2 ){
470 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
471 printf(" ERROR: CaloEnergy variables not filled \n");
472 return false;
473 };
474 //
475 Bool_t newentry = false;
476 //
477 if ( L2->IsORB() ){
478 if ( L2->GetOrbitalInfo()->pkt_num != APKT || L2->GetOrbitalInfo()->OBT != AOBT || L2->GetOrbitalInfo()->absTime != aatime || strcmp(ntr.Data(),asntr.Data()) ){
479 newentry = true;
480 AOBT = L2->GetOrbitalInfo()->OBT;
481 APKT = L2->GetOrbitalInfo()->pkt_num;
482 aatime = L2->GetOrbitalInfo()->absTime;
483 asntr = ntr;
484 };
485 } else {
486 newentry = true;
487 };
488 //
489 // if we have already called this method for this event and no input changed then return fSel and exit
490 //
491 if ( !newentry ) return fSel;
492 //
493 // process the event
494 //
495 if ( debug ) printf(" ########## IsInsideAcceptance ######### \n");
496 //
497 // clear variables
498 //
499 this->Clear();
500 //
501 section.ToUpper();
502 //
503 // Count the number of section(s) given as input
504 //
505 Int_t fNumSec = Int_t(section.Contains("XO"))+Int_t(section.Contains("XE"))+Int_t(section.Contains("YO"))+Int_t(section.Contains("YE"));
506 if ( !fNumSec ){
507 printf(" ERROR: section must be XO or XE or YO or YE while it is %s \n",section.Data());
508 return false;
509 };
510 //
511 // If the presampler object exists then use the presampler output instead of the level2 output
512 //
513 CaloLevel2 *cl2 = NULL;
514 if ( cp ){
515 cl2 = cp->GetCaloLevel2();
516 } else {
517 cl2 = L2->GetCaloLevel2();
518 };
519 //
520 // get the energy for every strip of the calorimeter
521 //
522 for (Int_t ch=0; ch< L2->GetCaloLevel1()->istrip; ch++){
523 en = L2->GetCaloLevel1()->DecodeEstrip(ch,view,plane,strip);
524 enstrip[view][plane][strip]=en;
525 };
526 //
527 if ( debug && ((fM1+0.122-0.244*(Float_t)fRad) < 0.) ) printf("Error: (fM1+0.122-0.244*(Float_t)fRad) < 0. fM1 %f fRad %i %f \n",fM1,fRad,(fM1+0.122-0.244*(Float_t)fRad));
528 //
529 // sum energy plane by plane for each sections
530 //
531 Float_t fen_xep[11];
532 Float_t fen_xop[11];
533 Float_t fen_yep[11];
534 Float_t fen_yop[11];
535 memset(fen_xep,0,11*sizeof(Float_t));
536 memset(fen_xop,0,11*sizeof(Float_t));
537 memset(fen_yep,0,11*sizeof(Float_t));
538 memset(fen_yop,0,11*sizeof(Float_t));
539 //
540 for (Int_t i=0;i<11;i++){
541 for(strip=0; strip<96; strip++) {
542 fen_xep[i] += enstrip[1][2*i][strip];
543 fen_yop[i] += enstrip[0][2*i][strip];
544 fen_xop[i] += enstrip[1][(2*i)+1][strip];
545 fen_yep[i] += enstrip[0][(2*i)+1][strip];
546 if ( fRad < 0 ){
547 //
548 // run over all the strips of the plane
549 //
550 en_xep[i] += enstrip[1][2*i][strip];
551 en_yop[i] += enstrip[0][2*i][strip];
552 en_xop[i] += enstrip[1][(2*i)+1][strip];
553 en_yep[i] += enstrip[0][(2*i)+1][strip];
554 } else {
555 //
556 // use only the strips inside a cylinder of given radius fRad
557 //
558 if ( cl2->cibar[2*i][1] >= 1 && cl2->cibar[2*i][1] <= 96 &&
559 (strip >= (cl2->cibar[2*i][1]-1-fRad)) && (strip <= (cl2->cibar[2*i][1]-1+fRad)) ) en_xep[i] += enstrip[1][2*i][strip];
560
561 if ( cl2->cibar[2*i][0] >= 1 && cl2->cibar[2*i][0] <= 96 &&
562 (strip >= (cl2->cibar[2*i][0]-1-fRad)) && (strip <= (cl2->cibar[2*i][0]-1+fRad)) ) en_yop[i] += enstrip[0][2*i][strip];
563
564 if ( cl2->cibar[(2*i)+1][1] >= 1 && cl2->cibar[(2*i)+1][1] <= 96 &&
565 (strip >= (cl2->cibar[(2*i)+1][1]-1-fRad)) && (strip <= (cl2->cibar[(2*i)+1][1]-1+fRad)) ) en_xop[i] += enstrip[1][(2*i)+1][strip];
566
567 if ( cl2->cibar[(2*i)+1][0] >= 1 && cl2->cibar[(2*i)+1][0] <= 96 &&
568 (strip >= (cl2->cibar[(2*i)+1][0]-1-fRad)) && (strip <= (cl2->cibar[(2*i)+1][0]-1+fRad)) ) en_yep[i] += enstrip[0][(2*i)+1][strip];
569 };
570 };
571 if ( debug ) printf(" ex_xep[%i] %f cibar %i \n",i,en_xep[i],cl2->cibar[2*i][1]);
572 if ( debug ) printf(" ex_xop[%i] %f cibar %i \n",i,en_xop[i],cl2->cibar[(2*i)+1][1]);
573 if ( debug ) printf(" ex_yep[%i] %f cibar %i \n",i,en_yep[i],cl2->cibar[(2*i)+1][0]);
574 if ( debug ) printf(" ex_yop[%i] %f cibar %i \n",i,en_yop[i],cl2->cibar[2*i][0]);
575 energyxe += en_xep[i];
576 energyyo += en_yop[i];
577 energyxo += en_xop[i];
578 energyye += en_yep[i];
579 };
580 //
581 // Find the plane of maximum for each section
582 //
583 //
584 Int_t xen = 0;
585 Int_t yon = 0;
586 Int_t xon = 0;
587 Int_t yen = 0;
588 Float_t en = 0.;
589 //
590 if ( section.Contains("XE") ){
591 yon++;
592 xon++;
593 yen++;
594 for (Int_t ipl =0; ipl < 11; ipl ++) {
595 en = fen_xep[ipl];
596 if ( !fAllpl ) en = en_xep[ipl];
597 if(en > xemax_en) {
598 xemax_en = en;
599 fMax_planexe = ipl;
600 };
601 };
602 };
603 //
604 if ( section.Contains("YO") ){
605 xon++;
606 yen++;
607 for (Int_t ipl =0; ipl < 11; ipl ++) {
608 en = fen_yop[ipl];
609 if ( !fAllpl ) en = en_yop[ipl];
610 if(en > yomax_en) {
611 yomax_en = en;
612 fMax_planeyo = ipl;
613 };
614 };
615 };
616 //
617 if ( section.Contains("XO") ){
618 yen++;
619 for (Int_t ipl =0; ipl < 11; ipl ++) {
620 en = fen_xop[ipl];
621 if ( !fAllpl ) en = en_xop[ipl];
622 if(en > xomax_en) {
623 xomax_en = en;
624 fMax_planexo = ipl;
625 };
626 };
627 };
628 //
629 if ( section.Contains("YE") ){
630 for (Int_t ipl =0; ipl < 11; ipl ++) {
631 en = fen_yep[ipl];
632 if ( !fAllpl ) en = en_yep[ipl];
633 if(en > yemax_en) {
634 yemax_en = en;
635 fMax_planeye = ipl;
636 };
637 };
638 };
639 //
640 // the maximum is given externally as X0, convert it to plane and section
641 //
642 if ( x0max > 0. ){
643 if ( debug ) printf(" CALCULATE MAX PLANE GIVEN X0 ASSUMING PERPENDICULAR TRACK \n");
644 Int_t wpl = (Int_t)roundf(x0max/0.76);
645 Bool_t isY = false;
646 if ( ((x0max/0.76)-(Float_t)wpl) > 0. ) isY = true;
647 xomax_en = 0.;
648 yemax_en = 0.;
649 xemax_en = 0.;
650 yomax_en = 0.;
651 //
652 if ( !(wpl%2) ){
653 // 0, 2, 4, ...
654 if ( isY ){
655 if ( section.Contains("YO") ) yomax_en = 1000.;
656 if ( section.Contains("XE") ) xemax_en = 500.;
657 fMax_planeyo=wpl/2;
658 fMax_planexe=wpl/2;
659 if ( section.Contains("XO") ) xomax_en = 10.;
660 if ( section.Contains("YE") ) yemax_en = 5.;
661 } else {
662 if ( section.Contains("YO") ) yomax_en = 500.;
663 if ( section.Contains("XE") ) xemax_en = 1000.;
664 fMax_planeyo=wpl/2;
665 fMax_planexe=wpl/2;
666 if ( section.Contains("XO") ) xomax_en = 5.;
667 if ( section.Contains("YE") ) yemax_en = 10.;
668 };
669 } else {
670 // 1, 3, 5, ...
671 if ( isY ){
672 if ( section.Contains("YE") ) yemax_en = 1000.;
673 if ( section.Contains("XO") ) xomax_en = 500.;
674 fMax_planeye=(wpl-1)/2;
675 fMax_planexo=(wpl-1)/2;
676 if ( section.Contains("XE") ) xemax_en = 10.;
677 if ( section.Contains("YO") ) yomax_en = 5.;
678 } else {
679 if ( section.Contains("YE") ) yemax_en = 500.;
680 if ( section.Contains("XO") ) xomax_en = 1000.;
681 fMax_planeye=(wpl-1)/2;
682 fMax_planexo=(wpl-1)/2;
683 if ( section.Contains("XE") ) xemax_en = 5.;
684 if ( section.Contains("YO") ) yomax_en = 10.;
685 };
686 };
687 if ( debug ) printf(" x0max %f x0max/0.76 %f wpl %i isY %i yomax_en %f xemax_en %f yemax_en %f xomax_en %f fMaxplane %i %i %i %i\n",x0max,(x0max/0.76),wpl,isY,yomax_en,xemax_en,yemax_en,xomax_en,fMax_planeyo,fMax_planexe,fMax_planeye,fMax_planexo);
688 };
689 //
690 Int_t nPl = fPl;
691 //
692 // Set the maximum in case of coherent mode was selected
693 //
694 if ( !indep ){
695 nPl = 0;
696 if ( debug ) printf(" A: Check maximum, coherent mode: xoen %f yoen %f xeen %f yeen %f xomax %i yomax %i xemax %i yemax %i\n",xomax_en,yomax_en,xemax_en,yemax_en,fMax_planexo,fMax_planeyo,fMax_planexe,fMax_planeye);
697 Int_t nummod = 0;
698 Int_t numexpl = 0;
699 if ( xomax_en > xemax_en && xomax_en > yemax_en && xomax_en > yomax_en ){
700 //
701 // Section XO contains the maximum energy release per plane of the whole calorimeter
702 //
703 if ( debug ) printf(" XO is MAX %i %i %i %i\n",xen,yon,xon,yen);
704 //
705 // fMax_plane is the plane of maximum + number of additional dE/dx measurement counting planes from 0 to 43
706 //
707 fMax_plane = (fNumSec * fMax_planexo) +(Float_t)xon + fPl;
708 //
709 // nummod is the integer part of the number of modules in which the maximum is contained
710 //
711 nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
712 //
713 // numexpl is the number of additional planes (0,1,2) inside the module
714 //
715 numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
716 //
717 };
718 if ( xemax_en > xomax_en && xemax_en > yemax_en && xemax_en > yomax_en ){
719 if ( debug ) printf(" XE is MAX %i %i %i %i\n",xen,yon,xon,yen);
720 fMax_plane = (fNumSec * fMax_planexe) +(Float_t)xen + fPl;
721 nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
722 numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
723 //
724 };
725
726 if ( yemax_en > xomax_en && yemax_en > xemax_en && yemax_en > yomax_en ){
727 if ( debug ) printf(" YE is MAX %i %i %i %i\n",xen,yon,xon,yen);
728 fMax_plane = (fNumSec * fMax_planeye) +(Float_t)yen + fPl;
729 nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
730 numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
731 //
732 };
733 if ( yomax_en > xemax_en && yomax_en > yemax_en && yomax_en > xomax_en ){
734 if ( debug ) printf(" YO is MAX %i %i %i %i\n",xen,yon,xon,yen);
735 fMax_plane = (fNumSec * fMax_planeyo) +(Float_t)yon + fPl;
736 nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
737 numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
738 //
739 };
740 //
741 // find the plane up to which is necessary to integrate the energy for each section
742 //
743 Int_t a = 0;
744 Int_t b = 0;
745 Int_t c = 0;
746 if ( numexpl > xen ) a = 1;
747 if ( numexpl > yon ) b = 1;
748 if ( numexpl > xon ) c = 1;
749 fMax_planexe = nummod;
750 fMax_planeyo = nummod - 1 + a;
751 fMax_planexo = nummod - 1 + b;
752 fMax_planeye = nummod - 1 + c;
753 if ( debug ) printf(" fMax_plane %f nummod %i numexpl %i a %i b %i c %i \n",fMax_plane,nummod,numexpl,a,b,c);
754 if ( debug ) printf(" DONE: Check maximum, coherent mode: xoen %f yoen %f xeen %f yeen %f xomax %i yomax %i xemax %i yemax %i\n",xomax_en,yomax_en,xemax_en,yemax_en,fMax_planexo,fMax_planeyo,fMax_planexe,fMax_planeye);
755 };
756 //
757 // for each plane of the calorimeter find the position of the track in the direction along the strip (where we do not have a measurement from the selected plane) by looking at the plane above/below of the other view and extrapolating the trajectory to the given plane
758 //
759 //
760 Float_t tgx_cl2;
761 Float_t tgy_cl2;
762 tgx_cl2 = cl2->tanx[0];
763 tgy_cl2 = cl2->tany[0];
764 //
765 for (Int_t p=0; p<22; p++){
766 track_coordy[p][1] = cl2->cbar[p][1];
767 track_coordx[p][1] = cl2->cbar[p][0] - fabs(trk_z[p][1]-trk_z[p][0])*tgx_cl2;
768 // track_coordx[p][1] = cl2->cbar[p][0] + fabs(trk_z[p][1]-trk_z[p][0])*tgx_cl2;
769 track_coordx[p][0] = cl2->cbar[p][0];
770 track_coordy[p][0] = cl2->cbar[p][1] - fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2;
771 // track_coordy[p][0] = cl2->cbar[p][1] + fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2;
772 if ( debug ) printf(" p %i track_coordy[p][1] %f track_coordx[p][1] %f track_coordx[p][0] %f track_coordy[p][0] %f \n",p,track_coordy[p][1],track_coordx[p][1],track_coordx[p][0],track_coordy[p][0]);
773 };
774 //
775 if ( debug ) printf(" acceptance fNumSec %i tgx %f tgy %f\n",fNumSec,tgx_cl2,tgy_cl2);
776 //
777 if ( section.Contains("XO") ){
778 //
779 // find the column hit in the first plane
780 //
781 Int_t ix = -1;
782 Int_t iy = -1;
783 if ( track_coordx[(2*0)+1][1] >= (-12.054+fM) && track_coordx[(2*0)+1][1] <= (-4.246-fM) ) ix = 0;
784 if ( track_coordx[(2*0)+1][1] >= ( -4.004+fM) && track_coordx[(2*0)+1][1] <= ( 3.804-fM) ) ix = 1;
785 if ( track_coordx[(2*0)+1][1] >= ( 4.046+fM) && track_coordx[(2*0)+1][1] <= (11.854-fM) ) ix = 2;
786 if ( cl2->cbar[(2*0)+1][1] >= (xo1 + fM1) && cl2->cbar[(2*0)+1][1] <= (xo2 - fM1) ) iy = 0;
787 if ( cl2->cbar[(2*0)+1][1] >= (xo3 + fM1) && cl2->cbar[(2*0)+1][1] <= (xo4 - fM1) ) iy = 1;
788 if ( cl2->cbar[(2*0)+1][1] >= (xo5 + fM1) && cl2->cbar[(2*0)+1][1] <= (xo6 - fM1) ) iy = 2;
789 if ( ix > -1 && iy > -1 ) fColXO = ix + iy*3;
790 //
791 // check event is inside XO acceptance, if multicol is false (SingleColumn mode) then the track must be contained in a column.
792 //
793 for (Int_t i=0; i<11; i++) {
794 if ((
795 ( track_coordx[(2*i)+1][1] >= (-12.054+fM) && track_coordx[(2*i)+1][1] <= (-4.246-fM) && (ix == 0 || multicol) ) ||
796 ( track_coordx[(2*i)+1][1] >= ( -4.004+fM) && track_coordx[(2*i)+1][1] <= ( 3.804-fM) && (ix == 1 || multicol) ) ||
797 ( track_coordx[(2*i)+1][1] >= ( 4.046+fM) && track_coordx[(2*i)+1][1] <= (11.854-fM) && (ix == 2 || multicol) )
798 ) && (
799 ( cl2->cbar[(2*i)+1][1] >= (xo1 + fM1) && cl2->cbar[(2*i)+1][1] <= (xo2 - fM1) && (iy == 0 || multicol) ) ||
800 ( cl2->cbar[(2*i)+1][1] >= (xo3 + fM1) && cl2->cbar[(2*i)+1][1] <= (xo4 - fM1) && (iy == 1 || multicol) ) ||
801 ( cl2->cbar[(2*i)+1][1] >= (xo5 + fM1) && cl2->cbar[(2*i)+1][1] <= (xo6 - fM1) && (iy == 2 || multicol) )
802 )){
803 fXosel = true;
804 fXoout = i;
805 } else {
806 fXosel = false;
807 break;
808 };
809 };
810 //
811 // if it goes out of the acceptance BUT the plane up to which we are integrating the energy is contained then the event is "partially" contained
812 //
813 if ( !fXosel && fXoout >= fXomin && fXoout >= (Int_t)(fMax_planexo+nPl) ){
814 if ( debug ) printf(" XO: this event is only partially contained: fXoout %i fXomin %i fMax_planexo + nPl %i \n",fXoout,fXomin,(Int_t)(fMax_planexo+nPl));
815 fPartsel = true;
816 fXosel = true;
817 };
818 //
819 // event is contained (or partially contained) hence we can integrate energy up to the maximum and calculate the energy as measured by this section
820 //
821 if ( fXosel ){
822 for (Int_t iplm=0; iplm<=TMath::Min(10,(Int_t)(fMax_planexo+nPl)); iplm++){
823 fXOen_maxplane += en_xop[iplm];
824 if ( debug ) printf(" XO iplm %i fXOen_maxplane %f en_xop[iplm] %f\n",iplm,fXOen_maxplane,en_xop[iplm]);
825 };
826 fEnergyxo = fXOen_maxplane/fConv_rxo;
827 //
828 for (Int_t i=0;i<11;i++){
829 for(strip=0; strip<96; strip++) {
830 //
831 // run over all the strips of the plane
832 //
833 if ( strip >= 0 && strip < 32 ) encol[1][0] += enstrip[1][(2*i)+1][strip];
834 if ( strip >= 32 && strip < 64 ) encol[1][1] += enstrip[1][(2*i)+1][strip];
835 if ( strip >= 64 && strip < 96 ) encol[1][2] += enstrip[1][(2*i)+1][strip];
836 entot[1] += enstrip[1][(2*i)+1][strip];
837 //
838 };
839 };
840 };
841 };
842 //
843 if ( section.Contains("XE") ){
844 //
845 // find the column hit in the first plane
846 //
847 Int_t ix = -1;
848 Int_t iy = -1;
849 if ( track_coordx[(2*0)][1] >= (-11.854+fM) && track_coordx[(2*0)][1] <= (-4.046-fM) ) ix = 0;
850 if ( track_coordx[(2*0)][1] >= ( -3.804+fM) && track_coordx[(2*0)][1] <= ( 4.004-fM) ) ix = 1;
851 if ( track_coordx[(2*0)][1] >= ( 4.246+fM) && track_coordx[(2*0)][1] <= (12.054-fM) ) ix = 2;
852 if ( cl2->cbar[(2*0)][1] >= (xe1 + fM1) && cl2->cbar[(2*0)][1] <= (xe2 - fM1) ) iy = 0;
853 if ( cl2->cbar[(2*0)][1] >= (xe3 + fM1) && cl2->cbar[(2*0)][1] <= (xe4 - fM1) ) iy = 1;
854 if ( cl2->cbar[(2*0)][1] >= (xe5 + fM1) && cl2->cbar[(2*0)][1] <= (xe6 - fM1) ) iy = 2;
855 if ( ix > -1 && iy > -1 ) fColXE = ix + iy*3;
856 //
857 // check event is inside XO acceptance
858 //
859 for (Int_t i=0; i<11; i++) {
860 if ((
861 ( track_coordx[(2*i)][1] >= (-11.854+fM) && track_coordx[(2*i)][1] <= (-4.046-fM) && (ix == 0 || multicol) ) ||
862 ( track_coordx[(2*i)][1] >= ( -3.804+fM) && track_coordx[(2*i)][1] <= ( 4.004-fM) && (ix == 1 || multicol) ) ||
863 ( track_coordx[(2*i)][1] >= ( 4.246+fM) && track_coordx[(2*i)][1] <= (12.054-fM) && (ix == 2 || multicol) )
864 ) && (
865 ( cl2->cbar[(2*i)][1] >= (xe1 + fM1) && cl2->cbar[(2*i)][1] <= (xe2 - fM1) && (iy == 0 || multicol) ) ||
866 ( cl2->cbar[(2*i)][1] >= (xe3 + fM1) && cl2->cbar[(2*i)][1] <= (xe4 - fM1) && (iy == 1 || multicol) ) ||
867 ( cl2->cbar[(2*i)][1] >= (xe5 + fM1) && cl2->cbar[(2*i)][1] <= (xe6 - fM1) && (iy == 2 || multicol) )
868 )){
869 fXesel = true;
870 fXeout = i;
871 } else {
872 fXesel = false;
873 break;
874 };
875 };
876 //
877 if ( !fXesel && fXeout >= fXemin && fXeout >= (Int_t)(fMax_planexe+nPl) ){
878 if ( debug ) printf(" XE: this event is only partially contained: fXeout %i fXemin %i fMax_planexe + nPl %i \n",fXeout,fXemin,(Int_t)(fMax_planexe+nPl));
879 fPartsel = true;
880 fXesel = true;
881 };
882 if ( fXesel ){
883 for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planexe+nPl)) ;iplm++){
884 fXEen_maxplane += en_xep[iplm];
885 if ( debug ) printf(" XE iplm %i fXOen_maxplane %f en_xop[iplm] %f\n",iplm,fXEen_maxplane,en_xep[iplm]);
886 };
887 fEnergyxe = fXEen_maxplane/fConv_rxe;
888 //
889 for (Int_t i=0;i<11;i++){
890 for(strip=0; strip<96; strip++) {
891 //
892 // run over all the strips of the plane
893 //
894 if ( strip >= 0 && strip < 32 ) encol[1][0] += enstrip[1][(2*i)][strip];
895 if ( strip >= 32 && strip < 64 ) encol[1][1] += enstrip[1][(2*i)][strip];
896 if ( strip >= 64 && strip < 96 ) encol[1][2] += enstrip[1][(2*i)][strip];
897 entot[1] += enstrip[1][(2*i)][strip];
898 //
899 };
900 };
901 };
902 };
903 //
904 if ( section.Contains("YE") ){
905 //
906 // find the column hit in the first plane
907 //
908 Int_t ix = -1;
909 Int_t iy = -1;
910 if ( track_coordy[(2*0)+1][0] >= (-12.154+fM) && track_coordy[(2*0)+1][0] <= (-4.346-fM) ) iy = 0;
911 if ( track_coordy[(2*0)+1][0] >= ( -4.104+fM) && track_coordy[(2*0)+1][0] <= ( 3.704-fM) ) iy = 1;
912 if ( track_coordy[(2*0)+1][0] >= ( 3.946+fM) && track_coordy[(2*0)+1][0] <= (11.754-fM) ) iy = 2;
913 if ( cl2->cbar[(2*0)+1][0] >= (ye1 + fM1) && cl2->cbar[(2*0)+1][0] <= (ye2 - fM1) ) ix = 0;
914 if ( cl2->cbar[(2*0)+1][0] >= (ye3 + fM1) && cl2->cbar[(2*0)+1][0] <= (ye4 - fM1) ) ix = 1;
915 if ( cl2->cbar[(2*0)+1][0] >= (ye5 + fM1) && cl2->cbar[(2*0)+1][0] <= (ye6 - fM1) ) ix = 2;
916 if ( ix > -1 && iy > -1 ) fColYE = ix + iy*3;
917 //
918 // check event is inside XO acceptance
919 //
920 for (Int_t i=0; i<11; i++) {
921 if ((
922 ( track_coordy[(2*i)+1][0] >= (-12.154+fM) && track_coordy[(2*i)+1][0] <= (-4.346-fM) && (iy == 0 || multicol) ) ||
923 ( track_coordy[(2*i)+1][0] >= ( -4.104+fM) && track_coordy[(2*i)+1][0] <= ( 3.704-fM) && (iy == 1 || multicol) ) ||
924 ( track_coordy[(2*i)+1][0] >= ( 3.946+fM) && track_coordy[(2*i)+1][0] <= (11.754-fM) && (iy == 2 || multicol) )
925 ) && (
926 ( cl2->cbar[(2*i)+1][0] >= (ye1 + fM1) && cl2->cbar[(2*i)+1][0] <= (ye2 - fM1) && (ix == 0 || multicol) ) ||
927 ( cl2->cbar[(2*i)+1][0] >= (ye3 + fM1) && cl2->cbar[(2*i)+1][0] <= (ye4 - fM1) && (ix == 1 || multicol) ) ||
928 ( cl2->cbar[(2*i)+1][0] >= (ye5 + fM1) && cl2->cbar[(2*i)+1][0] <= (ye6 - fM1) && (ix == 2 || multicol) )
929 )){
930 fYesel = true;
931 fYeout = i;
932 } else {
933 fYesel = false;
934 break;
935 };
936 };
937 //
938 if ( !fYesel && fYeout >= fYemin && fYeout >= (Int_t)(fMax_planeye+nPl) ){
939 if ( debug ) printf(" YE: this event is only partially contained: fYeout %i fYemin %i fMax_planeye + nPl %i \n",fYeout,fYemin,(Int_t)(fMax_planeye+nPl));
940 fPartsel = true;
941 fYesel = true;
942 };
943 if ( fYesel ){
944 for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeye+nPl)) ;iplm++) fYEen_maxplane += en_yep[iplm];
945 fEnergyye = fYEen_maxplane/fConv_rye;
946 //
947 for (Int_t i=0;i<11;i++){
948 for(strip=0; strip<96; strip++) {
949 //
950 // run over all the strips of the plane
951 //
952 if ( strip >= 0 && strip < 32 ) encol[0][0] += enstrip[0][(2*i)+1][strip];
953 if ( strip >= 32 && strip < 64 ) encol[0][1] += enstrip[0][(2*i)+1][strip];
954 if ( strip >= 64 && strip < 96 ) encol[0][2] += enstrip[0][(2*i)+1][strip];
955 entot[0] += enstrip[0][(2*i)+1][strip];
956 //
957 };
958 };
959 //
960 };
961 };
962 //
963 if ( section.Contains("YO") ){
964 //
965 // find the column hit in the first plane
966 //
967 Int_t ix = -1;
968 Int_t iy = -1;
969 if ( track_coordy[(2*0)][0] >= (-11.954+fM) && track_coordy[(2*0)][0] <= (-4.146-fM) ) iy = 0;
970 if ( track_coordy[(2*0)][0] >= ( -3.904+fM) && track_coordy[(2*0)][0] <= ( 3.904-fM) ) iy = 1;
971 if ( track_coordy[(2*0)][0] >= ( 4.146+fM) && track_coordy[(2*0)][0] <= (11.954-fM) ) iy = 2;
972 if ( cl2->cbar[(2*0)][0] >= (yo1 + fM1) && cl2->cbar[(2*0)][0] <= (yo2 - fM1) ) ix = 0;
973 if ( cl2->cbar[(2*0)][0] >= (yo3 + fM1) && cl2->cbar[(2*0)][0] <= (yo4 - fM1) ) ix = 1;
974 if ( cl2->cbar[(2*0)][0] >= (yo5 + fM1) && cl2->cbar[(2*0)][0] <= (yo6 - fM1) ) ix = 2;
975 if ( ix > -1 && iy > -1 ) fColYO = ix + iy*3;
976 //
977 // check event is inside XO acceptance
978 //
979 for (Int_t i=0; i<11; i++) {
980 if ((
981 ( track_coordy[(2*i)][0] >= (-11.954+fM) && track_coordy[(2*i)][0] <= (-4.146-fM) && (iy == 0 || multicol) ) ||
982 ( track_coordy[(2*i)][0] >= ( -3.904+fM) && track_coordy[(2*i)][0] <= ( 3.904-fM) && (iy == 1 || multicol) ) ||
983 ( track_coordy[(2*i)][0] >= ( 4.146+fM) && track_coordy[(2*i)][0] <= (11.954-fM) && (iy == 2 || multicol) )
984 ) && (
985 ( cl2->cbar[(2*i)][0] >= (yo1 + fM1) && cl2->cbar[(2*i)][0] <= (yo2 - fM1) && (ix == 0 || multicol) ) ||
986 ( cl2->cbar[(2*i)][0] >= (yo3 + fM1) && cl2->cbar[(2*i)][0] <= (yo4 - fM1) && (ix == 1 || multicol) ) ||
987 ( cl2->cbar[(2*i)][0] >= (yo5 + fM1) && cl2->cbar[(2*i)][0] <= (yo6 - fM1) && (ix == 2 || multicol) )
988 )){
989 fYosel = true;
990 fYoout = i;
991 } else {
992 fYosel = false;
993 break;
994 };
995 };
996 //
997 if ( !fYosel && fYoout >= fYomin && fYoout >= (Int_t)(fMax_planeyo+nPl) ){
998 if ( debug ) printf(" YO: this event is only partially contained: fYoout %i fYomin %i fMax_planeyo + nPl %i \n",fYoout,fYomin,(Int_t)(fMax_planeyo+nPl));
999 fPartsel = true;
1000 fYosel = true;
1001 };
1002 if ( fYosel ){
1003 for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeyo+nPl)) ;iplm++) fYOen_maxplane += en_yop[iplm];
1004 fEnergyyo = fYOen_maxplane/fConv_ryo;
1005 //
1006 for (Int_t i=0;i<11;i++){
1007 for(strip=0; strip<96; strip++) {
1008 //
1009 // run over all the strips of the plane
1010 //
1011 if ( strip >= 0 && strip < 32 ) encol[0][0] += enstrip[0][(2*i)][strip];
1012 if ( strip >= 32 && strip < 64 ) encol[0][1] += enstrip[0][(2*i)][strip];
1013 if ( strip >= 64 && strip < 96 ) encol[0][2] += enstrip[0][(2*i)][strip];
1014 entot[0] += enstrip[0][(2*i)][strip];
1015 //
1016 };
1017 };
1018 };
1019 };
1020 //
1021 // Count the number of sections in which the event is contained
1022 //
1023 fCount = (Float_t)((Int_t)fXesel+(Int_t)fXosel+(Int_t)fYesel+(Int_t)fYosel);
1024 //
1025 if ( indep ){
1026 //
1027 // independent mode, average the energy measurement and max plane of the contained sections
1028 //
1029 fSel = ( fXesel || fYesel || fXosel || fYosel );
1030 fMax_plane = (Float_t)(fMax_planeyo+fMax_planeye+fMax_planexo+fMax_planexe)/fCount;
1031 fEnergy = (fEnergyxe+fEnergyyo+fEnergyye+fEnergyxo)/fCount;
1032 //
1033 } else {
1034 //
1035 // coherent mode, sum the energy [MIP] of the given sections and convert using fConv_rxo. **** NB: it is assumed that the conversion factor is unique and the method SetConvertionFactor(Float_t) has been used**** The event is selected only if it is contained in all the given sections
1036 //
1037 if ( fCount != fNumSec ){
1038 fSel = false;
1039 } else {
1040 fSel = true;
1041 };
1042 fEnergy = (fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane)/fConv_rxo;
1043 if ( fSel ){
1044 if ( fXesel ) fColumn = fColXE;
1045 if ( fXosel ){
1046 if ( fColXO != fColumn && fColumn > -1 ){
1047 printf(" ERROR! mismatch in column number between different sections! fColumn %i fColXO %i \n",fColumn,fColXO);
1048 } else {
1049 fColumn = fColXO;
1050 };
1051 };
1052 if ( fYesel ){
1053 if ( fColYE != fColumn && fColumn > -1 ){
1054 printf(" ERROR! mismatch in column number between different sections! fColumn %i fColYE %i \n",fColumn,fColYE);
1055 } else {
1056 fColumn = fColYE;
1057 };
1058 };
1059 if ( fYosel ){
1060 if ( fColYO != fColumn && fColumn > -1 ){
1061 printf(" ERROR! mismatch in column number between different sections! fColumn %i fColYO %i \n",fColumn,fColYO);
1062 } else {
1063 fColumn = fColYO;
1064 };
1065 };
1066 };
1067 };
1068 //
1069 if ( debug ) printf("sel %i indep %i fMax_plane %f conv_r %f en_maxplane %f encalo %f \n",fSel,indep,fMax_plane,fConv_rxo,fXOen_maxplane,fEnergy);
1070 if ( debug ) printf(" IsInside XE %i XO %i YE %i YO %i => SEL %i \n",fXesel,fXosel,fYesel,fYosel,fSel);
1071 //
1072 // finish exit
1073 //
1074 return fSel;
1075 //
1076 }
1077
1078 void CaloEnergy::Process(){
1079 TString xo = "XO";
1080 this->Process(xo);
1081 }
1082
1083
1084 void CaloEnergy::Process(TString section){
1085 //
1086 // process the event
1087 //
1088 TString ntr = section;
1089 if ( !L2 ){
1090 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
1091 printf(" ERROR: CaloEnergy variables not filled \n");
1092 return;
1093 };
1094 //
1095 Bool_t newentry = false;
1096 //
1097 if ( L2->IsORB() ){
1098 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || strcmp(ntr.Data(),sntr.Data()) ){
1099 newentry = true;
1100 OBT = L2->GetOrbitalInfo()->OBT;
1101 PKT = L2->GetOrbitalInfo()->pkt_num;
1102 atime = L2->GetOrbitalInfo()->absTime;
1103 sntr = ntr;
1104 };
1105 } else {
1106 newentry = true;
1107 };
1108 //
1109 // if we have already called this method for this event and no input changed then return fSel and exit
1110 //
1111 if ( !newentry ) return;
1112 //
1113 // process the event
1114 //
1115 if ( debug ) printf(" Processing event at OBT %u PKT %u time %u section %s\n",OBT,PKT,atime,section.Data());
1116 //
1117 // check if the cylinder of integration can go out of the sensor given the frame which has been set (if we use all the calorimeter fRad is < 0 and the printout is suppressed)
1118 //
1119 if ( (fM1+0.122-0.244*(Float_t)fRad) < 0. ) printf("Error: (fM1+0.122-0.244*(Float_t)fRad) < 0. fM1 %f fRad %i %f \n",fM1,fRad,(fM1+0.122-0.244*(Float_t)fRad));
1120 //
1121 if ( fLong ){
1122 if ( debug ) printf(" ==================================================================> LONGITUDINAL FIT! \n");
1123 //
1124 // use long fit to measure energy
1125 //
1126 if ( this->IsInsideAcceptance(section) ){
1127 //
1128 if ( debug ) printf(" ==================================================================> LONG INSIDE! \n");
1129 //
1130 Float_t myene[2][22];
1131 memset(myene,0,(sizeof(Float_t))*2*22);
1132 for (Int_t j=0; j<11; j++){
1133 if ( section.Contains("XE") ) myene[1][2*j] = en_xep[j];
1134 if ( section.Contains("YO") ) myene[0][2*j] = en_yop[j];
1135 if ( section.Contains("XO") ) myene[1][(2*j)+1] = en_xop[j];
1136 if ( section.Contains("YE") ) myene[0][(2*j)+1] = en_yep[j];
1137 };
1138 clong->UnMaskSections();
1139 if ( !(section.Contains("YE")) ) clong->MaskSection("YE");
1140 if ( !(section.Contains("YO")) ) clong->MaskSection("YO");
1141 if ( !(section.Contains("XO")) ) clong->MaskSection("XO");
1142 if ( !(section.Contains("XE")) ) clong->MaskSection("XE");
1143 clong->ForceNextFit();
1144 clong->SetEnergies(myene);
1145 if ( debug ){
1146 clong->Fit(true);
1147 } else {
1148 clong->Fit();
1149 };
1150 if ( clong->GetLowerLimit() != 0. || clong->GetUpperLimit() != 0. ){
1151 fXOen_maxplane = clong->Get_defE0();
1152 } else {
1153 fXOen_maxplane = clong->Get_E0();
1154 };
1155 fMax_plane = clong->Get_tmax();
1156 fYOen_maxplane = 0.;
1157 fYEen_maxplane = 0.;
1158 fXEen_maxplane = 0.;
1159 fEnergy=fXOen_maxplane/fConv_rxo;
1160 if ( fEnergy != fEnergy || clong->Get_fitresult() != 0 ) fEnergy = -1.;
1161 // if ( fEnergy != fEnergy ) fEnergy = 1.;
1162 //
1163 } else {
1164 //
1165 // if the event is not in the acceptance, return a negative energy.
1166 //
1167 if ( debug ) printf(" Outside acceptance \n");
1168 fEnergy *= -1.;
1169 //
1170 };
1171 //
1172 } else {
1173 //
1174 // use the energy measurement
1175 //
1176 if ( this->IsInsideAcceptance(section) ){
1177 //
1178 // the event is good
1179 //
1180 if ( debug ) printf(" XE %i XO %i YE %i YO %i \n",fXesel,fXosel,fYesel,fYosel);
1181 //
1182 } else {
1183 //
1184 // if the event is not in the acceptance, return a negative energy.
1185 //
1186 if ( debug ) printf(" Outside acceptance \n");
1187 fEnergy *= -1.;
1188 //
1189 };
1190 };
1191 //
1192 }

  ViewVC Help
Powered by ViewVC 1.1.23