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

  ViewVC Help
Powered by ViewVC 1.1.23