/[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.7 - (show annotations) (download)
Wed Aug 12 14:44:10 2009 UTC (15 years, 4 months ago) by mocchiut
Branch: MAIN
Changes since 1.6: +23 -13 lines
New small features added

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

  ViewVC Help
Powered by ViewVC 1.1.23