--- calo/flight/CaloEnergy/src/CaloEnergy.cpp 2009/07/13 16:10:23 1.3 +++ calo/flight/CaloEnergy/src/CaloEnergy.cpp 2009/08/21 07:06:37 1.14 @@ -35,7 +35,51 @@ // } +void CaloEnergy::Delete(){ + Clear(); + delete this; +} + +void CaloEnergy::UseLevel2(){ + if ( clong ){ + delete clong; + clong = new CaloLong(L2); + clong->SplitInto(0,22); + }; + if ( cp ) delete cp; + cp = NULL; +} + +void CaloEnergy::UseCaloPreSampler(){ + // + // 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 + // + if ( !cp ) cp = new CaloPreSampler(L2); + cp->SplitInto(0,22); + cp->SetForceFitMode(1000); +// cp->UseTracker(false); +// cp->ForceCaloFit(); +// cp->SetDebug(true); +// cp->Process(); + if ( clong ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer()); +} + + +void CaloEnergy::UseLongFit(){ + fPl = 0; + fLong = true; + if ( !clong ){ + clong = new CaloLong(L2); + if ( cp ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer()); + clong->SplitInto(0,22); + }; + // +} + void CaloEnergy::Set(){ + // + // set default values, NB default conversion factor for energy is just very approximated! + // OBT = 0; PKT = 0; atime = 0; @@ -48,9 +92,17 @@ // debug = false; // + indep = false; + // + fAllpl = true; fLong = false; fPl = 1; fRad = -1; + cp = NULL; + clong = NULL; + x0max = -1.; + // + multicol = false; // this->DefineGeometry(); fXosel =true; @@ -65,147 +117,99 @@ fXemin = 1000; fYomin = 1000; fYemin = 1000; + // + this->UseCaloPreSampler(); // use it by default, to go back to "standard" mode use CaloEnergy::UseLevel2(). + // } - -void CaloEnergy::SetMinimumContainment(TString section, Int_t plane){ - section.ToUpper(); - if ( section.Contains("XO") ) fXomin = plane; - if ( section.Contains("XE") ) fXemin = plane; - if ( section.Contains("YO") ) fYomin = plane; - if ( section.Contains("YE") ) fYemin = plane; -} - -Int_t CaloEnergy::GetMinimumContainment(TString section){ - section.ToUpper(); - if ( section.Contains("XO") ) return(fXomin); - if ( section.Contains("XE") ) return(fXemin); - if ( section.Contains("YE") ) return(fYemin); - if ( section.Contains("YO") ) return(fYomin); - printf(" ERROR: section not recognized \n"); - return(-1000); -} - -void CaloEnergy::SetConversionFactor(TString section, Float_t conv){ - section.ToUpper(); - if ( section.Contains("XO") ) fConv_rxo = conv; - if ( section.Contains("XE") ) fConv_rxe = conv; - if ( section.Contains("YO") ) fConv_ryo = conv; - if ( section.Contains("YE") ) fConv_rye = conv; -} - -Float_t CaloEnergy::GetConversionFactor(TString section){ - section.ToUpper(); - if ( section.Contains("XO") ) return(fConv_rxo); - if ( section.Contains("XE") ) return(fConv_rxe); - if ( section.Contains("YO") ) return(fConv_ryo); - if ( section.Contains("YE") ) return(fConv_rye); - printf(" ERROR: section not recognized \n"); - return(-1000.); -} - -Int_t CaloEnergy::GetMaxplane(TString section){ - section.ToUpper(); - if ( section.Contains("XO") ) return fMax_planexo; - if ( section.Contains("XE") ) return fMax_planexe; - if ( section.Contains("YO") ) return fMax_planeyo; - if ( section.Contains("YE") ) return fMax_planeye; - return(-1); -} - void CaloEnergy::DefineGeometry(){ // - fM = 2. + 0.096; // real position from cbar - fM1 = 2. - 0.122 - 0.096; // due to calculation of xe1 etc. + // Use CaloStrip to determine once the position of border strips for each section + // + // fM = 2. + 0.096; // real position from cbar BUG the 0.096 is already taken into account in the border calculation made by Giovanna + fM = 2. ; // real position from cbar + // 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 + fM1 = 2. - 0.122 - 0.096 + 0.096; // due to calculation of xe1 etc. if ( fM1 < 0. ) fM1 = 0.; // CaloStrip *cs = new CaloStrip(fSimu); - - // vista y piano 0 striscia 0 + // + // view y plane 0 strip 0 cs->Set(1,0,0); xe1= cs->GetY(); - z1= cs->GetZ(); - // vista y piano 0 striscia 31 + // view y plane 0 strip 31 cs->Set(1,0,31); xe2= cs->GetY(); - // vista y piano 0 striscia 32 + // view y plane 0 strip 32 cs->Set(1,0,32); xe3= cs->GetY(); - // vista y piano 0 striscia 63 + // view y plane 0 strip 63 cs->Set(1,0,63); xe4= cs->GetY(); - // vista y piano 0 striscia 64 + // view y plane 0 strip 64 cs->Set(1,0,64); xe5= cs->GetY(); - // vista y piano 0 striscia 95 + // view y plane 0 strip 95 cs->Set(1,0,95); xe6= cs->GetY(); - - // vista x piano 0 striscia 0 + // view x plane 0 strip 0 cs->Set(0,0,0); yo1= cs->GetX(); - z2= cs->GetZ(); - // vista x piano 0 striscia 31 + // view x plane 0 strip 31 cs->Set(0,0,31); yo2= cs->GetX(); - // vista x piano 0 striscia 32 + // view x plane 0 strip 32 cs->Set(0,0,32); yo3= cs->GetX(); - // vista x piano 0 striscia 63 + // view x plane 0 strip 63 cs->Set(0,0,63); yo4= cs->GetX(); - // vista x piano 0 striscia 64 + // view x plane 0 strip 64 cs->Set(0,0,64); yo5= cs->GetX(); - // vista x piano 0 striscia 95 + // view x plane 0 strip 95 cs->Set(0,0,95); yo6= cs->GetX(); - - - // vista y piano 1 striscia 0 + // view y plane 1 strip 0 cs->Set(1,1,0); xo1= cs->GetY(); - z3= cs->GetZ(); - // vista y piano 1 striscia 31 + // view y plane 1 strip 31 cs->Set(1,1,31); xo2= cs->GetY(); - // vista y piano 1 striscia 32 + // view y plane 1 strip 32 cs->Set(1,1,32); xo3= cs->GetY(); - // vista y piano 1 striscia 63 + // view y plane 1 strip 63 cs->Set(1,1,63); xo4= cs->GetY(); - // vista y piano 1 striscia 64 + // view y plane 1 strip 64 cs->Set(1,1,64); xo5= cs->GetY(); - // vista y piano 1 striscia 95 + // view y plane 1 strip 95 cs->Set(1,1,95); - xo6= cs->GetY(); - - // vista x piano 1 striscia 0 + xo6= cs->GetY(); + // view x plane 1 strip 0 cs->Set(0,1,0); ye1= cs->GetX(); - z4= cs->GetZ(); - // vista x piano 1 striscia 31 + // view x plane 1 strip 31 cs->Set(0,1,31); ye2= cs->GetX(); - // vista x piano 1 striscia 32 + // view x plane 1 strip 32 cs->Set(0,1,32); ye3= cs->GetX(); - // vista x piano 1 striscia 63 + // view x plane 1 strip 63 cs->Set(0,1,63); ye4= cs->GetX(); - // vista x piano 1 striscia 64 + // view x plane 1 strip 64 cs->Set(0,1,64); ye5= cs->GetX(); - // vista x piano 1 striscia 95 + // view x plane 1 strip 95 cs->Set(0,1,95); ye6= cs->GetX(); - // for (Int_t p = 0; p<22; p ++){ for (Int_t v = 0; v<2; v++ ){ cs->Set(v,p,0); - trk_z[p][v]= cs->GetZ(); // << cooordinata Z piano + trk_z[p][v]= cs->GetZ(); // Z coord for each plane }; }; // @@ -215,6 +219,8 @@ void CaloEnergy::Clear(){ // + // clear variables + // fPartsel = false; fSel = false; fXosel = false; @@ -231,7 +237,12 @@ fMax_planexo = 0; fMax_planexe = 0; fMax_planeyo = 0; - fMax_planeye = 0; + fMax_planeye = 0; + xomax_en= 0.; + xemax_en= 0.; + yomax_en= 0.; + yemax_en= 0.; + // memset(enstrip,0,2*22*96*(sizeof(Float_t))); en = 0.; view = 0; @@ -245,11 +256,24 @@ fYeout = 0; fXoout = 0; fXeout = 0; + fXEen_maxplane = 0.; + fXOen_maxplane = 0.; + fYEen_maxplane = 0.; + fYOen_maxplane = 0.; memset(en_xep,0,11*sizeof(Float_t)); memset(en_yep,0,11*sizeof(Float_t)); memset(en_xop,0,11*sizeof(Float_t)); memset(en_yop,0,11*sizeof(Float_t)); // + fColumn = -1; + fColXE = -1; + fColXO = -1; + fColYE = -1; + fColYO = -1; + memset(encol,0,2*3*sizeof(Float_t)); + entot[0] = 0.; + entot[1] = 0.; + // } void CaloEnergy::Print(){ @@ -282,6 +306,16 @@ printf(" fM1 :.............. %f \n",fM1); printf(" fRad :.............. %i \n",fRad); printf(" fPl :.............. %i \n",fPl); + printf(" fColumn :.............. %i \n",fColumn); + printf(" multicol:.............. %i \n",multicol); + printf(" encol x :.............. %f \n",this->GetEncol(0)); + printf(" encol y :.............. %f \n",this->GetEncol(1)); + printf(" entot x :.............. %f \n",this->GetEntot(0)); + printf(" entot y :.............. %f \n",this->GetEntot(1)); + printf(" fColXE :.............. %i \n",fColXE); + printf(" fColXO :.............. %i \n",fColXO); + printf(" fColYE :.............. %i \n",fColYE); + printf(" fColYO :.............. %i \n",fColYO); printf(" fConv_rxe ............. %f \n",fConv_rxe); printf(" fConv_rxo ............. %f \n",fConv_rxo); printf(" fConv_ryo ............. %f \n",fConv_ryo); @@ -291,19 +325,147 @@ printf(" energyxo:.............. %f \n",energyxo); printf(" energyye:.............. %f \n",energyye); printf(" energyyo:.............. %f \n",energyyo); + printf(" fXEen_maxplane:........ %f \n",fXEen_maxplane); + printf(" fXOen_maxplane:........ %f \n",fXOen_maxplane); + printf(" fYEen_maxplane:........ %f \n",fYEen_maxplane); + printf(" fYOen_maxplane:........ %f \n",fYOen_maxplane); + printf(" x0max :.............. %f \n",x0max); printf(" debug :.............. %i \n",debug); printf("========================================================================\n"); // } -void CaloEnergy::Delete(){ - Clear(); +void CaloEnergy::SetMinimumContainment(TString section, Int_t plane){ + section.ToUpper(); + if ( section.Contains("XO") ) fXomin = plane; + if ( section.Contains("XE") ) fXemin = plane; + if ( section.Contains("YO") ) fYomin = plane; + if ( section.Contains("YE") ) fYemin = plane; +} + +void CaloEnergy::SetMinimumContainment(Int_t plane){ + this->SetMinimumContainment("XEXOYEYO",plane); +} + +void CaloEnergy::SetConversionFactor(TString section, Float_t conv){ + section.ToUpper(); + if ( section.Contains("XO") ) fConv_rxo = conv; + if ( section.Contains("XE") ) fConv_rxe = conv; + if ( section.Contains("YO") ) fConv_ryo = conv; + if ( section.Contains("YE") ) fConv_rye = conv; +} + +void CaloEnergy::SetConversionFactor(Float_t conv){ + this->SetConversionFactor("XEXOYEYO",conv); +} + +Int_t CaloEnergy::GetMinimumContainment(TString section){ + section.ToUpper(); + if ( section.Contains("XO") ) return(fXomin); + if ( section.Contains("XE") ) return(fXemin); + if ( section.Contains("YE") ) return(fYemin); + if ( section.Contains("YO") ) return(fYomin); + printf(" ERROR: section not recognized \n"); + return(-1000); +} + +Float_t CaloEnergy::GetConversionFactor(TString section){ + section.ToUpper(); + if ( section.Contains("XO") ) return(fConv_rxo); + if ( section.Contains("XE") ) return(fConv_rxe); + if ( section.Contains("YO") ) return(fConv_ryo); + if ( section.Contains("YE") ) return(fConv_rye); + printf(" ERROR: section not recognized \n"); + return(-1000.); +} + +Int_t CaloEnergy::GetMaxplane(TString section){ + section.ToUpper(); + if ( section.Contains("XO") ) return fMax_planexo; + if ( section.Contains("XE") ) return fMax_planexe; + if ( section.Contains("YO") ) return fMax_planeyo; + if ( section.Contains("YE") ) return fMax_planeye; + return(-1); +} + +Int_t CaloEnergy::GetColumn(TString section){ + section.ToUpper(); + if ( section.Contains("XO") ) return fColXO; + if ( section.Contains("XE") ) return fColXE; + if ( section.Contains("YO") ) return fColYO; + if ( section.Contains("YE") ) return fColYE; + return(-1); +} + +Float_t CaloEnergy::GetMipEnergyAtMaxplane(TString section){ + printf(" WARNING: OBSOLETE METHOD, use GetMipEnergyAtMaxplane(TString) instead! \n"); + return (this->GetEnergyAtMaxplane(section)); +} + +Float_t CaloEnergy::GetEnergyAtMaxplane(TString section){ + section.ToUpper(); + if ( section.Contains("XO") ) return xomax_en; + if ( section.Contains("XE") ) return xemax_en; + if ( section.Contains("YO") ) return yomax_en; + if ( section.Contains("YE") ) return yemax_en; + return(-1); +} + +Float_t CaloEnergy::GetMaxEnergy(TString section){ + printf(" WARNING: OBSOLETE METHOD, use GetMipEnergy(TString) instead! \n"); + return (this->GetMipEnergy(section)); +} + +Float_t CaloEnergy::GetMipEnergy(TString section){ + section.ToUpper(); + if ( fLong ){ + this->Process(section); + return fXOen_maxplane; + } else { + if ( section.Contains("XO") ) return fXOen_maxplane; + if ( section.Contains("XE") ) return fXEen_maxplane; + if ( section.Contains("YO") ) return fYOen_maxplane; + if ( section.Contains("YE") ) return fYEen_maxplane; + }; + return(-1); +} + +Float_t CaloEnergy::GetEncol(Int_t i){ + if ( fColumn > -1 && (((fXesel || fXosel)&&i==1) || ((fYesel || fYosel)&&i==0)) ){ + Int_t t = -1; + if ( i == 0 ){ + if ( fColumn == 0 || fColumn == 3 || fColumn == 6 ) t = 0; + if ( fColumn == 1 || fColumn == 4 || fColumn == 7 ) t = 1; + if ( fColumn == 2 || fColumn == 5 || fColumn == 8 ) t = 2; + } else { + if ( fColumn == 0 || fColumn == 1 || fColumn == 2 ) t = 0; + if ( fColumn == 3 || fColumn == 4 || fColumn == 5 ) t = 1; + if ( fColumn == 6 || fColumn == 7 || fColumn == 8 ) t = 2; + }; + if ( debug ) printf(" encol: i %i t %i encol %f \n",i,t,encol[i][t]); + return encol[i][t]; + }; + return(-1.); +} + +Float_t CaloEnergy::GetMaxEnergy(){ + printf(" WARNING: OBSOLETE METHOD, use GetMipEnergy() instead! \n"); + return (this->GetMipEnergy()); +} + +Float_t CaloEnergy::GetMipEnergy(){ + if ( fLong ){ + if ( debug ) printf(" oh! call process! with asntr %s and sntr %s \n",asntr.Data(),sntr.Data()); + this->Process(asntr); + }; + return((fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane)); } Bool_t CaloEnergy::IsInsideAcceptance(TString section){ // + // check if the event is inside the acceptance of the given section(s) // TString ntr = section; if ( !L2 ){ @@ -326,127 +488,320 @@ newentry = true; }; // + // if we have already called this method for this event and no input changed then return fSel and exit + // if ( !newentry ) return fSel; // + // process the event + // if ( debug ) printf(" ########## IsInsideAcceptance ######### \n"); // + // clear variables + // this->Clear(); // section.ToUpper(); // - Int_t ss = Int_t(section.Contains("XO"))+Int_t(section.Contains("XE"))+Int_t(section.Contains("YO"))+Int_t(section.Contains("YE")); - if ( !ss ){ + // Count the number of section(s) given as input + // + Int_t fNumSec = Int_t(section.Contains("XO"))+Int_t(section.Contains("XE"))+Int_t(section.Contains("YO"))+Int_t(section.Contains("YE")); + if ( !fNumSec ){ printf(" ERROR: section must be XO or XE or YO or YE while it is %s \n",section.Data()); return false; }; // + // If the presampler object exists then use the presampler output instead of the level2 output + // + CaloLevel2 *cl2 = NULL; + if ( cp ){ + cl2 = cp->GetCaloLevel2(); + } else { + cl2 = L2->GetCaloLevel2(); + }; + // + // get the energy for every strip of the calorimeter // for (Int_t ch=0; ch< L2->GetCaloLevel1()->istrip; ch++){ en = L2->GetCaloLevel1()->DecodeEstrip(ch,view,plane,strip); enstrip[view][plane][strip]=en; }; + // + 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)); + // + // sum energy plane by plane for each sections + // + Float_t fen_xep[11]; + Float_t fen_xop[11]; + Float_t fen_yep[11]; + Float_t fen_yop[11]; + memset(fen_xep,0,11*sizeof(Float_t)); + memset(fen_xop,0,11*sizeof(Float_t)); + memset(fen_yep,0,11*sizeof(Float_t)); + memset(fen_yop,0,11*sizeof(Float_t)); + // for (Int_t i=0;i<11;i++){ for(strip=0; strip<96; strip++) { + fen_xep[i] += enstrip[1][2*i][strip]; + fen_yop[i] += enstrip[0][2*i][strip]; + fen_xop[i] += enstrip[1][(2*i)+1][strip]; + fen_yep[i] += enstrip[0][(2*i)+1][strip]; if ( fRad < 0 ){ + // + // run over all the strips of the plane + // en_xep[i] += enstrip[1][2*i][strip]; en_yop[i] += enstrip[0][2*i][strip]; - en_xop[i] += enstrip[1][2*i+1][strip]; - en_yep[i] += enstrip[0][2*i+1][strip]; + 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]; + // + // use only the strips inside a cylinder of given radius fRad + // + 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]; + + 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]; + + 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]; + + 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]; }; }; + if ( debug ) printf(" ex_xep[%i] %f cibar %i \n",i,en_xep[i],cl2->cibar[2*i][1]); + if ( debug ) printf(" ex_xop[%i] %f cibar %i \n",i,en_xop[i],cl2->cibar[(2*i)+1][1]); + if ( debug ) printf(" ex_yep[%i] %f cibar %i \n",i,en_yep[i],cl2->cibar[(2*i)+1][0]); + if ( debug ) printf(" ex_yop[%i] %f cibar %i \n",i,en_yop[i],cl2->cibar[2*i][0]); energyxe += en_xep[i]; energyyo += en_yop[i]; energyxo += en_xop[i]; energyye += en_yep[i]; }; - Float_t en_maxplane = 0.; - Float_t 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; + // Find the plane of maximum for each section + // + // + Int_t xen = 0; + Int_t yon = 0; + Int_t xon = 0; + Int_t yen = 0; + Float_t en = 0.; + // + if ( section.Contains("XE") ){ + yon++; + xon++; + yen++; + for (Int_t ipl =0; ipl < 11; ipl ++) { + en = fen_xep[ipl]; + if ( !fAllpl ) en = en_xep[ipl]; + if(en > xemax_en) { + xemax_en = en; + fMax_planexe = ipl; }; + }; }; // - 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; + if ( section.Contains("YO") ){ + xon++; + yen++; + for (Int_t ipl =0; ipl < 11; ipl ++) { + en = fen_yop[ipl]; + if ( !fAllpl ) en = en_yop[ipl]; + if(en > yomax_en) { + yomax_en = en; + fMax_planeyo = ipl; + }; }; }; // - 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; + if ( section.Contains("XO") ){ + yen++; + for (Int_t ipl =0; ipl < 11; ipl ++) { + en = fen_xop[ipl]; + if ( !fAllpl ) en = en_xop[ipl]; + if(en > xomax_en) { + xomax_en = en; + fMax_planexo = ipl; + }; }; }; // - max_en = 0.; - for (Int_t ipl =0; ipl < 11; ipl ++) { - if(en_yep[ipl] > max_en) { - max_en = en_yep[ipl]; - fMax_planeye = ipl; + if ( section.Contains("YE") ){ + for (Int_t ipl =0; ipl < 11; ipl ++) { + en = fen_yep[ipl]; + if ( !fAllpl ) en = en_yep[ipl]; + if(en > yemax_en) { + yemax_en = en; + fMax_planeye = ipl; + }; }; }; // - Float_t track_coordx[22][2]; - Float_t track_coordy[22][2]; + // the maximum is given externally as X0, convert it to plane and section + // + if ( x0max > 0. ){ + if ( debug ) printf(" CALCULATE MAX PLANE GIVEN X0 ASSUMING PERPENDICULAR TRACK \n"); + Int_t wpl = (Int_t)roundf(x0max/0.76); + Bool_t isY = false; + if ( ((x0max/0.76)-(Float_t)wpl) > 0. ) isY = true; + xomax_en = 0.; + yemax_en = 0.; + xemax_en = 0.; + yomax_en = 0.; + // + if ( !(wpl%2) ){ + // 0, 2, 4, ... + if ( isY ){ + if ( section.Contains("YO") ) yomax_en = 1000.; + if ( section.Contains("XE") ) xemax_en = 500.; + fMax_planeyo=wpl/2; + fMax_planexe=wpl/2; + if ( section.Contains("XO") ) xomax_en = 10.; + if ( section.Contains("YE") ) yemax_en = 5.; + } else { + if ( section.Contains("YO") ) yomax_en = 500.; + if ( section.Contains("XE") ) xemax_en = 1000.; + fMax_planeyo=wpl/2; + fMax_planexe=wpl/2; + if ( section.Contains("XO") ) xomax_en = 5.; + if ( section.Contains("YE") ) yemax_en = 10.; + }; + } else { + // 1, 3, 5, ... + if ( isY ){ + if ( section.Contains("YE") ) yemax_en = 1000.; + if ( section.Contains("XO") ) xomax_en = 500.; + fMax_planeye=(wpl-1)/2; + fMax_planexo=(wpl-1)/2; + if ( section.Contains("XE") ) xemax_en = 10.; + if ( section.Contains("YO") ) yomax_en = 5.; + } else { + if ( section.Contains("YE") ) yemax_en = 500.; + if ( section.Contains("XO") ) xomax_en = 1000.; + fMax_planeye=(wpl-1)/2; + fMax_planexo=(wpl-1)/2; + if ( section.Contains("XE") ) xemax_en = 5.; + if ( section.Contains("YO") ) yomax_en = 10.; + }; + }; + 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); + }; + // + Int_t nPl = fPl; + // + // Set the maximum in case of coherent mode was selected + // + if ( !indep ){ + nPl = 0; + 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); + Int_t nummod = 0; + Int_t numexpl = 0; + if ( xomax_en > xemax_en && xomax_en > yemax_en && xomax_en > yomax_en ){ + // + // Section XO contains the maximum energy release per plane of the whole calorimeter + // + if ( debug ) printf(" XO is MAX %i %i %i %i\n",xen,yon,xon,yen); + // + // fMax_plane is the plane of maximum + number of additional dE/dx measurement counting planes from 0 to 43 + // + fMax_plane = (fNumSec * fMax_planexo) +(Float_t)xon + fPl; + // + // nummod is the integer part of the number of modules in which the maximum is contained + // + nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec); + // + // numexpl is the number of additional planes (0,1,2) inside the module + // + numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod); + // + }; + if ( xemax_en > xomax_en && xemax_en > yemax_en && xemax_en > yomax_en ){ + if ( debug ) printf(" XE is MAX %i %i %i %i\n",xen,yon,xon,yen); + fMax_plane = (fNumSec * fMax_planexe) +(Float_t)xen + fPl; + nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec); + numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod); + // + }; + + if ( yemax_en > xomax_en && yemax_en > xemax_en && yemax_en > yomax_en ){ + if ( debug ) printf(" YE is MAX %i %i %i %i\n",xen,yon,xon,yen); + fMax_plane = (fNumSec * fMax_planeye) +(Float_t)yen + fPl; + nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec); + numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod); + // + }; + if ( yomax_en > xemax_en && yomax_en > yemax_en && yomax_en > xomax_en ){ + if ( debug ) printf(" YO is MAX %i %i %i %i\n",xen,yon,xon,yen); + fMax_plane = (fNumSec * fMax_planeyo) +(Float_t)yon + fPl; + nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec); + numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod); + // + }; + // + // find the plane up to which is necessary to integrate the energy for each section + // + Int_t a = 0; + Int_t b = 0; + Int_t c = 0; + if ( numexpl > xen ) a = 1; + if ( numexpl > yon ) b = 1; + if ( numexpl > xon ) c = 1; + fMax_planexe = nummod; + fMax_planeyo = nummod - 1 + a; + fMax_planexo = nummod - 1 + b; + fMax_planeye = nummod - 1 + c; + if ( debug ) printf(" fMax_plane %f nummod %i numexpl %i a %i b %i c %i \n",fMax_plane,nummod,numexpl,a,b,c); + 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); + }; + // + // 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 + // // Float_t tgx_cl2; Float_t tgy_cl2; - tgx_cl2 = L2->GetCaloLevel2()->tanx[0]; - tgy_cl2 = L2->GetCaloLevel2()->tany[0]; + tgx_cl2 = cl2->tanx[0]; + tgy_cl2 = cl2->tany[0]; // for (Int_t p=0; p<22; p++){ - track_coordy[p][1] = L2->GetCaloLevel2()->cbar[p][1]; - track_coordx[p][1] = L2->GetCaloLevel2()->cbar[p][0] - fabs(trk_z[p][1]-trk_z[p][0])*tgx_cl2; - track_coordx[p][0] = L2->GetCaloLevel2()->cbar[p][0]; - track_coordy[p][0] = L2->GetCaloLevel2()->cbar[p][1] - fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2; + track_coordy[p][1] = cl2->cbar[p][1]; + track_coordx[p][1] = cl2->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; + track_coordx[p][0] = cl2->cbar[p][0]; + track_coordy[p][0] = cl2->cbar[p][1] - fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2; + // track_coordy[p][0] = cl2->cbar[p][1] + fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2; 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(" acceptance ss %i tgx %f tgy %f\n",ss,tgx_cl2,tgy_cl2); + if ( debug ) printf(" acceptance fNumSec %i tgx %f tgy %f\n",fNumSec,tgx_cl2,tgy_cl2); // if ( section.Contains("XO") ){ + // + // find the column hit in the first plane + // + Int_t ix = -1; + Int_t iy = -1; + if ( track_coordx[(2*0)+1][1] >= (-12.054+fM) && track_coordx[(2*0)+1][1] <= (-4.246-fM) ) ix = 0; + if ( track_coordx[(2*0)+1][1] >= ( -4.004+fM) && track_coordx[(2*0)+1][1] <= ( 3.804-fM) ) ix = 1; + if ( track_coordx[(2*0)+1][1] >= ( 4.046+fM) && track_coordx[(2*0)+1][1] <= (11.854-fM) ) ix = 2; + if ( cl2->cbar[(2*0)+1][1] >= (xo1 + fM1) && cl2->cbar[(2*0)+1][1] <= (xo2 - fM1) ) iy = 0; + if ( cl2->cbar[(2*0)+1][1] >= (xo3 + fM1) && cl2->cbar[(2*0)+1][1] <= (xo4 - fM1) ) iy = 1; + if ( cl2->cbar[(2*0)+1][1] >= (xo5 + fM1) && cl2->cbar[(2*0)+1][1] <= (xo6 - fM1) ) iy = 2; + if ( ix > -1 && iy > -1 ) fColXO = ix + iy*3; + // + // check event is inside XO acceptance, if multicol is false (SingleColumn mode) then the track must be contained in a column. + // for (Int_t i=0; i<11; i++) { - if - ( - (((track_coordx[2*i+1][1]>=(-12.054+fM))&& - (track_coordx[2*i+1][1]<=(-4.246-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo6 - fM1 ))))|| - - (((track_coordx[2*i+1][1]>=(-4.004+fM))&& - (track_coordx[2*i+1][1]<=(3.804-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo4 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo6 - fM1 ))))|| - - (((track_coordx[2*i+1][1]>=(4.046+fM))&& - (track_coordx[2*i+1][1]<=(11.854-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo2 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][1]>=xo5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][1]<=xo6 - fM1 )))) - ){ + if (( + ( track_coordx[(2*i)+1][1] >= (-12.054+fM) && track_coordx[(2*i)+1][1] <= (-4.246-fM) && (ix == 0 || multicol) ) || + ( track_coordx[(2*i)+1][1] >= ( -4.004+fM) && track_coordx[(2*i)+1][1] <= ( 3.804-fM) && (ix == 1 || multicol) ) || + ( track_coordx[(2*i)+1][1] >= ( 4.046+fM) && track_coordx[(2*i)+1][1] <= (11.854-fM) && (ix == 2 || multicol) ) + ) && ( + ( cl2->cbar[(2*i)+1][1] >= (xo1 + fM1) && cl2->cbar[(2*i)+1][1] <= (xo2 - fM1) && (iy == 0 || multicol) ) || + ( cl2->cbar[(2*i)+1][1] >= (xo3 + fM1) && cl2->cbar[(2*i)+1][1] <= (xo4 - fM1) && (iy == 1 || multicol) ) || + ( cl2->cbar[(2*i)+1][1] >= (xo5 + fM1) && cl2->cbar[(2*i)+1][1] <= (xo6 - fM1) && (iy == 2 || multicol) ) + )){ fXosel = true; fXoout = i; } else { @@ -454,49 +809,65 @@ break; }; }; - if ( !fXosel && fXoout >= fXomin && fXoout >= (Int_t)(fMax_planexo+fPl) ){ - if ( debug ) printf(" XO: this event is only partially contained: fXoout %i fXomin %i fMax_planexo + fPl %i \n",fXoout,fXomin,(Int_t)(fMax_planexo+fPl)); + // + // 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 + // + if ( !fXosel && fXoout >= fXomin && fXoout >= (Int_t)(fMax_planexo+nPl) ){ + 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)); fPartsel = true; fXosel = true; }; + // + // event is contained (or partially contained) hence we can integrate energy up to the maximum and calculate the energy as measured by this section + // if ( fXosel ){ - en_maxplane = 0.; - 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; + for (Int_t iplm=0; iplm<=TMath::Min(10,(Int_t)(fMax_planexo+nPl)); iplm++){ + fXOen_maxplane += en_xop[iplm]; + if ( debug ) printf(" XO iplm %i fXOen_maxplane %f en_xop[iplm] %f\n",iplm,fXOen_maxplane,en_xop[iplm]); + }; + fEnergyxo = fXOen_maxplane/fConv_rxo; + // + for (Int_t i=0;i<11;i++){ + for(strip=0; strip<96; strip++) { + // + // run over all the strips of the plane + // + if ( strip >= 0 && strip < 32 ) encol[1][0] += enstrip[1][(2*i)+1][strip]; + if ( strip >= 32 && strip < 64 ) encol[1][1] += enstrip[1][(2*i)+1][strip]; + if ( strip >= 64 && strip < 96 ) encol[1][2] += enstrip[1][(2*i)+1][strip]; + entot[1] += enstrip[1][(2*i)+1][strip]; + // + }; + }; }; }; // if ( section.Contains("XE") ){ + // + // find the column hit in the first plane + // + Int_t ix = -1; + Int_t iy = -1; + if ( track_coordx[(2*0)][1] >= (-11.854+fM) && track_coordx[(2*0)][1] <= (-4.046-fM) ) ix = 0; + if ( track_coordx[(2*0)][1] >= ( -3.804+fM) && track_coordx[(2*0)][1] <= ( 4.004-fM) ) ix = 1; + if ( track_coordx[(2*0)][1] >= ( 4.246+fM) && track_coordx[(2*0)][1] <= (12.054-fM) ) ix = 2; + if ( cl2->cbar[(2*0)][1] >= (xe1 + fM1) && cl2->cbar[(2*0)][1] <= (xe2 - fM1) ) iy = 0; + if ( cl2->cbar[(2*0)][1] >= (xe3 + fM1) && cl2->cbar[(2*0)][1] <= (xe4 - fM1) ) iy = 1; + if ( cl2->cbar[(2*0)][1] >= (xe5 + fM1) && cl2->cbar[(2*0)][1] <= (xe6 - fM1) ) iy = 2; + if ( ix > -1 && iy > -1 ) fColXE = ix + iy*3; + // + // check event is inside XO acceptance + // for (Int_t i=0; i<11; i++) { - if - ( - (((track_coordx[2*i][1]>=(-11.854+fM))&& - (track_coordx[2*i][1]<=(-4.046-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i][1]>=xe1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][1]>=xe3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][1]>=xe5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe6 - fM1 ))))|| - - (((track_coordx[2*i][1]>=(-3.804+fM))&& - (track_coordx[2*i][1]<=(4.004-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i][1]>=xe1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][1]>=xe3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe4 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i][1]>=xe5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe6 - fM1 ))))|| - - (((track_coordx[2*i][1]>=(4.246+fM))&& - (track_coordx[2*i][1]<=(12.054-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i][1]>=xe1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe2 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i][1]>=xe3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][1]>=xe5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][1]<=xe6 - fM1 )))) - ){ + if (( + ( track_coordx[(2*i)][1] >= (-11.854+fM) && track_coordx[(2*i)][1] <= (-4.046-fM) && (ix == 0 || multicol) ) || + ( track_coordx[(2*i)][1] >= ( -3.804+fM) && track_coordx[(2*i)][1] <= ( 4.004-fM) && (ix == 1 || multicol) ) || + ( track_coordx[(2*i)][1] >= ( 4.246+fM) && track_coordx[(2*i)][1] <= (12.054-fM) && (ix == 2 || multicol) ) + ) && ( + ( cl2->cbar[(2*i)][1] >= (xe1 + fM1) && cl2->cbar[(2*i)][1] <= (xe2 - fM1) && (iy == 0 || multicol) ) || + ( cl2->cbar[(2*i)][1] >= (xe3 + fM1) && cl2->cbar[(2*i)][1] <= (xe4 - fM1) && (iy == 1 || multicol) ) || + ( cl2->cbar[(2*i)][1] >= (xe5 + fM1) && cl2->cbar[(2*i)][1] <= (xe6 - fM1) && (iy == 2 || multicol) ) + )){ fXesel = true; fXeout = i; } else { @@ -504,49 +875,60 @@ break; }; }; - if ( !fXesel && fXeout >= fXemin && fXeout >= (Int_t)(fMax_planexe+fPl) ){ - if ( debug ) printf(" XE: this event is only partially contained: fXeout %i fXemin %i fMax_planexe + fPl %i \n",fXeout,fXemin,(Int_t)(fMax_planexe+fPl)); + // + if ( !fXesel && fXeout >= fXemin && fXeout >= (Int_t)(fMax_planexe+nPl) ){ + 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)); fPartsel = true; fXesel = true; }; if ( fXesel ){ - en_maxplane = 0.; - 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; + for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planexe+nPl)) ;iplm++){ + fXEen_maxplane += en_xep[iplm]; + if ( debug ) printf(" XE iplm %i fXOen_maxplane %f en_xop[iplm] %f\n",iplm,fXEen_maxplane,en_xep[iplm]); + }; + fEnergyxe = fXEen_maxplane/fConv_rxe; + // + for (Int_t i=0;i<11;i++){ + for(strip=0; strip<96; strip++) { + // + // run over all the strips of the plane + // + if ( strip >= 0 && strip < 32 ) encol[1][0] += enstrip[1][(2*i)][strip]; + if ( strip >= 32 && strip < 64 ) encol[1][1] += enstrip[1][(2*i)][strip]; + if ( strip >= 64 && strip < 96 ) encol[1][2] += enstrip[1][(2*i)][strip]; + entot[1] += enstrip[1][(2*i)][strip]; + // + }; + }; }; }; // if ( section.Contains("YE") ){ + // + // find the column hit in the first plane + // + Int_t ix = -1; + Int_t iy = -1; + if ( track_coordy[(2*0)+1][0] >= (-12.154+fM) && track_coordy[(2*0)+1][0] <= (-4.346-fM) ) iy = 0; + if ( track_coordy[(2*0)+1][0] >= ( -4.104+fM) && track_coordy[(2*0)+1][0] <= ( 3.704-fM) ) iy = 1; + if ( track_coordy[(2*0)+1][0] >= ( 3.946+fM) && track_coordy[(2*0)+1][0] <= (11.754-fM) ) iy = 2; + if ( cl2->cbar[(2*0)+1][0] >= (ye1 + fM1) && cl2->cbar[(2*0)+1][0] <= (ye2 - fM1) ) ix = 0; + if ( cl2->cbar[(2*0)+1][0] >= (ye3 + fM1) && cl2->cbar[(2*0)+1][0] <= (ye4 - fM1) ) ix = 1; + if ( cl2->cbar[(2*0)+1][0] >= (ye5 + fM1) && cl2->cbar[(2*0)+1][0] <= (ye6 - fM1) ) ix = 2; + if ( ix > -1 && iy > -1 ) fColYE = ix + iy*3; + // + // check event is inside XO acceptance + // for (Int_t i=0; i<11; i++) { - if - ( - (((track_coordy[2*i+1][0]>=(-12.154+fM))&& - (track_coordy[2*i+1][0]<=(-4.346-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye6 - fM1 ))))|| - - (((track_coordy[2*i+1][0]>=(-4.104+fM))&& - (track_coordy[2*i+1][0]<=(3.704-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye4 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye6 - fM1 ))))|| - - (((track_coordy[2*i+1][0]>=(3.946+fM))&& - (track_coordy[2*i+1][0]<=(11.754-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye2 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i+1][0]>=ye5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i+1][0]<=ye6 - fM1 )))) - ){ + if (( + ( track_coordy[(2*i)+1][0] >= (-12.154+fM) && track_coordy[(2*i)+1][0] <= (-4.346-fM) && (iy == 0 || multicol) ) || + ( track_coordy[(2*i)+1][0] >= ( -4.104+fM) && track_coordy[(2*i)+1][0] <= ( 3.704-fM) && (iy == 1 || multicol) ) || + ( track_coordy[(2*i)+1][0] >= ( 3.946+fM) && track_coordy[(2*i)+1][0] <= (11.754-fM) && (iy == 2 || multicol) ) + ) && ( + ( cl2->cbar[(2*i)+1][0] >= (ye1 + fM1) && cl2->cbar[(2*i)+1][0] <= (ye2 - fM1) && (ix == 0 || multicol) ) || + ( cl2->cbar[(2*i)+1][0] >= (ye3 + fM1) && cl2->cbar[(2*i)+1][0] <= (ye4 - fM1) && (ix == 1 || multicol) ) || + ( cl2->cbar[(2*i)+1][0] >= (ye5 + fM1) && cl2->cbar[(2*i)+1][0] <= (ye6 - fM1) && (ix == 2 || multicol) ) + )){ fYesel = true; fYeout = i; } else { @@ -554,51 +936,58 @@ break; }; }; - if ( !fYesel && fYeout >= fYemin && fYeout >= (Int_t)(fMax_planeye+fPl) ){ - if ( debug ) printf(" YE: this event is only partially contained: fYeout %i fYemin %i fMax_planeye + fPl %i \n",fYeout,fYemin,(Int_t)(fMax_planeye+fPl)); + // + if ( !fYesel && fYeout >= fYemin && fYeout >= (Int_t)(fMax_planeye+nPl) ){ + 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)); fPartsel = true; fYesel = true; }; if ( fYesel ){ - en_maxplane = 0.; - for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeye+fPl)) ;iplm++) en_maxplane += en_yep[iplm]; - fEnergyye = en_maxplane/fConv_rye; + for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeye+nPl)) ;iplm++) fYEen_maxplane += en_yep[iplm]; + fEnergyye = fYEen_maxplane/fConv_rye; + // + for (Int_t i=0;i<11;i++){ + for(strip=0; strip<96; strip++) { + // + // run over all the strips of the plane + // + if ( strip >= 0 && strip < 32 ) encol[0][0] += enstrip[0][(2*i)+1][strip]; + if ( strip >= 32 && strip < 64 ) encol[0][1] += enstrip[0][(2*i)+1][strip]; + if ( strip >= 64 && strip < 96 ) encol[0][2] += enstrip[0][(2*i)+1][strip]; + entot[0] += enstrip[0][(2*i)+1][strip]; + // + }; + }; + // }; }; // if ( section.Contains("YO") ){ + // + // find the column hit in the first plane + // + Int_t ix = -1; + Int_t iy = -1; + if ( track_coordy[(2*0)][0] >= (-11.954+fM) && track_coordy[(2*0)][0] <= (-4.146-fM) ) iy = 0; + if ( track_coordy[(2*0)][0] >= ( -3.904+fM) && track_coordy[(2*0)][0] <= ( 3.904-fM) ) iy = 1; + if ( track_coordy[(2*0)][0] >= ( 4.146+fM) && track_coordy[(2*0)][0] <= (11.954-fM) ) iy = 2; + if ( cl2->cbar[(2*0)][0] >= (yo1 + fM1) && cl2->cbar[(2*0)][0] <= (yo2 - fM1) ) ix = 0; + if ( cl2->cbar[(2*0)][0] >= (yo3 + fM1) && cl2->cbar[(2*0)][0] <= (yo4 - fM1) ) ix = 1; + if ( cl2->cbar[(2*0)][0] >= (yo5 + fM1) && cl2->cbar[(2*0)][0] <= (yo6 - fM1) ) ix = 2; + if ( ix > -1 && iy > -1 ) fColYO = ix + iy*3; + // + // check event is inside XO acceptance + // for (Int_t i=0; i<11; i++) { - 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 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 - ( - (((track_coordy[2*i][0]>=(-11.954+fM))&& - (track_coordy[2*i][0]<=(-4.146-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i][0]>=yo1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][0]>=yo3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][0]>=yo5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo6 - fM1 ))))|| - - (((track_coordy[2*i][0]>=(-3.904+fM))&& - (track_coordy[2*i][0]<=(+3.904-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i][0]>=yo1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo2 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][0]>=yo3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo4 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i][0]>=yo5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo6 - fM1 ))))|| - - (((track_coordy[2*i][0]>=(4.146+fM))&& - (track_coordy[2*i][0]<=(11.954-fM)))&& - (((L2->GetCaloLevel2()->cbar[2*i][0]>=yo1 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo2 - fM1))|| - ((L2->GetCaloLevel2()->cbar[2*i][0]>=yo3 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo4 - fM1 ))|| - ((L2->GetCaloLevel2()->cbar[2*i][0]>=yo5 + fM1)&& - (L2->GetCaloLevel2()->cbar[2*i][0]<=yo6 - fM1 )))) - ){ + if (( + ( track_coordy[(2*i)][0] >= (-11.954+fM) && track_coordy[(2*i)][0] <= (-4.146-fM) && (iy == 0 || multicol) ) || + ( track_coordy[(2*i)][0] >= ( -3.904+fM) && track_coordy[(2*i)][0] <= ( 3.904-fM) && (iy == 1 || multicol) ) || + ( track_coordy[(2*i)][0] >= ( 4.146+fM) && track_coordy[(2*i)][0] <= (11.954-fM) && (iy == 2 || multicol) ) + ) && ( + ( cl2->cbar[(2*i)][0] >= (yo1 + fM1) && cl2->cbar[(2*i)][0] <= (yo2 - fM1) && (ix == 0 || multicol) ) || + ( cl2->cbar[(2*i)][0] >= (yo3 + fM1) && cl2->cbar[(2*i)][0] <= (yo4 - fM1) && (ix == 1 || multicol) ) || + ( cl2->cbar[(2*i)][0] >= (yo5 + fM1) && cl2->cbar[(2*i)][0] <= (yo6 - fM1) && (ix == 2 || multicol) ) + )){ fYosel = true; fYoout = i; } else { @@ -606,24 +995,83 @@ break; }; }; - if ( !fYosel && fYoout >= fYomin && fYoout >= (Int_t)(fMax_planeyo+fPl) ){ - if ( debug ) printf(" YO: this event is only partially contained: fYoout %i fYomin %i fMax_planeyo + fPl %i \n",fYoout,fYomin,(Int_t)(fMax_planeyo+fPl)); + // + if ( !fYosel && fYoout >= fYomin && fYoout >= (Int_t)(fMax_planeyo+nPl) ){ + 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)); fPartsel = true; fYosel = true; }; if ( fYosel ){ - en_maxplane = 0.; - 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; + for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeyo+nPl)) ;iplm++) fYOen_maxplane += en_yop[iplm]; + fEnergyyo = fYOen_maxplane/fConv_ryo; + // + for (Int_t i=0;i<11;i++){ + for(strip=0; strip<96; strip++) { + // + // run over all the strips of the plane + // + if ( strip >= 0 && strip < 32 ) encol[0][0] += enstrip[0][(2*i)][strip]; + if ( strip >= 32 && strip < 64 ) encol[0][1] += enstrip[0][(2*i)][strip]; + if ( strip >= 64 && strip < 96 ) encol[0][2] += enstrip[0][(2*i)][strip]; + entot[0] += enstrip[0][(2*i)][strip]; + // + }; + }; }; }; // - fSel = ( fXesel || fYesel || fXosel || fYosel ); + // Count the number of sections in which the event is contained + // fCount = (Float_t)((Int_t)fXesel+(Int_t)fXosel+(Int_t)fYesel+(Int_t)fYosel); + // + if ( indep ){ + // + // independent mode, average the energy measurement and max plane of the contained sections + // + fSel = ( fXesel || fYesel || fXosel || fYosel ); + fMax_plane = (Float_t)(fMax_planeyo+fMax_planeye+fMax_planexo+fMax_planexe)/fCount; + fEnergy = (fEnergyxe+fEnergyyo+fEnergyye+fEnergyxo)/fCount; + // + } else { + // + // 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 + // + if ( fCount != fNumSec ){ + fSel = false; + } else { + fSel = true; + }; + fEnergy = (fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane)/fConv_rxo; + if ( fSel ){ + if ( fXesel ) fColumn = fColXE; + if ( fXosel ){ + if ( fColXO != fColumn && fColumn > -1 ){ + printf(" ERROR! mismatch in column number between different sections! fColumn %i fColXO %i \n",fColumn,fColXO); + } else { + fColumn = fColXO; + }; + }; + if ( fYesel ){ + if ( fColYE != fColumn && fColumn > -1 ){ + printf(" ERROR! mismatch in column number between different sections! fColumn %i fColYE %i \n",fColumn,fColYE); + } else { + fColumn = fColYE; + }; + }; + if ( fYosel ){ + if ( fColYO != fColumn && fColumn > -1 ){ + printf(" ERROR! mismatch in column number between different sections! fColumn %i fColYO %i \n",fColumn,fColYO); + } else { + fColumn = fColYO; + }; + }; + }; + }; + // + 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); if ( debug ) printf(" IsInside XE %i XO %i YE %i YO %i => SEL %i \n",fXesel,fXosel,fYesel,fYosel,fSel); - fMax_plane = (Float_t)(fMax_planeyo+fMax_planeye+fMax_planexo+fMax_planexe)/fCount; - fEnergy = (fEnergyxe+fEnergyyo+fEnergyye+fEnergyxo)/fCount; - if ( debug ) printf(" fMax_plane %f conv_r %f en_maxplane %f encalo %f \n",fMax_plane,fConv_rxo,en_maxplane,fEnergy); + // + // finish exit // return fSel; // @@ -637,6 +1085,8 @@ void CaloEnergy::Process(TString section){ // + // process the event + // TString ntr = section; if ( !L2 ){ printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); @@ -658,25 +1108,86 @@ newentry = true; }; // + // if we have already called this method for this event and no input changed then return fSel and exit + // if ( !newentry ) return; // + // process the event + // if ( debug ) printf(" Processing event at OBT %u PKT %u time %u section %s\n",OBT,PKT,atime,section.Data()); // - // this->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) // 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 ( fLong ){ - // use long fit to measure energy + if ( debug ) printf(" ==================================================================> LONGITUDINAL FIT! \n"); + // + // use long fit to measure energy + // + if ( this->IsInsideAcceptance(section) ){ + // + if ( debug ) printf(" ==================================================================> LONG INSIDE! \n"); + // + Float_t myene[2][22]; + memset(myene,0,(sizeof(Float_t))*2*22); + for (Int_t j=0; j<11; j++){ + if ( section.Contains("XE") ) myene[1][2*j] = en_xep[j]; + if ( section.Contains("YO") ) myene[0][2*j] = en_yop[j]; + if ( section.Contains("XO") ) myene[1][(2*j)+1] = en_xop[j]; + if ( section.Contains("YE") ) myene[0][(2*j)+1] = en_yep[j]; + }; + clong->UnMaskSections(); + if ( !(section.Contains("YE")) ) clong->MaskSection("YE"); + if ( !(section.Contains("YO")) ) clong->MaskSection("YO"); + if ( !(section.Contains("XO")) ) clong->MaskSection("XO"); + if ( !(section.Contains("XE")) ) clong->MaskSection("XE"); + clong->ForceNextFit(); + clong->SetEnergies(myene); + if ( debug ){ + clong->Fit(true); + } else { + clong->Fit(); + }; + if ( clong->GetLowerLimit() != 0. || clong->GetUpperLimit() != 0. ){ + fXOen_maxplane = clong->Get_defE0(); + } else { + fXOen_maxplane = clong->Get_E0(); + }; + fMax_plane = clong->Get_tmax(); + fYOen_maxplane = 0.; + fYEen_maxplane = 0.; + fXEen_maxplane = 0.; + fEnergy=fXOen_maxplane/fConv_rxo; + if ( fEnergy != fEnergy || clong->Get_fitresult() != 0 ) fEnergy = -1.; + // if ( fEnergy != fEnergy ) fEnergy = 1.; + // + } else { + // + // if the event is not in the acceptance, return a negative energy. + // + if ( debug ) printf(" Outside acceptance \n"); + fEnergy *= -1.; + // + }; + // } else { // + // use the energy measurement + // if ( this->IsInsideAcceptance(section) ){ // + // the event is good + // if ( debug ) printf(" XE %i XO %i YE %i YO %i \n",fXesel,fXosel,fYesel,fYosel); // } else { + // + // if the event is not in the acceptance, return a negative energy. + // if ( debug ) printf(" Outside acceptance \n"); - fEnergy = -1.; + fEnergy *= -1.; + // }; }; //