/[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.11 - (show annotations) (download)
Mon Aug 17 15:00:53 2009 UTC (15 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.10: +18 -12 lines
New features added + bugs fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23