/[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.21 - (show annotations) (download)
Wed Jun 30 12:20:09 2010 UTC (14 years, 5 months ago) by mocchiut
Branch: MAIN
CVS Tags: HEAD
Changes since 1.20: +1 -1 lines
Warning suppressed

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

  ViewVC Help
Powered by ViewVC 1.1.23