#include <EmcMCBhaEvent.h>
|
|
00051 {m_oneProng=1, m_twoProng=2};
|
|
|
00051 {m_oneProng=1, m_twoProng=2};
|
|
||||||||||||
|
00072 :Algorithm(name, pSvcLocator), 00073 m_digiRangeCut(6), 00074 m_ShEneThreshCut(0.02), 00075 m_ShEneLeptonCut(1.), 00076 m_minNrXtalsShowerCut(10), 00077 m_maxNrXtalsShowerCut(70), 00078 m_phiDiffMinCut(0.05), 00079 m_phiDiffMaxCut(0.2), 00080 m_nrShThreshCut(20), 00081 m_eNormCut(0.5), 00082 m_thetaDiffCut(0.04), 00083 m_LATCut(0.8), 00084 m_minAngCut(0.3), 00085 m_showersAccepted(0), 00086 //-------------------- 00087 m_writeMVToFile(true), 00088 m_fileExt(""), 00089 m_fileDir("/ihepbatch/besdata/public/liucx/Calib/"), 00090 m_nXtals(6240), 00091 m_sigmaCut(1.), 00092 m_beamEnergy(1.548), 00093 m_MsgFlag(0) 00094 00095 { 00096 00097 // Declare the properties 00098 00099 declareProperty("digiRangeCut", m_digiRangeCut); 00100 00101 declareProperty("ShEneThreshCut", m_ShEneThreshCut); 00102 declareProperty("ShEneLeptonCut", m_ShEneLeptonCut); 00103 00104 declareProperty("minNrXtalsShowerCut", m_minNrXtalsShowerCut); 00105 declareProperty("maxNrXtalsShowerCut", m_maxNrXtalsShowerCut); 00106 declareProperty("phiDiffMinCut", m_phiDiffMinCut); 00107 declareProperty("phiDiffMaxCut", m_phiDiffMaxCut); 00108 declareProperty("nrShThreshCut", m_nrShThreshCut); 00109 declareProperty("eNormCut", m_eNormCut); 00110 declareProperty("thetaDiffCut", m_thetaDiffCut); 00111 declareProperty("LATCut", m_LATCut); 00112 declareProperty("minAngCut", m_minAngCut); 00113 00114 //-------------------- 00115 declareProperty("writeMVToFile", m_writeMVToFile); 00116 declareProperty("fileExt", m_fileExt); 00117 declareProperty("fileDir", m_fileDir); 00118 declareProperty("sigmaCut", m_sigmaCut); 00119 declareProperty("beamEnergy", m_beamEnergy); 00120 00121 declareProperty("MsgFlag", m_MsgFlag); 00122 00123 00124 int j = 0; 00125 m_index = new int*[56]; 00126 for (j=0;j<56;j++ ) { 00127 m_index[j] = new int[120]; 00128 for ( int k=0; k<120; k++) { 00129 m_index[j][k]=-1; 00130 } 00131 } 00132 00133 m_iGeoSvc=0; 00134 00135 for (int i=0; i<6240;i++) 00136 { 00137 m_inputConst[i] = 1.0; 00138 } 00139 00140 }
|
|
|
00145 {
00146
00147 if ( m_index != 0 ) {
00148 for (int i =0; i<56; i++ )
00149 delete[] m_index[i];
00150 delete[] m_index;
00151 m_index = 0;
00152 }
00153
00154
00155 }
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
00718 {
00719
00720 MsgStream log(msgSvc(), name());
00721
00722 //check if the Bhabhas were found
00723 if ( 0 != myBhaEvt->positron()->found() ||
00724 0 != myBhaEvt->electron()->found() ) {
00725
00726 m_taken++;
00727 //fill the EmcBhabhas
00728 double calibEnergy=0.;
00729 double energyError=0.;
00730
00731 //------------- electron found --------------------------------------
00732 if (myBhaEvt->electron()->found() ) {
00733
00734 calibEnergy = myBhaEvt->
00735 getDepoMCShowerEnergy(myBhaEvt->electron()->thetaIndex(),
00736 myBhaEvt->electron()->phiIndex(),
00737 m_corFun[myBhaEvt->electron()->thetaIndex()],
00738 m_beamEnergy);
00739
00740 if ( calibEnergy > 0 ) {
00741
00742 energyError = myBhaEvt->
00743 getErrorDepoMCShowerEnergy(myBhaEvt->electron()->thetaIndex(),
00744 myBhaEvt->electron()->phiIndex(),
00745 m_eSigma[myBhaEvt->electron()->thetaIndex()]);
00746
00747 } else
00748 log << MSG::WARNING << " Did not find MC deposited cluster energy "
00749 <<" for this cluster: thetaInd: "
00750 << myBhaEvt->electron()->thetaIndex()
00751 << " phiInd: "
00752 << myBhaEvt->electron()->phiIndex()
00753 << endl
00754 << "Will not use this cluster for the Emc "
00755 << "Bhabha calibration !"
00756 << endreq;
00757
00758 //set all that stuff in an EmcBhabha
00759 myBhaEvt->setElectron()->setErrorOnCalibEnergy(energyError);
00760 myBhaEvt->setElectron()->setCalibEnergy(calibEnergy);
00761
00762 //myBhaEvt->electron()->print();
00763
00764 } else
00765 log << MSG::INFO<< name()
00766 << ": Electron was not selected ! "
00767 << endreq;
00768
00769 calibEnergy=0.;
00770 energyError=0.;
00771
00772 //---------------- positron found ----------------------------------
00773 if (myBhaEvt->positron()->found() ) {
00774
00775 calibEnergy = myBhaEvt->
00776 getDepoMCShowerEnergy(myBhaEvt->positron()->thetaIndex(),
00777 myBhaEvt->positron()->phiIndex(),
00778 m_corFun[myBhaEvt->electron()->thetaIndex()],
00779 m_beamEnergy);
00780
00781 if ( calibEnergy > 0. ) {
00782
00783 energyError = myBhaEvt->
00784 getErrorDepoMCShowerEnergy(myBhaEvt->positron()->thetaIndex(),
00785 myBhaEvt->positron()->phiIndex(),
00786 m_eSigma[myBhaEvt->electron()->thetaIndex()]);
00787
00788 } else
00789 log << MSG::WARNING << " Did not find MC deposited cluster energy "
00790 << "for this cluster: thetaInd: "
00791 << myBhaEvt->positron()->thetaIndex()
00792 << " phiInd: "
00793 << myBhaEvt->positron()->phiIndex()
00794 << endl
00795 << "Will not use this cluster for the Emc "
00796 << "Bhabha calibration !"
00797 << endreq;
00798
00799
00800 //set all that stuff in an EmcBhabha
00801 myBhaEvt->setPositron()->setErrorOnCalibEnergy(energyError);
00802 myBhaEvt->setPositron()->setCalibEnergy(calibEnergy);
00803
00804 //myBhaEvt->positron()->print();
00805
00806 } else
00807 log << MSG::INFO << name()
00808 << ": Positron was not selected ! "
00809 << endreq;
00810
00811 //calculate elements of Matrix M and vector R from Bhabha event,
00812 //M is in SLAP triad format
00813 fillMatrix();
00814
00815 } else {
00816 log << MSG::WARNING << " No Bhabha data for calibration found in event !"
00817 << endreq;
00818 m_rejected++;
00819 }
00820
00821 }
|
|
|
|
|
|
00196 {
00197
00198 MsgStream log(msgSvc(), name());
00199 log << MSG::DEBUG << "in execute()" << endreq;
00200
00201 //create the object that store the needed data of the Bhabha event
00202
00203 int event, run;
00204
00205 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00206 if (!eventHeader) {
00207 log << MSG::FATAL << "Could not find Event Header" << endreq;
00208 return( StatusCode::FAILURE);
00209 }
00210 run=eventHeader->runNumber();
00211 event=eventHeader->eventNumber();
00212 cout<<"---------"<<event<<"---------"<<run<<endl;
00213 m_event = event;
00214 m_run = run;
00215
00216 myBhaEvt = new EmcBhabhaEvent();
00217
00218 //retreive shower list of event
00219 RetreiveShowerList();
00220 //select bhabha event
00221 SelectBhabha();
00222 //collect bhabha event for digi-calibration of EMC
00223 //and fill matrix and vector of system of linear equations
00224 CollectBhabha();
00225
00226 delete myBhaEvt;
00227 myBhaEvt=0;
00228
00229 return StatusCode::SUCCESS;
00230 }
|
|
|
|
|
|
00258 {
00259
00260 double E_tot = 2.0*m_beamEnergy;
00261
00262 double EP_tot_2= (E_tot*E_tot) /2.;
00263 double ene = EP_tot_2/E_tot;
00264 return ene;
00265 }
|
|
|
|
|
|
01094 {
01095
01096 ofstream calibOut;
01097 std::string calibFile = m_fileDir;
01098 calibFile += m_fileExt;
01099 calibFile += "CalibInput.dat";
01100 calibOut.open(calibFile.c_str());
01101 for(int i=0;i<6240;i++)
01102 calibOut<<RandGauss::shoot(1.0, 0.05)<< " ";
01103 }
|
|
|
|
|
|
00825 {
00826
00827 EmcBhabha _theBhabha=EmcBhabha::EmcBhabha();
00828 EmcShower _theShower=EmcShower::EmcShower();
00829 EmcShDigi _DigiMax=EmcShDigi::EmcShDigi();
00830 double _sigmaBha;
00831
00832 // ---- get the current calibration constants
00833
00834
00835 // ---- loop over the two particles
00836 for ( int i = 1; i <= 2; i++ ) {
00837
00838 if ( i == 2 ) _theBhabha = *(myBhaEvt->electron());
00839 else _theBhabha = *(myBhaEvt->positron());
00840 //-----------------------------------------------------------
00841 // ---- fill the matrix only if we found the particle and ---
00842 // ---- a energy to calibrate on -----
00843 double _sig=_theBhabha.errorOnCalibEnergy();
00844 double _calibEne=_theBhabha.calibEnergy();
00845 double _bhaEne=_theBhabha.shower().energy();
00846
00847 if ( _theBhabha.found()!=0 && _theBhabha.calibEnergy() > 0.
00848 && _bhaEne >= _calibEne - m_sigmaCut * _sig && _bhaEne <= _calibEne + m_sigmaCut * _sig) {
00849
00850 //the error (measurement + error on the energy to calibrate on)
00851 _sigmaBha = _theBhabha.sigma2();
00852 //cout<<"sigma2="<<_sigmaBha<<endl;
00853
00854 _theShower = _theBhabha.shower();
00855
00856 //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55)
00857 //in Emc Bhabha Calibration
00858
00859 _DigiMax = _theShower.maxima();
00860
00861 unsigned int newThetaIndMax;
00862 if (_DigiMax.module()==0) newThetaIndMax = _DigiMax.thetaIndex();
00863 if (_DigiMax.module()==1) newThetaIndMax = _DigiMax.thetaIndex() + 6;
00864 if (_DigiMax.module()==2) newThetaIndMax = 55 - _DigiMax.thetaIndex();
00865
00866 int maximaIndex=0;
00867 maximaIndex = index(newThetaIndMax,_DigiMax.phiIndex());
00868
00869 list<EmcShDigi> _DigiList=_theShower.digiList();
00870
00871 list<EmcShDigi>::const_iterator _Digi1,_Digi2;
00872
00873 //------------------------------------------------------
00874 //---- double loop over the digis to fill the matrix ---
00875
00876 //first loop over Digis
00877 for (_Digi1 = _DigiList.begin();
00878 _Digi1 != _DigiList.end(); _Digi1++) {
00879
00880 //get Xtal index
00881
00882 //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55)
00883 //in Emc Bhabha Calibration
00884 unsigned int newThetaInd1;
00885 if (_Digi1->module()==0) newThetaInd1 = _Digi1->thetaIndex();
00886 if (_Digi1->module()==1) newThetaInd1 = _Digi1->thetaIndex() + 6;
00887 if (_Digi1->module()==2) newThetaInd1 = 55 - _Digi1->thetaIndex();
00888
00889 int Digi1Index = index(newThetaInd1,_Digi1->phiIndex());
00890
00891 //calculate the vector with MC data
00892 double dvec = ( (static_cast<double>(_Digi1->energy())) * m_inputConst[Digi1Index] *
00893 _theBhabha.calibEnergy() / _sigmaBha);
00894
00895 //fill the vector
00896 if ( m_writeMVToFile )
00897 (myCalibData->vectorR(Digi1Index)) += dvec;
00898
00899 //count hits in Xtals and keep theta and phi
00900 if ( m_writeMVToFile )
00901 (myCalibData->xtalHits(Digi1Index))++;
00902
00903 //if ( *(_theShower.maxima()) == *(_Digi1) ) {
00904
00905 if ( Digi1Index == maximaIndex ) {
00906 if ( m_writeMVToFile ){
00907 (myCalibData->xtalHitsDir(Digi1Index))++;
00908 //cout<<"xtalHitsDir="<<myCalibData->xtalHitsDir(Digi1Index)<<endl;
00909 }
00910 }
00911
00912 //second loop over Digis
00913 for (_Digi2 = _Digi1;
00914 _Digi2 != _DigiList.end(); _Digi2++) {
00915
00916 unsigned int newThetaInd2;
00917 if (_Digi2->module()==0) newThetaInd2 = _Digi2->thetaIndex();
00918 if (_Digi2->module()==1) newThetaInd2 = _Digi2->thetaIndex() + 6;
00919 if (_Digi2->module()==2) newThetaInd2 = 55 - _Digi2->thetaIndex();
00920
00921 int Digi2Index = index(newThetaInd2, _Digi2->phiIndex());
00922
00923 //calculate the matrix element with MC data
00924 double val =
00925 static_cast<double>((( (_Digi1->energy() * m_inputConst[Digi1Index]) *
00926 (_Digi2->energy() * m_inputConst[Digi2Index])
00927 ))/_sigmaBha);
00928 if ( m_writeMVToFile )
00929 myCalibData->matrixMEle( Digi1Index, Digi2Index) += val;
00930
00931 }
00932 }
00933
00934 } //if paricle found and calibration energy
00935
00936 }//loop over particles
00937
00938 }
|
|
|
|
|
|
00233 {
00234
00235 MsgStream log(msgSvc(), name());
00236
00237 log << MSG::INFO << "in finalize()" << endreq;
00238
00239 //output the data of Matrix and vector to files
00240 OutputMV();
00241
00242 if ( m_writeMVToFile ) {
00243 delete myCalibData;
00244 myCalibData = 0;
00245 }
00246
00247 cout<<"Event="<<m_events<<endl;
00248 cout<<"m_taken="<<m_taken<<endl;
00249 cout<<"m_rejected="<<m_rejected<<endl;
00250 return StatusCode::SUCCESS;
00251 }
|
|
||||||||||||
|
|
|
||||||||||||
|
01083 {
01084 double diff;
01085 diff = phi1 - phi2; //rad
01086
01087 while( diff> pi ) diff -= twopi;
01088 while( diff< -pi ) diff += twopi;
01089
01090 return diff;
01091 }
|
|
||||||||||||
|
00071 {
00072 int val = ((m_index)[theta][phi]);
00073 return (val); }
|
|
||||||||||||
|
00071 {
00072 int val = ((m_index)[theta][phi]);
00073 return (val); }
|
|
|
|
|
|
00942 {
00943
00944 MsgStream log(msgSvc(), name());
00945
00946 int numberOfXtalsRing;
00947
00948 EmcStructure* theEmcStruc=new EmcStructure() ;
00949 theEmcStruc->setEmcStruc();
00950
00951 //number Of Theta Rings
00952 int nrOfTheRings = theEmcStruc->getNumberOfTheRings();
00953
00954 m_nXtals = theEmcStruc->getNumberOfXtals();
00955 //init geometry
00956 for ( int the = theEmcStruc->startingTheta(); the< nrOfTheRings; the++) {
00957
00958 numberOfXtalsRing = theEmcStruc->crystalsInRing((unsigned int) the );
00959
00960 for ( int phi=0; phi < numberOfXtalsRing; phi++) {
00961
00962 m_index[the][phi] = theEmcStruc->getIndex( (unsigned int)the , (unsigned int)phi);
00963
00964 }
00965
00966 }
00967
00968 log << MSG::INFO << " Emc geometry for Bhabha calibration initialized !! "
00969 << endl
00970 << "Number of Xtals: " << m_nXtals << endreq;
00971 delete theEmcStruc;
00972
00973 //readout inputConst from file, only for MC data
00974 std::string constFile = m_fileDir;
00975 constFile += m_fileExt;
00976 constFile += "CalibInput.dat";
00977 ifstream inputConstfile(constFile.c_str());
00978 int ind=0;
00979 while( !inputConstfile.eof() ) {
00980 inputConstfile>>m_inputConst[ind];
00981 //cout<<"ind="<<ind<<" "<<m_inputConst[ind]<<endl;
00982 ind++;
00983 }
00984 inputConstfile.close();
00985
00986 }
|
|
|
|
|
|
00157 {
00158
00159 MsgStream log(msgSvc(), name());
00160 log << MSG::INFO << "in initialize()" << endreq;
00161
00162 //---------------------------------------<<<<<<<<<<
00163 m_taken=0;
00164 m_rejected=0;
00165 m_events=0;
00166
00167 //--------------------------------------->>>>>>>>>>>
00168 //initialize Emc geometry to convert between index <=> theta,phi
00169 initGeom();
00170
00171 //create the object that holds the calibration data
00172 if ( m_writeMVToFile )
00173 myCalibData = new EmcBhaCalibData(m_nXtals,m_MsgFlag);
00174 else
00175 myCalibData = 0;
00176
00177 //get EmcRecGeoSvc
00178
00179 ISvcLocator* svcLocator = Gaudi::svcLocator();
00180 StatusCode sc = svcLocator->service("EmcRecGeoSvc",m_iGeoSvc);
00181 if(sc!=StatusCode::SUCCESS) {
00182 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00183 }
00184
00185 // read correction function from the file of 'cor.dat'
00186 readCorFun();
00187 // read Esigma(Itheta) from the file of 'sigma.dat'
00188 readEsigma();
00189 //generate fake calibration constants
00190 //fakeCalibConst();
00191
00192 return StatusCode::SUCCESS;
00193 }
|
|
|
|
|
|
00269 {
00270
00271 double value = 0.;
00272
00273 //find the largest digis
00274 double ene1=0, ene2=0;
00275
00276 list<EmcShDigi>::const_iterator digi, digi1=0,digi2 = 0;
00277
00278 //find the two digis with the largest energy
00279 for ( digi = theShower->digiList().begin();
00280 digi != theShower->digiList().end(); digi++) {
00281
00282 if (digi->energy() > ene1) {
00283
00284 digi2 = digi1;
00285 ene2 = ene1;
00286 digi1 = digi;
00287 ene1 = digi->energy();
00288
00289 } else
00290
00291 if (digi->energy() > ene2) {
00292 digi2 = digi;
00293 ene2 = digi->energy();
00294 }
00295 }
00296
00297 //now calculate LAT
00298 for ( digi = theShower->digiList().begin();
00299 digi != theShower->digiList().end(); digi++) {
00300
00301 if ( digi != digi1 && digi != digi2) {
00302
00303 Hep3Vector shower_pos(theShower->where().x(),
00304 theShower->where().y(),
00305 theShower->where().z());
00306 Hep3Vector digi_pos(digi->where().x(),
00307 digi->where().y(),
00308 digi->where().z());
00309 Hep3Vector r = digi_pos - shower_pos;
00310
00311 shower_pos *= 1.0/shower_pos.mag();
00312 r = r - ((r.dot(shower_pos)) *
00313 ( shower_pos ) );
00314 value+=(r.mag()*r.mag())*digi->energy();
00315 }
00316 }
00317
00318 value /= ((value +
00319 25.0*(digi1->energy() + digi2->energy()) ));
00320
00321 return value;
00322
00323 }
|
|
|
|
|
|
00990 {
00991
00992 MsgStream log(msgSvc(), name());
00993
00994 //check this first because I sometimes got two endJob transitions
00995 if ( myCalibData != 0 )
00996
00997 //if set write the matrix and vector to a file
00998 if ( m_writeMVToFile ) {
00999
01000 int count;
01001 char cnum[10];
01002 if (m_run<0){
01003 count = sprintf(cnum,"MC%d",-m_run);
01004 }
01005 if (m_run>=0){
01006 count = sprintf(cnum,"%d",m_run);
01007 }
01008
01009 std::string runnum="";
01010 runnum.append(cnum);
01011
01012 ofstream runnumberOut;
01013 std::string runnumberFile = m_fileDir;
01014 runnumberFile += m_fileExt;
01015 runnumberFile +="runnumbers.dat";
01016 runnumberOut.open(runnumberFile.c_str(),ios::out|ios::app);
01017
01018 ifstream runnumberIn;
01019 runnumberIn.open(runnumberFile.c_str());
01020 bool status=false;
01021 while( !runnumberIn.eof() ) {
01022
01023 std::string num;
01024 runnumberIn >> num;
01025 if (runnum==num) {
01026 status=true;
01027 log << MSG::INFO<< " the runnumber: "<<runnum
01028 <<" exists in the file runnumbers.dat" <<endreq;
01029 break;
01030 }else{
01031 status=false;
01032 log << MSG::INFO<< " the runnumber: "<<runnum
01033 <<" does not exist in the file runnumbers.dat" <<endreq;
01034 }
01035 }
01036 runnumberIn.close();
01037
01038 ofstream vectorOut;
01039 std::string vectorFile = m_fileDir;
01040 vectorFile += m_fileExt;
01041 vectorFile += runnum;
01042 vectorFile += "CalibVector.dat";
01043 vectorOut.open(vectorFile.c_str());
01044
01045 ofstream matrixOut;
01046 std::string matrixFile = m_fileDir;
01047 matrixFile += m_fileExt;
01048 matrixFile += runnum;
01049 matrixFile += "CalibMatrix.dat";
01050 matrixOut.open(matrixFile.c_str());
01051
01052 if ( vectorOut.good() && matrixOut.good() &&runnumberOut.good()) {
01053
01054 myCalibData->writeOut(matrixOut, vectorOut);
01055
01056 log << MSG::INFO<< " Wrote files "
01057 << (vectorFile) << " and "
01058 << (matrixFile) <<endreq;
01059
01060
01061 if ( !status ) {
01062 runnumberOut<<runnum<<"\n";
01063 log << MSG::INFO<< "Wrote files "<<runnumberFile<< endreq;
01064 }
01065
01066 } else
01067 log << MSG::WARNING << " Could not open vector and/or matrix file !"
01068 << endl
01069 << "matrix file : " << matrixFile << endl
01070 << "vector file : " << vectorFile
01071 << endreq;
01072
01073 vectorOut.close();
01074 matrixOut.close();
01075 runnumberOut.close();
01076 }
01077
01078 }
|
|
|
00066 { return m_passed;}
|
|
|
00066 { return m_passed;}
|
|
|
|
|
|
01106 {
01107
01108 ifstream corFunIn;
01109 std::string corFunFile = m_fileDir;
01110 corFunFile += m_fileExt;
01111 corFunFile += "cor.conf";
01112 corFunIn.open(corFunFile.c_str());
01113 for(int i=0;i<56;i++) {
01114 corFunIn>>m_corFun[i];
01115 cout<<"energy corFun="<<m_corFun[i]<<endl;
01116 }
01117 corFunIn.close();
01118 }
|
|
|
|
|
|
01121 {
01122
01123 ifstream EsigmaIn;
01124 std::string EsigmaFile = m_fileDir;
01125 EsigmaFile += m_fileExt;
01126 EsigmaFile += "sigma.conf";
01127 EsigmaIn.open(EsigmaFile.c_str());
01128 for(int i=0;i<56;i++) {
01129 EsigmaIn>>m_eSigma[i];
01130 cout<<"Sigma="<<m_eSigma[i]<<endl;
01131 }
01132 EsigmaIn.close();
01133 }
|
|
|
|
|
|
00327 {
00328
00329 MsgStream log(msgSvc(), name());
00330
00331 // Retreive RecEmcShowerCol
00332 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00333 if (!aShowerCol) {
00334 log << MSG::FATAL << "Could not find RecEmcShowerCol" << endreq;
00335
00336 }
00337
00338 // * * * * * * * * * * * * * * * * * * * * * * * * * ** * ** *
00339 //loop all showers of an event set EmcShower and EmcShDigi
00340
00341 EmcShower* aShower =new EmcShower();
00342
00343 double ene,theta,phi,eseed;
00344 double showerX,showerY,showerZ;
00345 long int thetaIndex,phiIndex,numberOfDigis;
00346
00347 RecEmcID showerId;
00348 unsigned int showerModule;
00349
00350 HepPoint3D showerPosition(0,0,0);
00351
00352 if ( ! m_showerList.empty()) m_showerList.clear();
00353
00354 RecEmcShowerCol::iterator iShowerCol;
00355 for(iShowerCol=aShowerCol->begin();
00356 iShowerCol!=aShowerCol->end();
00357 iShowerCol++){
00358
00359 //ene=(*iShowerCol)->energy(); //corrected shower energy unit GeV
00360 ene=(*iShowerCol)->e5x5(); //uncorrected 5x5 energy unit GeV
00361 eseed=(*iShowerCol)->eSeed(); //unit GeV
00362 //cout<<"eseed="<<eseed<<endl;
00363
00364 showerPosition = (*iShowerCol)->position();
00365 theta = showerPosition.theta();
00366 phi= showerPosition.phi();
00367 showerX=showerPosition.x();
00368 showerY=showerPosition.y();
00369 showerZ=showerPosition.z();
00370
00371 showerId = (*iShowerCol)->getShowerId();
00372 showerModule = EmcID::barrel_ec(showerId);
00373
00374 //The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5)
00375 //module is defined by Endcap_east(0),Barrel(1),Endcap_west(2)
00376
00377 thetaIndex=EmcID::theta_module(showerId);
00378
00379 phiIndex=EmcID::phi_module(showerId);
00380
00381 //-------------------
00382
00383 EmcShDigi* aShDigi= new EmcShDigi();
00384 EmcShDigi maxima =EmcShDigi::EmcShDigi();
00385 list<EmcShDigi> shDigiList;
00386 RecEmcID cellId;
00387 unsigned int module;
00388 double digiEne, time, fraction, digiTheta, digiPhi;
00389 double digiX, digiY, digiZ;
00390 long int digiThetaIndex,digiPhiIndex;
00391 HepPoint3D digiPos(0,0,0);
00392
00393 RecEmcFractionMap::const_iterator ciFractionMap;
00394
00395 if ( ! shDigiList.empty()) shDigiList.clear();
00396
00397 for(ciFractionMap=(*iShowerCol)->getFractionMap5x5().begin();
00398 ciFractionMap!=(*iShowerCol)->getFractionMap5x5().end();
00399 ciFractionMap++){
00400
00401 digiEne = ciFractionMap->second.getEnergy(); //digit energy unit GeV
00402
00403 time = ciFractionMap->second.getTime();
00404 fraction = ciFractionMap->second.getFraction();
00405
00406 cellId=ciFractionMap->second.getCellId();
00407
00408 digiPos=m_iGeoSvc->GetCFrontCenter(cellId);
00409 digiTheta = digiPos.theta();
00410 digiPhi = digiPos.phi();
00411 digiX = digiPos.x();
00412 digiY = digiPos.y();
00413 digiZ = digiPos.z();
00414
00415 module=EmcID::barrel_ec(cellId);
00416 //The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5)
00417 //module is defined by Endcap_east(0),Barrel(1),Endcap_west(2)
00418
00419 digiThetaIndex=EmcID::theta_module(cellId);
00420
00421 digiPhiIndex = EmcID::phi_module(cellId);
00422 /*
00423 cout<<"digiEne="<<digiEne<<"digiTheta= "<<digiTheta
00424 <<"digiPhi= "<<digiPhi<<"module= "<<module
00425 <<"digiThetaIndex= "<<digiThetaIndex
00426 <<"digiPhiIndex= "<<digiPhiIndex
00427 <<"time="<<time <<"fraction="<<fraction
00428 <<"digiPos="<<digiPos<<"digiX= "<<digiX
00429 <<"digiY= "<<digiY<<"digiZ "<<digiZ<<endl;
00430 */
00431 //set EmcShDigi
00432 aShDigi->setEnergy(digiEne);
00433 aShDigi->setTheta(digiTheta);
00434 aShDigi->setPhi(digiPhi);
00435 aShDigi->setModule(module);
00436 aShDigi->setThetaIndex(digiThetaIndex);
00437 aShDigi->setPhiIndex(digiPhiIndex);
00438 aShDigi->setTime(time);
00439 aShDigi->setFraction(fraction);
00440 aShDigi->setWhere(digiPos);
00441 aShDigi->setY(digiX);
00442 aShDigi->setY(digiY);
00443 aShDigi->setZ(digiZ);
00444
00445 shDigiList.push_back(*aShDigi);
00446
00447 }
00448 shDigiList.sort(); //sort energy from small to large
00449
00450 numberOfDigis = shDigiList.size();
00451
00452 maxima = *(--shDigiList.end());
00453 //cout<<"maxima="<<maxima.energy()<<endl;
00454 //set EmcShower
00455 aShower->setEnergy(ene);
00456 aShower->setTheta(theta);
00457 aShower->setPhi(phi);
00458 aShower->setModule(showerModule);
00459 aShower->setThetaIndex(thetaIndex);
00460 aShower->setPhiIndex(phiIndex);
00461 aShower->setNumberOfDigis(numberOfDigis);
00462 aShower->setDigiList(shDigiList);
00463 aShower->setMaxima(maxima);
00464 aShower->setWhere(showerPosition);
00465 aShower->setX(showerX);
00466 aShower->setY(showerY);
00467 aShower->setZ(showerZ);
00468 m_showerList.push_back(*aShower);
00469 delete aShDigi;
00470
00471 }
00472
00473 m_showerList.sort(); //sort energy from small to large
00474
00475 delete aShower;
00476
00477 }
|
|
|
|
|
|
00482 {
00483
00484 MsgStream log(msgSvc(), name());
00485 log << MSG::INFO << "Beam data: " << endl
00486 << " e- energy: " << m_beamEnergy
00487 << endl
00488 << " e+ energy: " << m_beamEnergy
00489 << endreq;
00490
00491
00492 setPassed(false);
00493
00494 ++m_events;
00495
00497 int nrClThresh = 0;
00498
00499 EmcShower cl1,cl2,cl_tmp;
00500 cl1=EmcShower::EmcShower();
00501 cl2=EmcShower::EmcShower();
00502 cl_tmp=EmcShower::EmcShower();
00503
00504 double newTheta1 = -1., newTheta2 = -1.;
00505 double thetaDiff_min=999.;
00506 double theta_cluster_min = 999., theta_cluster_max=0.;
00507
00508 bool bcl1=false, bcl2=false, bcl_tmp=false;
00509
00510 double newTheta_tmp = -1., theta_cluster_min_tmp = 999.;
00511 int onlyOneFound = -1;
00512
00513 //Shower list
00514
00515 if (m_showerList.size() > 0) {
00516
00517 //first take off all digis far away from the maxima digi of a shower
00518 //make a new list with cutt of showers
00519 list<EmcShower> goodShowersList;
00520
00521 //---------------------------->
00522 //iterator for the bump list
00523 list<EmcShower>::const_iterator theShower;
00524
00525 int nrGoodShowers = 0;
00526
00527 //loop over bumps
00528 for (theShower = m_showerList.begin();
00529 theShower != m_showerList.end(); theShower++) {
00530
00531 log << MSG::INFO << "E: " << theShower->energy() << " "
00532 << "the: " << theShower->theta() << " "
00533 << "theI: " <<theShower->thetaIndex() << " "
00534 << "phi: " << theShower->phi() << " "
00535 << "phiI: " <<theShower->phiIndex() << " "
00536 << "nDig: " << theShower->numberOfDigis()
00537 << endreq;
00538
00539 // cluster energy threshold cut
00540
00541 if ( theShower->energy() > m_ShEneThreshCut ) {
00542 ++nrClThresh;
00543 goodShowersList.push_back( *theShower );
00544 nrGoodShowers++;
00545 }
00546 }
00547
00548 //-------------------------------------------------------------------
00549 //---- loop over the good bumps and find one/two for calibration ----
00550 //
00551
00552 if ( goodShowersList.size() > 0 ) {
00553
00554 //iterator for the shower list
00555 list<EmcShower>::const_iterator theShower,theShower1;
00556
00557 //loop over bumps
00558 for (theShower = goodShowersList.begin();
00559 theShower != goodShowersList.end(); theShower++) {
00560
00561 //high energy clusters
00562 if ( theShower->energy() >m_ShEneLeptonCut ) {
00563
00564 //---- sort out the two high energy clusters that are most
00565 // back to back
00566
00567 double newtheta = theShower->theta();
00568 unsigned int newthetaInd = theShower->thetaIndex();
00569
00570 if ( newthetaInd < theta_cluster_min )
00571 theta_cluster_min = newthetaInd;
00572
00573 //keep the most backward bump in case we find only one or no
00574 //back-to-back bumps
00575 if ( newtheta > theta_cluster_max ) {
00576
00577 cl_tmp = *theShower;
00578 bcl_tmp=true;
00579 theta_cluster_max = newtheta;
00580 theta_cluster_min_tmp = newthetaInd;
00581 newTheta_tmp = newtheta;
00582 if ( onlyOneFound == -1 ) onlyOneFound = 1;
00583 }
00584
00585
00586 for (theShower1 = goodShowersList.begin();
00587 theShower1 != goodShowersList.end(); theShower1++) {
00588
00589 //not the same high energy bump
00590 if ( theShower != theShower1 &&
00591 theShower1->energy() >m_ShEneLeptonCut ) {
00592
00593 double newtheta1 = theShower1->theta();
00594
00595 //difference of the two clusters in phi
00596 double phiDiff;
00597 phiDiff = ( theShower->phi() - theShower1->phi() -pi );
00598 if (phiDiff < -pi) phiDiff += twopi;
00599 if (phiDiff > pi) phiDiff -= twopi;
00600
00601 //difference of the two clusters in theta
00602 double thetaDiff = fabs( theShower->theta() +
00603 theShower1->theta() -pi );
00604
00605 // if ( thetaDiff < thetaDiff_min
00606 // &&
00607 // fabs( phiDiff ) > m_phiDiffMinCut
00608 // &&
00609 // fabs( phiDiff ) < m_phiDiffMaxCut ) {
00610
00611 thetaDiff_min = thetaDiff;
00612 onlyOneFound = 0;
00613
00614 if ( newtheta > newtheta1 ) {
00615
00616 cl1 = *theShower;
00617 bcl1=true;
00618 newTheta1 = newtheta;
00619
00620 cl2 = *theShower1;
00621 bcl2=true;
00622 newTheta2 = newtheta1;
00623 } else {
00624
00625 cl1 = *theShower1;
00626 bcl1=true;
00627 newTheta1 = newtheta1;
00628
00629 cl2 = *theShower;
00630 bcl2=true;
00631 newTheta2 = newtheta;
00632
00633 }
00634
00635 // }
00636
00637 }
00638 }
00639 }
00640 }
00641
00642 //if we found only one bump or no back-to-back bumps
00643 if ( bcl_tmp!=0 && onlyOneFound == 1 ) {
00644 cl1 = cl_tmp;
00645 newTheta1 = newTheta_tmp;
00646 theta_cluster_min = theta_cluster_min_tmp;
00647 cout << "Only one cluster found !" << endl;
00648 }
00649
00650 if ( bcl1==0 && bcl2==0 ) {
00651 cout << " No matched cluster found !"
00652 << endl;
00653 }
00654
00655 //fill the Bhabha event
00656 if ( bcl1!=0 ) {
00657 //set properties
00658 myBhaEvt->setPositron()->setShower(cl1);
00659 myBhaEvt->setPositron()->setTheta(newTheta1);
00660 myBhaEvt->setPositron()->setPhi(cl1.phi());
00661 unsigned int newthetaInd ;
00662 //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55)
00663 //in Emc Bhabha Calibration
00664 if (cl1.module()==0) newthetaInd = cl1.thetaIndex();
00665 if (cl1.module()==1) newthetaInd = cl1.thetaIndex() + 6;
00666 if (cl1.module()==2) newthetaInd = 55 - cl1.thetaIndex();
00667
00668 myBhaEvt->setPositron()->setThetaIndex(newthetaInd);
00669
00670 unsigned int newphiInd=cl1.phiIndex();
00671 myBhaEvt->setPositron()->setPhiIndex(newphiInd);
00672 myBhaEvt->setPositron()->setFound(true);
00673 m_showersAccepted++;
00674
00675 log << MSG::INFO << name() << ": Positron: theta,phi energy "
00676 << myBhaEvt->positron()->theta() << ","
00677 << myBhaEvt->positron()->shower().phi() << " "
00678 << myBhaEvt->positron()->shower().energy()
00679 << endreq;
00680
00681 }
00682
00683 if ( bcl2!=0) {
00684 //set properties including vertex corrected theta
00685 myBhaEvt->setElectron()->setShower(cl2);
00686 myBhaEvt->setElectron()->setTheta(newTheta2);
00687 myBhaEvt->setElectron()->setPhi(cl2.phi());
00688 unsigned int newthetaInd;
00689 //The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55)
00690 //in Emc Bhabha Calibration
00691 if (cl2.module()==0) newthetaInd = cl2.thetaIndex();
00692 if (cl2.module()==1) newthetaInd = cl2.thetaIndex() + 6;
00693 if (cl2.module()==2) newthetaInd = 55 - cl2.thetaIndex();
00694
00695 myBhaEvt->setElectron()->setThetaIndex(newthetaInd);
00696 unsigned int newphiInd=cl2.phiIndex();
00697 myBhaEvt->setElectron()->setPhiIndex(newphiInd);
00698 myBhaEvt->setElectron()->setFound(true);
00699 m_showersAccepted++;
00700
00701 log << MSG::INFO << name() << ": Electron: theta,phi energy "
00702 << myBhaEvt->electron()->theta() << ","
00703 << myBhaEvt->electron()->shower().phi() << " "
00704 << myBhaEvt->electron()->shower().energy()
00705 << endreq;
00706 }
00707
00708 } else //goodShowersList length > 0
00709 log << MSG::INFO << " No good showers found !"
00710 << endreq;
00711 } //shower list length > 0
00712 else
00713 log << MSG::WARNING << " No Showers found in event !" << endreq;
00714
00715 }
|
|
|
00067 { m_passed = passed;}
|
|
|
00067 { m_passed = passed;}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1