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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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