/[PAMELA software]/calo/flight/CaloEnergy/src/CaloEnergy.cpp
ViewVC logotype

Annotation of /calo/flight/CaloEnergy/src/CaloEnergy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.12 - (hide annotations) (download)
Tue Aug 18 09:24:51 2009 UTC (15 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.11: +5 -0 lines
_incomplete_

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

  ViewVC Help
Powered by ViewVC 1.1.23