/[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.20 - (hide annotations) (download)
Mon Dec 14 14:56:37 2009 UTC (15 years ago) by mocchiut
Branch: MAIN
Changes since 1.19: +4 -0 lines
Do not use plane 18x by default

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

  ViewVC Help
Powered by ViewVC 1.1.23