EmcMCfethe Class Reference

#include <EmcMCfethe.h>

List of all members.

Public Types
enum	{ m_oneProng = 1, m_twoProng = 2 }
enum	{ m_oneProng = 1, m_twoProng = 2 }
Public Member Functions
void	CollectBhabha ()
void	CollectBhabha ()
double	Ecorr (double En, double costhe)
double	Ecorr (double En, double costhe)
	EmcMCfethe (const std::string &name, ISvcLocator *pSvcLocator)
	EmcMCfethe (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	execute ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	finalize ()
double	findPhiDiff (double phi1, double phi2)
double	findPhiDiff (double phi1, double phi2)
int	index (int theta, int phi) const
int	index (int theta, int phi) const
void	initGeom ()
void	initGeom ()
StatusCode	initialize ()
StatusCode	initialize ()
bool	passed ()
bool	passed ()
void	RetreiveShowerList ()
void	RetreiveShowerList ()
void	SelectBhabha ()
void	SelectBhabha ()
void	setPassed (bool passed)
void	setPassed (bool passed)
	~EmcMCfethe ()
	~EmcMCfethe ()
Private Attributes
NTuple::Item< float >	En
NTuple::Item< float >	En
NTuple::Item< float >	EnA
NTuple::Item< float >	EnA
NTuple::Item< float >	EnB
NTuple::Item< float >	EnB
double	m_CalibConst [6240]
int	m_digiRangeCut
IEmcCalibConstSvc *	m_emcCalibConstSvc
IEmcCalibConstSvc *	m_emcCalibConstSvc
double	m_eNormCut
int	m_event
long int	m_events
std::string	m_fileDir
std::string	m_fileExt
IEmcRecGeoSvc *	m_iGeoSvc
IEmcRecGeoSvc *	m_iGeoSvc
int **	m_index
int **	m_index
double	m_inputConst [6240]
double	m_LATCut
int	m_maxNrXtalsShowerCut
double	m_minAngCut
int	m_minNrXtalsShowerCut
int	m_MsgFlag
int	m_nrShThreshCut
int	m_nXtals
int	m_oneProngsSelelected
bool	m_passed
double	m_phiDiffMaxCut
double	m_phiDiffMinCut
long int	m_rejected
int	m_run
bool	m_setAlwaysPassed
double	m_ShEneLeptonCut
double	m_ShEneThreshCut
list< EmcShower >	m_showerList
list< EmcShower >	m_showerList
long	m_showersAccepted
long int	m_taken
double	m_thetaBackwardCut
double	m_thetaDiffCut
double	m_thetaForwardMinAngCut
double	m_thetaIndexForwardCut
NTuple::Tuple *	m_tuple1
NTuple::Tuple *	m_tuple1
int	m_twoProngsSelected
bool	m_writeMVToFile
EmcBhabhaEvent *	myBhaEvt
EmcBhabhaEvent *	myBhaEvt
EmcBhaCalibData *	myCalibData
EmcBhaCalibData *	myCalibData
NTuple::Item< float >	Phi
NTuple::Item< float >	Phi
NTuple::Item< float >	PhiInd
NTuple::Item< float >	PhiInd
NTuple::Item< float >	Seed
NTuple::Item< float >	Seed
NTuple::Item< float >	The
NTuple::Item< float >	The
NTuple::Item< float >	TheInd
NTuple::Item< float >	TheInd

---

Member Enumeration Documentation

anonymous enum

Enumeration values: m_oneProng  m_twoProng  {m_oneProng=1, m_twoProng=2};

anonymous enum

Enumeration values: m_oneProng  m_twoProng  {m_oneProng=1, m_twoProng=2};

---

Constructor & Destructor Documentation

EmcMCfethe::EmcMCfethe (  const std::string &  name, ISvcLocator *  pSvcLocator )

:Algorithm(name, pSvcLocator), m_digiRangeCut(6), m_ShEneThreshCut(0.02), m_ShEneLeptonCut(1.), m_minNrXtalsShowerCut(10), m_maxNrXtalsShowerCut(70), m_phiDiffMinCut(0.05), m_phiDiffMaxCut(0.2), m_nrShThreshCut(20), m_thetaIndexForwardCut(0.), m_eNormCut(0.5), m_thetaDiffCut(0.04), m_LATCut(0.8), m_thetaForwardMinAngCut(0.35), m_minAngCut(0.3), m_thetaBackwardCut(2.3), m_showersAccepted(0), m_setAlwaysPassed(false), //-------------------- m_writeMVToFile(true), m_fileExt(""), m_fileDir("/ihepbatch/besdata/public/liucx/Calib/"), m_nXtals(6240), m_MsgFlag(0) { // Declare the properties declareProperty("digiRangeCut", m_digiRangeCut); declareProperty("ShEneThreshCut", m_ShEneThreshCut); declareProperty("ShEneLeptonCut", m_ShEneLeptonCut); declareProperty("minNrXtalsShowerCut", m_minNrXtalsShowerCut); declareProperty("maxNrXtalsShowerCut", m_maxNrXtalsShowerCut); declareProperty("phiDiffMinCut", m_phiDiffMinCut); declareProperty("phiDiffMaxCut", m_phiDiffMaxCut); declareProperty("nrShThreshCut", m_nrShThreshCut); declareProperty("thetaIndexForwardCut", m_thetaIndexForwardCut); declareProperty("eNormCut", m_eNormCut); declareProperty("thetaDiffCut", m_thetaDiffCut); declareProperty("LATCut", m_LATCut); declareProperty("thetaForwardMinAngCut", m_thetaForwardMinAngCut); declareProperty("minAngCut", m_minAngCut); declareProperty("thetaBackCut", m_thetaBackwardCut); declareProperty("setAlwaysPassed", m_setAlwaysPassed); //-------------------- declareProperty("writeMVToFile", m_writeMVToFile); declareProperty("fileExt", m_fileExt); declareProperty("fileDir", m_fileDir); declareProperty("MsgFlag", m_MsgFlag); int j = 0; m_index = new int*[56]; for (j=0;j<56;j++ ) { m_index[j] = new int[120]; for ( int k=0; k<120; k++) { m_index[j][k]=-1; } } m_iGeoSvc=0; for (int i=0; i<6240;i++) { m_inputConst[i] = 1.0; } }

EmcMCfethe::~EmcMCfethe (   )

{ if ( m_index != 0 ) { for (int i =0; i<56; i++ ) delete[] m_index[i]; delete[] m_index; m_index = 0; } }

EmcMCfethe::EmcMCfethe (  const std::string &  name, ISvcLocator *  pSvcLocator )

EmcMCfethe::~EmcMCfethe (   )

---

Member Function Documentation

void EmcMCfethe::CollectBhabha (   )

void EmcMCfethe::CollectBhabha (   )

{ MsgStream log(msgSvc(), name()); //check if the Bhabhas were found if ( 0 != myBhaEvt->positron()->found() || 0 != myBhaEvt->electron()->found() ) { m_taken++; //------------- electron found -------------------------------------- if (myBhaEvt->electron()->found() ) { En=myBhaEvt->electron()->shower().energy(); Seed=myBhaEvt->electron()->shower().maxima().energy(); TheInd= myBhaEvt->electron()->thetaIndex(); PhiInd= myBhaEvt->electron()->phiIndex(); The= myBhaEvt->electron()->theta(); Phi= myBhaEvt->electron()->phi(); EmcShower _theShower=EmcShower::EmcShower(); _theShower = myBhaEvt->electron()->shower(); list<EmcShDigi> _DigiList=_theShower.digiList(); list<EmcShDigi>::const_iterator _Digi; //------------------------------------------------------ //---- double loop over the digis to fill the matrix --- //first loop over Digis double enCalibBefore=0; double enCalibAfter=0; for (_Digi = _DigiList.begin(); _Digi != _DigiList.end(); _Digi++) { //get Xtal index //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55) //in Emc Bhabha Calibration unsigned int newThetaInd; if (_Digi->module()==0) newThetaInd = _Digi->thetaIndex(); if (_Digi->module()==1) newThetaInd = _Digi->thetaIndex() + 6; if (_Digi->module()==2) newThetaInd = 55 - _Digi->thetaIndex(); int DigiIndex = index(newThetaInd,_Digi->phiIndex()); enCalibBefore += (static_cast<double>(_Digi->energy())) * m_inputConst[DigiIndex]; enCalibAfter +=(static_cast<double>(_Digi->energy())) * m_inputConst[DigiIndex] *m_CalibConst[DigiIndex]; } EnB = enCalibBefore; EnA = enCalibAfter; m_tuple1->write(); } else log << MSG::WARNING << " Did not find MC deposited cluster energy " <<" for this cluster: thetaInd: " << myBhaEvt->electron()->thetaIndex() << " phiInd: " << myBhaEvt->electron()->phiIndex() << endl << "Will not use this cluster for the Emc " << "Bhabha calibration !" << endreq; //---------------- positron found ---------------------------------- if (myBhaEvt->positron()->found() ) { En=myBhaEvt->positron()->shower().energy(); Seed=myBhaEvt->positron()->shower().maxima().energy(); TheInd= myBhaEvt->positron()->thetaIndex(); PhiInd= myBhaEvt->positron()->phiIndex(); The= myBhaEvt->positron()->theta(); Phi= myBhaEvt->positron()->phi(); EmcShower _theShower=EmcShower::EmcShower(); _theShower = myBhaEvt->positron()->shower(); list<EmcShDigi> _DigiList=_theShower.digiList(); list<EmcShDigi>::const_iterator _Digi; //------------------------------------------------------ //---- double loop over the digis to fill the matrix --- //first loop over Digis double enCalibBefore=0; double enCalibAfter=0; for (_Digi = _DigiList.begin(); _Digi != _DigiList.end(); _Digi++) { //get Xtal index //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55) //in Emc Bhabha Calibration unsigned int newThetaInd; if (_Digi->module()==0) newThetaInd = _Digi->thetaIndex(); if (_Digi->module()==1) newThetaInd = _Digi->thetaIndex() + 6; if (_Digi->module()==2) newThetaInd = 55 - _Digi->thetaIndex(); int DigiIndex = index(newThetaInd,_Digi->phiIndex()); enCalibBefore += (static_cast<double>(_Digi->energy())) * m_inputConst[DigiIndex]; enCalibAfter +=(static_cast<double>(_Digi->energy())) * m_inputConst[DigiIndex] *m_CalibConst[DigiIndex]; } EnB = enCalibBefore; EnA = enCalibAfter; m_tuple1->write(); } else log << MSG::WARNING << " Did not find MC deposited cluster energy " << "for this cluster: thetaInd: " << myBhaEvt->positron()->thetaIndex() << " phiInd: " << myBhaEvt->positron()->phiIndex() << endl << "Will not use this cluster for the Emc " << "Bhabha calibration !" << endreq; } else { log << MSG::WARNING << " No Bhabha data for calibration found in event !" << endreq; m_rejected++; } }

double EmcMCfethe::Ecorr (  double  En, double  costhe )

double EmcMCfethe::Ecorr (  double  En, double  costhe )

{ /* double par[6]={0.0443492, -0.0104957, -0.0284476, -0.00817278, -6.48242e-06, 0.0480414}; */ double par[6]={ 0.047828, 0.0090778, -0.00254, -0.000847, -0.0001967, 0.044144 }; double ecor,lne; lne=log(En); ecor = En* exp(par[0] +par[1]*lne +par[2]*lne*lne +par[3]*lne*lne*lne +par[4]*costhe +par[5]*costhe*costhe); return ecor; }

StatusCode EmcMCfethe::execute (   )

StatusCode EmcMCfethe::execute (   )

{ MsgStream log(msgSvc(), name()); log << MSG::DEBUG << "in execute()" << endreq; //create the object that store the needed data of the Bhabha event int event, run; SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); if (!eventHeader) { log << MSG::FATAL << "Could not find Event Header" << endreq; return( StatusCode::FAILURE); } run=eventHeader->runNumber(); event=eventHeader->eventNumber(); cout<<"---------"<<event<<"---------"<<run<<endl; m_event = event; m_run = run; myBhaEvt = new EmcBhabhaEvent(); //retreive shower list of event RetreiveShowerList(); //select bhabha event SelectBhabha(); //collect bhabha event for digi-calibration of EMC //and fill matrix and vector of system of linear equations CollectBhabha(); delete myBhaEvt; myBhaEvt=0; return StatusCode::SUCCESS; }

StatusCode EmcMCfethe::finalize (   )

StatusCode EmcMCfethe::finalize (   )

{ MsgStream log(msgSvc(), name()); log << MSG::INFO << "in finalize()" << endreq; if ( m_writeMVToFile ) { delete myCalibData; myCalibData = 0; } cout<<"Event="<<m_events<<endl; cout<<"m_taken="<<m_taken<<endl; cout<<"m_rejected="<<m_rejected<<endl; return StatusCode::SUCCESS; }

double EmcMCfethe::findPhiDiff (  double  phi1, double  phi2 )

double EmcMCfethe::findPhiDiff (  double  phi1, double  phi2 )

{ double diff; diff = phi1 - phi2; //rad while( diff> pi ) diff -= twopi; while( diff< -pi ) diff += twopi; return diff; }

int EmcMCfethe::index (  int  theta, int  phi )  const [inline]

{ int val = ((m_index)[theta][phi]); return (val); }

int EmcMCfethe::index (  int  theta, int  phi )  const [inline]

{ int val = ((m_index)[theta][phi]); return (val); }

void EmcMCfethe::initGeom (   )

void EmcMCfethe::initGeom (   )

{ MsgStream log(msgSvc(), name()); int numberOfXtalsRing; EmcStructure* theEmcStruc=new EmcStructure() ; //number Of Theta Rings int nrOfTheRings = theEmcStruc->getNumberOfTheRings(); m_nXtals = theEmcStruc->getNumberOfXtals(); //init geometry for ( int the = theEmcStruc->startingTheta(); the< nrOfTheRings; the++) { numberOfXtalsRing = theEmcStruc->crystalsInRing((unsigned int) the ); for ( int phi=0; phi < numberOfXtalsRing; phi++) { m_index[the][phi] = theEmcStruc->getIndex( (unsigned int)the , (unsigned int)phi); } } log << MSG::INFO << " Emc geometry for Bhabha calibration initialized !! " << endl << "Number of Xtals: " << m_nXtals << endreq; delete theEmcStruc; //readout inputConst from file, only for MC data std::string constFile = m_fileDir; constFile += m_fileExt; constFile += "CalibInput.dat"; ifstream inputConstfile(constFile.c_str()); int ind=0; while( !inputConstfile.eof() ) { inputConstfile>>m_inputConst[ind]; //cout<<"ind="<<ind<<" "<<m_inputConst[ind]<<endl; ind++; } inputConstfile.close(); }

StatusCode EmcMCfethe::initialize (   )

StatusCode EmcMCfethe::initialize (   )

{ MsgStream log(msgSvc(), name()); log << MSG::INFO << "in initialize()" << endreq; //---------------------------------------<<<<<<<<<< m_taken=0; m_rejected=0; m_events=0; m_oneProngsSelelected=0; m_twoProngsSelected=0; //--------------------------------------->>>>>>>>>>> StatusCode status1; NTuplePtr nt1(ntupleSvc(),"FILE302/n1"); if ( nt1 ) m_tuple1 = nt1; else { m_tuple1=ntupleSvc()->book("FILE302/n1",CLID_ColumnWiseTuple,"elepos"); if( m_tuple1 ) { status1 = m_tuple1->addItem ("En",En); status1 = m_tuple1->addItem ("EnB",EnB); status1 = m_tuple1->addItem ("EnA",EnA); status1 = m_tuple1->addItem ("Seed",Seed); status1 = m_tuple1->addItem ("The",The); status1 = m_tuple1->addItem ("Phi",Phi); status1 = m_tuple1->addItem ("TheInd",TheInd); status1 = m_tuple1->addItem ("PhiInd",PhiInd); } else { // did not manage to book the N tuple.... log << MSG::ERROR <<"Cannot book N-tuple:" << long(m_tuple1) << endmsg; return StatusCode::FAILURE; } } //initialize Emc geometry to convert between index <=> theta,phi initGeom(); //create the object that holds the calibration data if ( m_writeMVToFile ) myCalibData = new EmcBhaCalibData(m_nXtals,m_MsgFlag); else myCalibData = 0; //get EmcRecGeoSvc ISvcLocator* svcLocator = Gaudi::svcLocator(); StatusCode sc = svcLocator->service("EmcRecGeoSvc",m_iGeoSvc); if(sc!=StatusCode::SUCCESS) { cout<<"Error: Can't get EmcRecGeoSvc"<<endl; } // use EmcCalibConstSvc StatusCode scCalib; scCalib = Gaudi::svcLocator() -> service("EmcCalibConstSvc", m_emcCalibConstSvc); if( scCalib != StatusCode::SUCCESS){ log << MSG::ERROR << "can not use EmcCalibConstSvc" << endreq; } else { std::cout << "Test EmcCalibConstSvc DigiCalibConst(0)= " << m_emcCalibConstSvc -> getDigiCalibConst(0) << std::endl; } int nConstEmc; nConstEmc= m_emcCalibConstSvc -> getDigiCalibConstNo() ; log << MSG::VERBOSE << " Have starting values for " << nConstEmc << " Xtals !" << endl << "Set the others to 1.0 ! " << endmsg; int numberXtal; //m_emcCalibConstSvc -> Dump(); for ( int i = 0; i< nConstEmc; i++ ) { numberXtal=i; m_CalibConst[numberXtal]= m_emcCalibConstSvc -> getDigiCalibConst(i); // cout<<"m_calibConst[numberXtal]"<<numberXtal //<<" "<<m_CalibConst[numberXtal]<<endl; } return StatusCode::SUCCESS; }

bool EmcMCfethe::passed (   )  [inline]

{ return m_passed;}

bool EmcMCfethe::passed (   )  [inline]

{ return m_passed;}

void EmcMCfethe::RetreiveShowerList (   )

void EmcMCfethe::RetreiveShowerList (   )

{ MsgStream log(msgSvc(), name()); // Retreive EmcRecShowerCol SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol"); if (!aShowerCol) { log << MSG::FATAL << "Could not find RecEmcShowerCol" << endreq; } // * * * * * * * * * * * * * * * * * * * * * * * * * ** * ** * //loop all showers of an event set EmcShower and EmcShDigi EmcShower* aShower =new EmcShower(); double ene,theta,phi,eseed; double showerX,showerY,showerZ; long int thetaIndex,phiIndex,numberOfDigis; RecEmcID showerId; unsigned int showerModule; HepPoint3D showerPosition(0,0,0); if ( ! m_showerList.empty()) m_showerList.clear(); RecEmcShowerCol::iterator iShowerCol; for(iShowerCol=aShowerCol->begin(); iShowerCol!=aShowerCol->end(); iShowerCol++){ //ene=(*iShowerCol)->energy(); //shower energy unit GeV ene=(*iShowerCol)->e5x5(); //shower5x5 energy unit GeV eseed=(*iShowerCol)->eSeed(); //unit GeV //cout<<"eseed="<<eseed<<endl; showerPosition = (*iShowerCol)->position(); theta = showerPosition.theta(); phi= showerPosition.phi(); showerX=showerPosition.x(); showerY=showerPosition.y(); showerZ=showerPosition.z(); //ene=Ecorr(ene, showerPosition.cosTheta())*ene; //shower energy correction showerId = (*iShowerCol)->getShowerId(); showerModule = EmcID::barrel_ec(showerId); //The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5) //module is defined by Endcap_east(0),Barrel(1),Endcap_west(2) thetaIndex=EmcID::theta_module(showerId); phiIndex=EmcID::phi_module(showerId); //------------------- EmcShDigi* aShDigi= new EmcShDigi(); EmcShDigi maxima =EmcShDigi::EmcShDigi(); list<EmcShDigi> shDigiList; RecEmcID cellId; unsigned int module; double digiEne, time, fraction, digiTheta, digiPhi; double digiX, digiY, digiZ; long int digiThetaIndex,digiPhiIndex; HepPoint3D digiPos(0,0,0); RecEmcFractionMap::const_iterator ciFractionMap; if ( ! shDigiList.empty()) shDigiList.clear(); double NewEn=0; for(ciFractionMap=(*iShowerCol)->getFractionMap5x5().begin(); ciFractionMap!=(*iShowerCol)->getFractionMap5x5().end(); ciFractionMap++){ digiEne = ciFractionMap->second.getEnergy(); //digit energy unit GeV //cout<<"digienergy="<<digiEne<<endl; time = ciFractionMap->second.getTime(); fraction = ciFractionMap->second.getFraction(); cellId=ciFractionMap->second.getCellId(); digiPos=m_iGeoSvc->GetCFrontCenter(cellId); digiTheta = digiPos.theta(); digiPhi = digiPos.phi(); digiX = digiPos.x(); digiY = digiPos.y(); digiZ = digiPos.z(); module=EmcID::barrel_ec(cellId); //The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5) //module is defined by Endcap_east(0),Barrel(1),Endcap_west(2) digiThetaIndex=EmcID::theta_module(cellId); digiPhiIndex = EmcID::phi_module(cellId); NewEn += digiEne*m_inputConst[m_index[digiThetaIndex][digiPhiIndex] ] *m_CalibConst[m_index[digiThetaIndex][digiPhiIndex] ]; //set EmcShDigi aShDigi->setEnergy(digiEne); aShDigi->setTheta(digiTheta); aShDigi->setPhi(digiPhi); aShDigi->setModule(module); aShDigi->setThetaIndex(digiThetaIndex); aShDigi->setPhiIndex(digiPhiIndex); aShDigi->setTime(time); aShDigi->setFraction(fraction); aShDigi->setWhere(digiPos); aShDigi->setY(digiX); aShDigi->setY(digiY); aShDigi->setZ(digiZ); shDigiList.push_back(*aShDigi); } shDigiList.sort(); //sort energy from small to large numberOfDigis = shDigiList.size(); maxima = *(--shDigiList.end()); //cout<<"maxima="<<maxima.energy()<<endl; //set EmcShower //ene=NewEn; aShower->setEnergy(ene); aShower->setTheta(theta); aShower->setPhi(phi); aShower->setModule(showerModule); aShower->setThetaIndex(thetaIndex); aShower->setPhiIndex(phiIndex); aShower->setNumberOfDigis(numberOfDigis); aShower->setDigiList(shDigiList); aShower->setMaxima(maxima); aShower->setWhere(showerPosition); aShower->setX(showerX); aShower->setY(showerY); aShower->setZ(showerZ); m_showerList.push_back(*aShower); /* list<EmcShDigi>::const_iterator thedigi; for (thedigi = shDigiList.begin(); thedigi != shDigiList.end(); thedigi++) { cout << "E:digi " << thedigi->energy() << endl; } */ // cluster energy threshold cut delete aShDigi; } m_showerList.sort(); //sort energy from small to large /* list<EmcShower>::const_iterator theshower; for (theshower = m_showerList.begin(); theshower != m_showerList.end(); theshower++) { cout<< "E:shower " << theshower->energy() << endl; } */ delete aShower; }

void EmcMCfethe::SelectBhabha (   )

void EmcMCfethe::SelectBhabha (   )

{ MsgStream log(msgSvc(), name()); //get the event-by-event channel status object //EmcStatus *EmcEventStatus =??? //if ( _debugging) // EmcEventStatus->print(); ???? setPassed(false); ++m_events; int nrClThresh = 0; EmcShower *cl1=new EmcShower(), *cl2=new EmcShower(); double newTheta1 = -1., newTheta2 = -1.; double thetaDiff_min=999.; double theta_cluster_min = 999., theta_cluster_max=0.; EmcShower *cl_tmp=new EmcShower(); bool bcl1=false, bcl2=false, bcl_tmp=false; double newTheta_tmp = -1., theta_cluster_min_tmp = 999.; int onlyOneFound = -1; //Shower list if (m_showerList.size() > 0) { //first take off all digis far away from the maxima digi of a shower //make a new list with cutt of showers list<EmcShower> goodShowersList; //----------------------------> //iterator for the bump list list<EmcShower>::const_iterator theShower; int nrGoodShowers = 0; //loop over bumps for (theShower = m_showerList.begin(); theShower != m_showerList.end(); theShower++) { log << MSG::INFO << "E: " << theShower->energy() << " " << "the: " << theShower->theta() << " " << "theI: " <<theShower->thetaIndex() << " " << "phi: " << theShower->phi() << " " << "phiI: " <<theShower->phiIndex() << " " << "nDig: " << theShower->numberOfDigis() << endreq; // cluster energy threshold cut if ( theShower->energy() > m_ShEneThreshCut ) { ++nrClThresh; goodShowersList.push_back( *theShower ); nrGoodShowers++; } } //------------------------------------------------------------------- //---- loop over the good bumps and find one/two for calibration ---- // if ( goodShowersList.size() > 0 ) { //iterator for the shower list list<EmcShower>::const_iterator theShower,theShower1; //loop over bumps for (theShower = goodShowersList.begin(); theShower != goodShowersList.end(); theShower++) { //high energy clusters if ( theShower->energy() >m_ShEneLeptonCut ) { //---- sort out the two high energy clusters that are most // back to back double newtheta = theShower->theta(); unsigned int newthetaInd = theShower->thetaIndex(); if ( newthetaInd < theta_cluster_min ) theta_cluster_min = newthetaInd; //keep the most backward bump in case we find only one or no //back-to-back bumps if ( newtheta > theta_cluster_max ) { *cl_tmp = *theShower; bcl_tmp=true; theta_cluster_max = newtheta; theta_cluster_min_tmp = newthetaInd; newTheta_tmp = newtheta; if ( onlyOneFound == -1 ) onlyOneFound = 1; } for (theShower1 = goodShowersList.begin(); theShower1 != goodShowersList.end(); theShower1++) { //not the same high energy bump if ( theShower != theShower1 && theShower1->energy() >m_ShEneLeptonCut ) { double newtheta1 = theShower1->theta(); //difference of the two clusters in phi double phiDiff; phiDiff = ( theShower->phi() - theShower1->phi() -pi ); if (phiDiff < -pi) phiDiff += twopi; if (phiDiff > pi) phiDiff -= twopi; //difference of the two clusters in theta double thetaDiff = fabs( theShower->theta() + theShower1->theta() -pi ); // if ( thetaDiff < thetaDiff_min // && // fabs( phiDiff ) > m_phiDiffMinCut // && // fabs( phiDiff ) < m_phiDiffMaxCut ) { thetaDiff_min = thetaDiff; onlyOneFound = 0; if ( newtheta > newtheta1 ) { *cl1 = *theShower; bcl1=true; newTheta1 = newtheta; *cl2 = *theShower1; bcl2=true; newTheta2 = newtheta1; } else { *cl1 = *theShower1; bcl1=true; newTheta1 = newtheta1; *cl2 = *theShower; bcl2=true; newTheta2 = newtheta; } // } } } } } //if we found only one bump or no back-to-back bumps if ( bcl_tmp!=0 && onlyOneFound == 1 ) { cl1 = cl_tmp; newTheta1 = newTheta_tmp; theta_cluster_min = theta_cluster_min_tmp; cout << "Only one cluster found !" << endl; } if ( bcl1==0 && bcl2==0 ) { cout << " No matched cluster found !" << endl; } //fill the Bhabha event if ( bcl1!=0 ) { //set properties myBhaEvt->setPositron()->setShower(*cl1); myBhaEvt->setPositron()->setTheta(newTheta1); myBhaEvt->setPositron()->setPhi(cl1->phi()); unsigned int newthetaInd ; //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55) //in Emc Bhabha Calibration if (cl1->module()==0) newthetaInd = cl1->thetaIndex(); if (cl1->module()==1) newthetaInd = cl1->thetaIndex() + 6; if (cl1->module()==2) newthetaInd = 55 - cl1->thetaIndex(); myBhaEvt->setPositron()->setThetaIndex(newthetaInd); unsigned int newphiInd=cl1->phiIndex(); myBhaEvt->setPositron()->setPhiIndex(newphiInd); myBhaEvt->setPositron()->setFound(true); m_showersAccepted++; log << MSG::INFO << name() << ": Positron: theta,phi energy " << myBhaEvt->positron()->theta() << "," << myBhaEvt->positron()->shower().phi() << " " << myBhaEvt->positron()->shower().energy() << endreq; } if ( bcl2!=0) { //set properties including vertex corrected theta myBhaEvt->setElectron()->setShower(*cl2); myBhaEvt->setElectron()->setTheta(newTheta2); myBhaEvt->setElectron()->setPhi(cl2->phi()); unsigned int newthetaInd; //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55) //in Emc Bhabha Calibration if (cl2->module()==0) newthetaInd = cl2->thetaIndex(); if (cl2->module()==1) newthetaInd = cl2->thetaIndex() + 6; if (cl2->module()==2) newthetaInd = 55 - cl2->thetaIndex(); myBhaEvt->setElectron()->setThetaIndex(newthetaInd); unsigned int newphiInd=cl2->phiIndex(); myBhaEvt->setElectron()->setPhiIndex(newphiInd); myBhaEvt->setElectron()->setFound(true); m_showersAccepted++; log << MSG::INFO << name() << ": Electron: theta,phi energy " << myBhaEvt->electron()->theta() << "," << myBhaEvt->electron()->shower().phi() << " " << myBhaEvt->electron()->shower().energy() << endreq; } } else //goodShowersList length > 0 log << MSG::INFO << " No good showers found !" << endreq; } //shower list length > 0 else log << MSG::WARNING << " No Showers found in event !" << endreq; if( bcl1!=0) delete cl1; if( bcl2!=0) delete cl2; if( bcl_tmp!=0) delete cl_tmp; }

void EmcMCfethe::setPassed (  bool  passed  )  [inline]

{ m_passed = passed;}

void EmcMCfethe::setPassed (  bool  passed  )  [inline]

{ m_passed = passed;}

---

Member Data Documentation

NTuple::Item<float> EmcMCfethe::En [private]

NTuple::Item<float> EmcMCfethe::En [private]

NTuple::Item<float> EmcMCfethe::EnA [private]

NTuple::Item<float> EmcMCfethe::EnA [private]

NTuple::Item<float> EmcMCfethe::EnB [private]

NTuple::Item<float> EmcMCfethe::EnB [private]

double EmcMCfethe::m_CalibConst [private]

int EmcMCfethe::m_digiRangeCut [private]

IEmcCalibConstSvc* EmcMCfethe::m_emcCalibConstSvc [private]

IEmcCalibConstSvc* EmcMCfethe::m_emcCalibConstSvc [private]

double EmcMCfethe::m_eNormCut [private]

int EmcMCfethe::m_event [private]

long int EmcMCfethe::m_events [private]

std::string EmcMCfethe::m_fileDir [private]

std::string EmcMCfethe::m_fileExt [private]

IEmcRecGeoSvc* EmcMCfethe::m_iGeoSvc [private]

IEmcRecGeoSvc* EmcMCfethe::m_iGeoSvc [private]

int** EmcMCfethe::m_index [private]

int** EmcMCfethe::m_index [private]

double EmcMCfethe::m_inputConst [private]

double EmcMCfethe::m_LATCut [private]

int EmcMCfethe::m_maxNrXtalsShowerCut [private]

double EmcMCfethe::m_minAngCut [private]

int EmcMCfethe::m_minNrXtalsShowerCut [private]

int EmcMCfethe::m_MsgFlag [private]

int EmcMCfethe::m_nrShThreshCut [private]

int EmcMCfethe::m_nXtals [private]

int EmcMCfethe::m_oneProngsSelelected [private]

bool EmcMCfethe::m_passed [private]

double EmcMCfethe::m_phiDiffMaxCut [private]

double EmcMCfethe::m_phiDiffMinCut [private]

long int EmcMCfethe::m_rejected [private]

int EmcMCfethe::m_run [private]

bool EmcMCfethe::m_setAlwaysPassed [private]

double EmcMCfethe::m_ShEneLeptonCut [private]

double EmcMCfethe::m_ShEneThreshCut [private]

list<EmcShower> EmcMCfethe::m_showerList [private]

list<EmcShower> EmcMCfethe::m_showerList [private]

long EmcMCfethe::m_showersAccepted [private]

long int EmcMCfethe::m_taken [private]

double EmcMCfethe::m_thetaBackwardCut [private]

double EmcMCfethe::m_thetaDiffCut [private]

double EmcMCfethe::m_thetaForwardMinAngCut [private]

double EmcMCfethe::m_thetaIndexForwardCut [private]

NTuple::Tuple* EmcMCfethe::m_tuple1 [private]

NTuple::Tuple* EmcMCfethe::m_tuple1 [private]

int EmcMCfethe::m_twoProngsSelected [private]

bool EmcMCfethe::m_writeMVToFile [private]

EmcBhabhaEvent* EmcMCfethe::myBhaEvt [private]

EmcBhabhaEvent* EmcMCfethe::myBhaEvt [private]

EmcBhaCalibData* EmcMCfethe::myCalibData [private]

EmcBhaCalibData* EmcMCfethe::myCalibData [private]

NTuple::Item<float> EmcMCfethe::Phi [private]

NTuple::Item<float> EmcMCfethe::Phi [private]

NTuple::Item<float> EmcMCfethe::PhiInd [private]

NTuple::Item<float> EmcMCfethe::PhiInd [private]

NTuple::Item<float> EmcMCfethe::Seed [private]

NTuple::Item<float> EmcMCfethe::Seed [private]

NTuple::Item<float> EmcMCfethe::The [private]

NTuple::Item<float> EmcMCfethe::The [private]

NTuple::Item<float> EmcMCfethe::TheInd [private]

NTuple::Item<float> EmcMCfethe::TheInd [private]

---

The documentation for this class was generated from the following files:

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Emc/EmcCalib/EmcBhaCalib/EmcBhaCalib-00-00-18/EmcBhaCalib/EmcMCfethe.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/InstallArea/include/EmcBhaCalib/EmcBhaCalib/EmcMCfethe.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Emc/EmcCalib/EmcBhaCalib/EmcBhaCalib-00-00-18/src/EmcMCfethe.cxx

---