/[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.19 - (show annotations) (download)
Thu Dec 3 10:43:19 2009 UTC (15 years, 4 months ago) by mocchiut
Branch: MAIN
Changes since 1.18: +1 -1 lines
Relaxed condition on error printout in case of low margin left

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

  ViewVC Help
Powered by ViewVC 1.1.23