#include <EmcBhaCalib.h>
The matrix and the vector can either be fetched from file(s). If they are read from files the file names are supposed to have the name CalibVector.dat and calibMatrix.dat with an extension at the beginning. Only these extensions (usually the runnumbers) need to be listed in the read-in-file that is specified with runNumberFile.
Definition at line 58 of file EmcBhaCalib.h.
| EmcBhaCalib::EmcBhaCalib | ( | const std::string & | name, | |
| ISvcLocator * | pSvcLocator | |||
| ) |
Definition at line 68 of file EmcBhaCalib.cxx.
References m_absoluteConstants, m_askForMatrixInversion, m_calibConst, m_calibConstUnred, m_convCrit, m_dirHitsCut, m_equationSolver, m_fileDir, m_fileExt, m_fitPolynom, m_MsgFlag, m_oldConstants, m_readDataFromDB, m_runNumberFile, m_tuple1, and m_writeToFile.
00069 :Algorithm(name, pSvcLocator), 00070 m_dirHitsCut(200), 00071 m_convCrit(0.000001), 00072 m_askForMatrixInversion(true), 00073 m_fitPolynom(true), //no used now 00074 m_writeToFile(false), 00075 m_readDataFromDB(false), 00076 m_equationSolver("GMR"), 00077 m_fileExt(""), 00078 m_fileDir("/ihepbatch/besdata/public/liucx/Calib/"), 00079 m_nrNonZeros(0), 00080 m_nrXtalsEnoughHits(0), 00081 m_runNumberFile("runnumbers.dat"), 00082 m_MsgFlag(0) 00083 00084 { 00085 00086 // Declare the properties 00087 declareProperty("equationSolver", m_equationSolver); 00088 declareProperty("dirHitsCut", m_dirHitsCut); 00089 declareProperty("writeToFile", m_writeToFile); 00090 declareProperty("fileExt", m_fileExt); 00091 declareProperty("fileDir", m_fileDir); 00092 declareProperty("readDataFromDB", m_readDataFromDB); 00093 declareProperty("askForMatrixInversion", m_askForMatrixInversion); 00094 declareProperty("fitPolynom", m_fitPolynom); 00095 declareProperty("convCrit", m_convCrit); 00096 declareProperty("runNumberFile", m_runNumberFile); 00097 declareProperty("MsgFlag", m_MsgFlag); 00098 00099 m_calibConst = new double[6240]; 00100 m_calibConstUnred = new double[6240]; 00101 m_absoluteConstants = new double[6240]; 00102 m_oldConstants = new double[6240]; 00103 00104 00105 //ntuple 00106 m_tuple1=0; 00107 00108 00109 00110 }
| EmcBhaCalib::~EmcBhaCalib | ( | ) | [virtual] |
Definition at line 115 of file EmcBhaCalib.cxx.
References m_absoluteConstants, m_calibConst, m_calibConstUnred, m_myCalibData, and m_oldConstants.
00115 { 00116 if ( 0 != m_calibConst) { 00117 delete [] m_calibConst; 00118 m_calibConst = 0; 00119 } 00120 if ( 0 != m_calibConstUnred) { 00121 delete [] m_calibConstUnred; 00122 m_calibConstUnred = 0; 00123 } 00124 if ( 0 != m_absoluteConstants) { 00125 delete [] m_absoluteConstants; 00126 m_absoluteConstants = 0; 00127 } 00128 if ( 0 != m_oldConstants) { 00129 delete [] m_oldConstants; 00130 m_oldConstants = 0; 00131 } 00132 if ( 0 != m_myCalibData) { 00133 delete m_myCalibData; 00134 m_myCalibData = 0; 00135 } 00136 00137 }
| void EmcBhaCalib::digiConstCor | ( | ) | [private] |
Definition at line 869 of file EmcBhaCalib.cxx.
References IEmcCalibConstSvc::getIndex(), genRecEmupikp::i, m_emcCalibConstSvc, m_myCalibData, m_oldConstants, EmcBhaCalibData::nXtals(), and deljobs::string.
00869 { 00870 00871 // ofstream Fuout; 00872 //std::string Fuoutfile; 00873 00874 //Fuoutfile="emccalibconst.txt"; 00875 //Fuout.open(Fuoutfile.c_str()); 00876 00877 double digiConst[6240]; 00878 00879 for(int ind=0; ind<m_myCalibData->nXtals(); ind++ ) { 00880 00881 digiConst[ind]=m_oldConstants[ind]; 00882 // Fuout<<m_oldConstants[ind]<<endl; 00883 } 00884 00885 // Fuout.close(); 00886 00887 ifstream ExpEneIn; 00888 std::string ExpEneFile; 00889 ExpEneFile = "cor.conf"; 00890 ExpEneIn.open(ExpEneFile.c_str()); 00891 00892 double ExpEne[56]; 00893 double Exp_barrel[44],Exp_east[6],Exp_west[6]; 00894 00895 for(int i=0;i<56;i++) { 00896 00897 ExpEneIn>>ExpEne[i]; 00898 ExpEne[i]=ExpEne[i]*1.843; 00899 if (i<6) Exp_east[i]=ExpEne[i]; 00900 if (i>=6&&i<50) Exp_barrel[i-6]=ExpEne[i]; 00901 if (i>=50) Exp_west[55-i]=ExpEne[i]; 00902 00903 //cout<<ExpEne[i]<<endl; 00904 } 00905 ExpEneIn.close(); 00906 00907 ifstream EDepEneIn; 00908 std::string EDepEneFile = "allxtal.dat"; 00909 EDepEneIn.open(EDepEneFile.c_str()); 00910 00911 double CorDigiConst[6240]; 00912 int ipart,ithe,iphi; 00913 double peak,sigma; 00914 double coeff; 00915 00916 for(int i=0;i<6240;i++) { 00917 00918 EDepEneIn>>ipart>>ithe>>iphi>>peak>>sigma; 00919 int ix=m_emcCalibConstSvc->getIndex(ipart,ithe,iphi); 00920 00921 if (ipart==0) coeff=peak/Exp_east[ithe]; 00922 if (ipart==1) coeff=peak/Exp_barrel[ithe]; 00923 if (ipart==2) coeff=peak/Exp_west[ithe]; 00924 cout<<coeff<<endl; 00925 //CorDigiConst[ix]=digiConst[ix]/coeff; 00926 CorDigiConst[ix]=coeff; 00927 } 00928 EDepEneIn.close(); 00929 00930 00931 //define tree fill the absolute digicalibconsts into the root-file 00932 00933 TTree* constr= new TTree("DigiCalibConst","DigiCalibConst"); 00934 constr->Branch("DigiCalibConst",CorDigiConst,"CorDigiConst[6240]/D"); 00935 constr->Fill(); 00936 00937 double aa[6240]; 00938 constr->SetBranchAddress("DigiCalibConst", aa); 00939 constr->GetEntry(0); 00940 00941 TFile fconst("EmcCalibConstCor.root", "recreate"); 00942 constr->Write(); 00943 00944 fconst.Close(); 00945 delete constr; 00946 00947 00948 }
| StatusCode EmcBhaCalib::execute | ( | ) |
Definition at line 190 of file EmcBhaCalib.cxx.
References Bes_Common::DEBUG, genRecEmupikp::i, Bes_Common::INFO, m_askForMatrixInversion, m_calibConst, m_calibConstUnred, m_calibrationSuccessful, m_dirHitsCut, m_myCalibData, m_nrXtalsEnoughHits, m_readDataFromDB, m_writeToFile, msgSvc(), ntupleOut(), EmcBhaCalibData::nXtals(), EmcBhaCalibData::nXtalsHit(), prepareConstants(), EmcBhaCalibData::printVec(), readInFromDB(), readInFromFile(), EmcBhaCalibData::reduce(), solveEqua(), Bes_Common::WARNING, writeOut(), writeOutConst(), and EmcBhaCalibData::xtalInd().
00190 { 00191 00192 MsgStream log(msgSvc(), name()); 00193 log << MSG::DEBUG << "in execute()" << endreq; 00194 //read in accumulated matrix/matrices 00195 if ( m_readDataFromDB ) { 00196 m_calibrationSuccessful = readInFromDB(); 00197 00198 log << MSG::INFO << "read in accumulated matrix from DB"<<endreq; 00199 00200 } else { 00201 m_calibrationSuccessful = readInFromFile(); 00202 00203 log << MSG::INFO <<"read in accumulated matrix from file"<<endreq; 00204 00205 } 00206 00207 if ( m_calibrationSuccessful == true ) { 00208 00209 //ask first if to do a matrix inversion 00210 if (m_askForMatrixInversion==true) { 00211 00212 m_calibrationSuccessful = false; 00213 log << MSG::INFO << " Well, we have " 00214 << m_nrXtalsEnoughHits << " crystals whith more " 00215 << "than " << m_dirHitsCut 00216 << " direct hits. Do you want do " 00217 << "a matrix inversion ? [N]" << endreq; 00218 00219 00220 m_calibrationSuccessful = true; 00221 00222 } 00223 00224 if ( m_calibrationSuccessful == true ) { 00225 00226 //write added matrix and vector to files 00227 if ( m_writeToFile == true) { 00228 writeOut(); 00229 } 00230 //cout <<"writeout"<<endl; 00231 m_myCalibData->printVec(2); 00232 00233 //reduce to continious array of only non zeros elements for SLAP 00234 if ( m_calibrationSuccessful = m_myCalibData->reduce() ) { 00235 00236 cout<<"m_calibrationSuccessful="<<m_calibrationSuccessful<<endl; 00237 00238 if ( m_myCalibData->nXtalsHit() != m_myCalibData->nXtals() ){ 00239 log << MSG::INFO << " Reduced number of Xtals for calibration: " 00240 << m_myCalibData->nXtalsHit() 00241 << endreq; 00242 } 00243 cout<<"m_myCalibData->nXtalsHit()="<<m_myCalibData->nXtalsHit() 00244 <<"m_myCalibData->nXtals()="<<m_myCalibData->nXtals()<<endl; 00245 m_myCalibData->printVec(10); 00246 00247 //solve matrix equation 00248 if ( m_calibrationSuccessful = solveEqua() ) { 00249 00250 for (int i=0; i<m_myCalibData->nXtalsHit(); i++){ 00251 //fill an array of constants for storage in database 00252 m_calibConstUnred[m_myCalibData->xtalInd(i)] 00253 = m_calibConst[i]; 00254 00255 // if (m_myCalibData->xtalHitsDir(i)>10) 00256 // cout<<"Index,calibconst="<<" "<<i <<" "<<m_myCalibData->xtalInd(i) 00257 // <<" "<<m_calibConstUnred[m_myCalibData->xtalInd(i)] 00258 // <<" "<<m_calibConstUnred[m_myCalibData->xtalInd(i)]* 00259 // m_oldConstants[m_myCalibData->xtalInd(i)]<<endl; 00260 00261 } 00262 //do validation, fit polynom if necessary, fill CalList 00263 prepareConstants(); 00264 00265 //if wanted write constants to plain ASCII file 00266 if ( m_writeToFile==true){ 00267 writeOutConst(); 00268 } 00269 00270 ntupleOut(); 00271 00272 } 00273 00274 } else { 00275 log << MSG::WARNING << " Reduction of the Emc Bhabha calibration matrix FAILED !" 00276 << endl 00277 << " Will NOT perform Emc Bhabha calibration !" 00278 << endreq; 00279 return( StatusCode::FAILURE); 00280 } 00281 } 00282 } else { 00283 log << MSG::WARNING << " ERROR in reading in Emc Bhabha calibration data !" 00284 << endl 00285 << " Will NOT perform Emc Bhabha calibration !" 00286 << endreq; 00287 return( StatusCode::FAILURE); 00288 } 00289 00290 return StatusCode::SUCCESS; 00291 }
| StatusCode EmcBhaCalib::finalize | ( | ) |
Definition at line 294 of file EmcBhaCalib.cxx.
References Bes_Common::INFO, and msgSvc().
00294 { 00295 00296 MsgStream log(msgSvc(), name()); 00297 00298 00299 log << MSG::INFO << "in endRun()" << endreq; 00300 00301 00302 return StatusCode::SUCCESS; 00303 }
| void EmcBhaCalib::help | ( | ) | [virtual] |
Definition at line 307 of file EmcBhaCalib.cxx.
References Bes_Common::INFO, and msgSvc().
00307 { 00308 00309 MsgStream log(msgSvc(), name()); 00310 00311 log << MSG::INFO<< "Performs the Chi square fit of the accumulated " 00312 << endreq; 00313 }
| void EmcBhaCalib::initCalibConst | ( | ) | [private] |
Definition at line 317 of file EmcBhaCalib.cxx.
References genRecEmupikp::i, m_absoluteConstants, m_calibConst, m_calibConstUnred, m_emcCalibConstSvc, m_myCalibData, m_oldConstants, msgSvc(), EmcBhaCalibData::nXtals(), and VERBOSE.
Referenced by initialize().
00317 { 00318 00319 00320 MsgStream log(msgSvc(), name()); 00321 00322 for ( int i = 0; i< m_myCalibData->nXtals(); i++ ) { 00323 m_calibConst[i] = 1.; 00324 } 00325 00326 ifstream inConst("EmcCalib.Const"); 00327 if ( !inConst ) { 00328 log << MSG::VERBOSE 00329 << " No starting values for calibration constants found ! " 00330 << "(EmcCalib.Const)" 00331 << endl 00332 << "Set them to 1.0 ! " << endreq; 00333 } else { 00334 log << MSG::VERBOSE << " Read in starting values of calibration constants !" 00335 << endreq; 00336 00337 int nConst; 00338 inConst >> nConst; 00339 log << MSG::VERBOSE << " Have starting values for " 00340 << nConst << " Xtals !" << endl 00341 << "Set the others to 1.0 ! " << endmsg; 00342 00343 int numberXtal; 00344 for ( int i = 0; i< nConst; i++ ) { 00345 inConst >> numberXtal >> m_calibConst[numberXtal]; 00346 } 00347 } 00348 00349 int nConstEmc; 00350 00351 nConstEmc= m_emcCalibConstSvc -> getDigiCalibConstNo() ; 00352 00353 if ( nConstEmc!=6240) cout<<"number of calibconst="<< nConstEmc<<endl; 00354 00355 // log << MSG::VERBOSE << " Have starting values for " 00356 // << nConstEmc << " Xtals !" << endl 00357 // << "Set the others to 1.0 ! " << endmsg; 00358 00359 for ( int i = 0; i< nConstEmc; i++ ) { 00360 //cout<<i<<" " 00361 //<<m_emcCalibConstSvc->getThetaIndex(i)<<" " 00362 //<<m_emcCalibConstSvc->getPhiIndex(i) 00363 //<<" "<<m_emcCalibConstSvc->getEmcID(i)<<endl; 00364 m_calibConst[i] = m_emcCalibConstSvc -> getDigiCalibConst(i); 00365 //m_calibConst[i] =5.0; 00366 m_oldConstants[i]=m_emcCalibConstSvc -> getDigiCalibConst(i); 00367 00368 // initialize m_absoluteConstants as m_oldConstants. 00369 //m_absoluteConstants[i] =m_emcCalibConstSvc -> getDigiCalibConst(i); 00370 m_absoluteConstants[i] =1.0; 00371 // initialize m_calibConstUnred as 1. 00372 m_calibConstUnred[i] =1.0; 00373 } 00374 00375 00376 }
| StatusCode EmcBhaCalib::initialize | ( | ) |
Definition at line 140 of file EmcBhaCalib.cxx.
References calibConst, err, calibUtil::ERROR, gi0, Bes_Common::INFO, initCalibConst(), ixtal, m_calibrationSuccessful, m_emcCalibConstSvc, m_MsgFlag, m_myCalibData, m_tuple1, msgSvc(), nhitxtal, and ntupleSvc().
00140 { 00141 00142 MsgStream log(msgSvc(), name()); 00143 log << MSG::INFO << "in initialize()" << endreq; 00144 00145 m_myCalibData = new EmcBhaCalibData(6240, m_MsgFlag); 00146 00147 m_calibrationSuccessful = false; 00148 00149 StatusCode status1; 00150 00151 NTuplePtr nt1(ntupleSvc(),"FILE302/n1"); 00152 if ( nt1 ) m_tuple1 = nt1; 00153 else { 00154 m_tuple1=ntupleSvc()->book("FILE302/n1",CLID_ColumnWiseTuple,"Calib"); 00155 if( m_tuple1 ) { 00156 00157 status1 = m_tuple1->addItem ("ixtal",ixtal); 00158 status1 = m_tuple1->addItem ("gi0",gi0); 00159 status1 = m_tuple1->addItem ("calibConst",calibConst); 00160 status1 = m_tuple1->addItem ("err",err); 00161 status1 = m_tuple1->addItem ("nhitxtal",nhitxtal); 00162 00163 } 00164 else { // did not manage to book the N tuple.... 00165 log << MSG::ERROR <<"Cannot book N-tuple:" << long(m_tuple1) << endmsg; 00166 return StatusCode::FAILURE; 00167 } 00168 } 00169 00170 // use EmcCalibConstSvc 00171 StatusCode scCalib; 00172 scCalib = Gaudi::svcLocator() -> service("EmcCalibConstSvc", m_emcCalibConstSvc); 00173 if( scCalib != StatusCode::SUCCESS){ 00174 log << MSG::ERROR << "can not use EmcCalibConstSvc" << endreq; 00175 } 00176 else { 00177 std::cout << "Test EmcCalibConstSvc DigiCalibConst(0)= " 00178 << m_emcCalibConstSvc -> getDigiCalibConst(0) << std::endl; 00179 } 00180 00181 //init starting values for calibration constants from file or set to 1 00182 initCalibConst(); 00183 00184 //digiConstCor(); 00185 00186 return StatusCode::SUCCESS; 00187 }
| void EmcBhaCalib::ntupleOut | ( | ) | [private] |
Definition at line 828 of file EmcBhaCalib.cxx.
References calibConst, err, gi0, genRecEmupikp::i, ixtal, m_calibConstUnred, m_myCalibData, m_oldConstants, m_tuple1, nhitxtal, EmcBhaCalibData::nXtalsHit(), EmcBhaCalibData::xtalHitsDir(), and EmcBhaCalibData::xtalInd().
Referenced by execute().
00828 { 00829 00830 00831 for (int i=0; i < m_myCalibData->nXtalsHit(); i++) { 00832 int xtalIndex=m_myCalibData->xtalInd(i); 00833 00834 gi0 = m_oldConstants[xtalIndex]; 00835 00836 if (gi0<0) cout<<"gi0="<<gi0<<endl; 00837 err = (m_calibConstUnred[xtalIndex]-gi0)/gi0; 00838 calibConst=m_calibConstUnred[xtalIndex]; 00839 ixtal=xtalIndex; 00840 nhitxtal=m_myCalibData->xtalHitsDir(i); 00841 m_tuple1->write(); 00842 } 00843 00844 }
| bool EmcBhaCalib::prepareConstants | ( | ) | [private] |
Definition at line 741 of file EmcBhaCalib.cxx.
References IEmcCalibConstSvc::getIndex(), genRecEmupikp::i, m_absoluteConstants, m_calibConstUnred, m_dirHitsCut, m_emcCalibConstSvc, m_myCalibData, m_oldConstants, EmcBhaCalibData::nXtals(), EmcBhaCalibData::nXtalsHit(), EmcBhaCalibData::xtalHitsDir(), and EmcBhaCalibData::xtalInd().
Referenced by execute().
00741 { 00742 00743 bool successfull=false; 00744 00745 //the old constant 00746 //float theCalConst = 1.0; 00747 int chanIndex; 00748 00749 for ( int i = 0; i< m_myCalibData->nXtalsHit(); i++ ) { 00750 00751 chanIndex=m_myCalibData->xtalInd(i); //xtalInd(i) array of xtal indices 00752 00754 // m_calibConstUnred[chanIndex]=1; 00755 //if (chanIndex==689) m_calibConstUnred[chanIndex]=1.47; 00756 //if (chanIndex==689) 00757 // m_oldConstants[chanIndex]=m_oldConstants[chanIndex]*1.47; 00758 //---- get the new absolute constant ---- 00759 //set it to 0 if we have not enough hits -> we'll keep the old constant 00760 if ( m_myCalibData->xtalHitsDir(i) < m_dirHitsCut ) 00761 m_absoluteConstants[chanIndex] = m_oldConstants[chanIndex]; 00762 else 00763 m_absoluteConstants[chanIndex] = 00764 m_oldConstants[chanIndex] * m_calibConstUnred[chanIndex]; 00765 } 00766 00767 00768 double DigiConst[6240]; 00769 int IxtalNumber[6240]; 00770 00771 for(int ind=0; ind<m_myCalibData->nXtals(); ind++ ) { 00772 00773 DigiConst[ind]=m_absoluteConstants[ind]; 00774 // cout<<"ind="<<ind<<"\t"<< DigiConst[ind]<<endl; 00775 } 00776 00777 int ixtal1,ixtal2,ixtal3; 00778 ixtal1=m_emcCalibConstSvc->getIndex(1,20,26); 00779 ixtal2=m_emcCalibConstSvc->getIndex(1,23,26); 00780 ixtal3=m_emcCalibConstSvc->getIndex(1,24,26); 00781 //cout<<ixtal1<<" "<<ixtal2<<" "<<ixtal3<<" "<<endl; 00782 for(int ind=0; ind<m_myCalibData->nXtals(); ind++ ) { 00783 00784 IxtalNumber[ind]=-1; 00785 /* 00786 if (ind==ixtal1) IxtalNumber[ind]=ixtal3; 00787 if (ind==ixtal2) IxtalNumber[ind]=ixtal1; 00788 if (ind==ixtal3) IxtalNumber[ind]=ixtal2; 00789 */ 00790 } 00791 00792 TFile fconst("EmcCalibConst.root", "recreate"); 00793 00794 //define tree fill the absolute digicalibconsts into the root-file 00795 TTree* constr= new TTree("DigiCalibConst","DigiCalibConst"); 00796 constr->Branch("DigiCalibConst",DigiConst,"DigiConst[6240]/D"); 00797 constr->Branch("IxtalNumber",IxtalNumber,"IxtalNumber[6240]/I"); 00798 00799 constr->Fill(); 00800 /* 00801 00802 double aa[6240]; 00803 int Iaa[6240]; 00804 constr->SetBranchAddress("DigiCalibConst", aa); 00805 constr->SetBranchAddress("IxtalNumber", Iaa); 00806 constr->GetEntry(0); 00807 cout<<Iaa[10]<<endl; 00808 cout<<aa[10]<<endl; 00809 00810 cout<<constr->GetEntry(0)<<endl; 00811 00812 */ 00813 00814 constr->Write(); 00815 00816 delete constr; 00817 00818 fconst.Close(); 00819 00820 // end tree 00821 00822 successfull=true; 00823 return successfull; 00824 00825 }
| void EmcBhaCalib::printInfo | ( | std::string | xtalHitsDirFile | ) | [private] |
Definition at line 848 of file EmcBhaCalib.cxx.
References genRecEmupikp::i, m_dirHitsCut, m_fileDir, m_fileExt, m_myCalibData, EmcBhaCalibData::nXtalsHit(), deljobs::string, EmcBhaCalibData::xtalHitsDir(), and EmcBhaCalibData::xtalInd().
00848 { 00849 00850 ofstream xtalHitsDirOut; 00851 std::string xtalHitsDirFile = m_fileDir; 00852 xtalHitsDirFile += m_fileExt; 00853 xtalHitsDirFile += fileName; 00854 xtalHitsDirOut.open(xtalHitsDirFile.c_str()); 00855 00856 //nXtalsHit() is number of xtals hit 00857 for (int i=0; i < m_myCalibData->nXtalsHit(); i++) { 00858 if ( m_myCalibData->xtalHitsDir(i) > m_dirHitsCut ) 00859 { 00860 int chanindex=m_myCalibData->xtalInd(i); 00861 xtalHitsDirOut<<chanindex<<endl; 00862 } 00863 } 00864 xtalHitsDirOut.close(); 00865 00866 }
| bool EmcBhaCalib::readInFromDB | ( | ) | [private] |
| bool EmcBhaCalib::readInFromFile | ( | ) | [private] |
Definition at line 627 of file EmcBhaCalib.cxx.
References calibUtil::ERROR, EmcBhaCalibData::getMatrixM(), genRecEmupikp::i, Bes_Common::INFO, m_dirHitsCut, m_fileDir, m_fileExt, m_myCalibData, m_nrNonZeros, m_nrXtalsEnoughHits, m_runNumberFile, msgSvc(), num, EmcLSSMatrix::num_nonZeros(), EmcBhaCalibData::nXtals(), EmcBhaCalibData::nXtalsHit(), EmcBhaCalibData::readIn(), deljobs::string, and EmcBhaCalibData::xtalHitsDir().
Referenced by execute().
00627 { 00628 00629 MsgStream log(msgSvc(), name()); 00630 00631 log << MSG::INFO << " Read Bhabha calibration data from files !" 00632 << endreq; 00633 00634 std::string runNumberFile = m_fileDir; 00635 runNumberFile += m_fileExt; 00636 runNumberFile += m_runNumberFile; 00637 00638 bool success = false; 00639 log << MSG::INFO << " Get file names from input file : " 00640 << runNumberFile 00641 << endreq; 00642 00643 std::string vectorFile; 00644 std::string matrixFile; 00645 00646 log << MSG::INFO << "Runnumber non-zeros xtals" 00647 << endreq; 00648 00649 ifstream datafile(runNumberFile.c_str()); 00650 00651 //cout<<"datafile="<<runNumberFile.c_str()<<""<<datafile.good()<<endl; 00652 00653 if (datafile.good() > 0 ) { 00654 00655 while( !datafile.eof() ) { 00656 00657 ifstream vectorIn; 00658 ifstream matrixIn; 00659 00660 std::string num; 00661 datafile >> num; 00662 00663 if ( num.length() > 0 ) { 00664 00665 vectorFile = m_fileDir; 00666 vectorFile += m_fileExt; 00667 vectorFile += num; 00668 vectorFile += "CalibVector.dat"; 00669 00670 matrixFile = m_fileDir; 00671 matrixFile += m_fileExt; 00672 matrixFile += num; 00673 matrixFile += "CalibMatrix.dat"; 00674 00675 vectorIn.open(vectorFile.c_str()); 00676 matrixIn.open(matrixFile.c_str()); 00677 00678 if ( vectorIn.good() > 0 && matrixIn.good() > 0 ) { 00679 00680 m_myCalibData->readIn(matrixIn,vectorIn); 00681 00682 log << MSG::INFO << num 00683 << " " 00684 << m_myCalibData->getMatrixM()->num_nonZeros() 00685 << " " 00686 << m_myCalibData->nXtalsHit() 00687 << endreq; 00688 00689 success = true; 00690 00691 } 00692 else { 00693 if ( !vectorIn ) 00694 log << MSG::ERROR << " ERROR: Vector input file " 00695 << vectorFile 00696 << " doesn't exist !" << endreq; 00697 if ( !matrixIn ) 00698 log << MSG::ERROR << " ERROR: Matrix input file " 00699 << matrixFile 00700 << " doesn't exist !" << endreq; 00701 } 00702 vectorIn.close(); 00703 matrixIn.close(); 00704 } 00705 } 00706 } 00707 00708 if ( success == true) { 00709 00710 for (int i=0; i < m_myCalibData->nXtals(); i++) { 00711 00712 if ( m_myCalibData->xtalHitsDir(i) > m_dirHitsCut ) 00713 { 00714 m_nrXtalsEnoughHits++; 00715 } 00716 } 00717 m_nrNonZeros = m_myCalibData->getMatrixM()->num_nonZeros(); 00718 log << MSG::INFO<< " Final matrix : " 00719 << "Number of non zero elements: " << m_nrNonZeros 00720 << " " 00721 << m_myCalibData->nXtalsHit() 00722 << endreq; 00723 00724 } 00725 00726 00727 return success; 00728 }
| bool EmcBhaCalib::solveEqua | ( | ) | [private] |
Definition at line 381 of file EmcBhaCalib.cxx.
References EmcLSSMatrix::column(), dsdcg_(), dsdgmr_(), calibUtil::ERROR, EmcBhaCalibData::getMatrixM(), EmcBhaCalibData::getVectorR(), m_calibConst, m_convCrit, m_equationSolver, m_myCalibData, m_nrNonZeros, EmcLSSMatrix::matrix(), msgSvc(), EmcBhaCalibData::nXtals(), EmcBhaCalibData::nXtalsHit(), EmcLSSMatrix::row(), and VERBOSE.
Referenced by execute().
00381 { 00382 00383 MsgStream log(msgSvc(), name()); 00384 //----------------------------------------------------- 00385 // solve system of equations with SLAP package 00386 //----------------------------------------------------- 00387 log << MSG::VERBOSE << " Solve Matrix equation with method " 00388 << m_equationSolver 00389 << endl 00390 << "Xtals hit: " << m_myCalibData->nXtalsHit() 00391 << endl 00392 << "Nr of non Zeros: " << m_nrNonZeros << endreq; 00393 00394 //input parameters for SLAP 00395 long int isym = 1; //matrix M is symmetric 00396 long int nsave = 20; 00397 long int itol = 0; //type of convergence criterion 00398 double convCriterion = m_convCrit; //convergence crtiterion 00399 long int maxIter = 1000; //maximum number of iterations 00400 long int errUnit = 6; //unit on which to write error estimation 00401 // at each iteration 00402 00403 //working arrays needed by SLAP 00404 long int lengthDoubleWork; 00405 double* doubleWork; 00406 long int lengthIntWork; 00407 long int* intWork; 00408 00409 //output parameters of slap 00410 long int iters=0; //number of iterations required 00411 double errorv=1000; //error estimate 00412 long int errorFlag=9999; 00413 00414 //sparse Ax=b solver. 00415 //uses the generalized minimum residual method 00416 //The preconditioner is diagonal scaling. 00417 if ( m_equationSolver == "GMR" ) { 00418 00419 cout<<m_equationSolver<<endl; 00420 00421 cout<<"errorFlag="<<errorFlag<<endl; 00422 00423 //workspaces 00424 lengthDoubleWork = (1 + m_myCalibData->nXtals()*(nsave+7) 00425 + nsave*(nsave+3)) 00426 + 50000; 00427 doubleWork = new double[lengthDoubleWork]; 00428 lengthIntWork = 50; 00429 intWork = new long int [lengthIntWork]; 00430 00431 dsdgmr_ ( &(m_myCalibData->nXtalsHit()), 00432 m_myCalibData->getVectorR(), 00433 m_calibConst, 00434 &m_nrNonZeros, 00435 m_myCalibData->getMatrixM()->row(), 00436 m_myCalibData->getMatrixM()->column(), 00437 m_myCalibData->getMatrixM()->matrix(), 00438 &isym, 00439 &nsave, 00440 &itol, 00441 &convCriterion, 00442 &maxIter, 00443 &iters, 00444 &errorv, 00445 &errorFlag, 00446 &errUnit, 00447 doubleWork, 00448 &lengthDoubleWork, 00449 intWork, 00450 &lengthIntWork 00451 ); 00452 cout<<"errorFlag="<<errorFlag<<endl; 00453 00454 00455 log << MSG::VERBOSE << " Error flag: " << errorFlag << endreq; 00456 if ( errorFlag < 0 ) errorFlag = labs(errorFlag) + 2; 00457 switch ( errorFlag ) { 00458 case 0: log << MSG::VERBOSE << " all went well" 00459 << endreq; break; 00460 case 1: log << MSG::ERROR << " insufficient storage allocated for _doubleWork or _intWork" 00461 << endreq; break; 00462 case 2: log << MSG::ERROR << " method failed to reduce norm of current residual" 00463 << endreq; break; 00464 case 3: log << MSG::ERROR << " insufficient length for _doubleWork" 00465 << endreq; break; 00466 case 4: log << MSG::ERROR << " inconsistent _itol and values" 00467 << endreq; break; 00468 default: log << MSG::ERROR << " unknown flag" << endreq; 00469 } 00470 log << MSG::VERBOSE << " Integer workspace used = " << intWork[8] << endl 00471 << " Double workspace used = " << intWork[9] << endreq; 00472 } 00473 00474 //Routine to solve a symmetric positive definite linear 00475 //system Ax = b using the Preconditioned Conjugate 00476 //Gradient method. The preconditioner is diagonal scaling. 00477 if ( m_equationSolver == "PCG" ) { 00478 00479 cout<<m_equationSolver<<endl; 00480 00481 itol = 1; 00482 00483 //workspaces 00484 lengthDoubleWork = 5 *m_myCalibData->nXtals() + 50000; 00485 doubleWork = new double[lengthDoubleWork]; 00486 lengthIntWork = 50; 00487 intWork = new long int [lengthIntWork]; 00488 00489 dsdcg_( &(m_myCalibData->nXtalsHit()), 00490 m_myCalibData->getVectorR(), 00491 m_calibConst, 00492 &m_nrNonZeros, 00493 m_myCalibData->getMatrixM()->row(), 00494 m_myCalibData->getMatrixM()->column(), 00495 m_myCalibData->getMatrixM()->matrix(), 00496 &isym, 00497 &itol, 00498 &convCriterion, 00499 &maxIter, 00500 &iters, 00501 &errorv, 00502 &errorFlag, 00503 &errUnit, 00504 doubleWork, 00505 &lengthDoubleWork, 00506 intWork, 00507 &lengthIntWork 00508 ); 00509 00510 switch ( errorFlag ) { 00511 case 0: log << MSG::VERBOSE << "all went well" << endreq; break; 00512 case 1: log << MSG::ERROR << " insufficient storage allocated for WORK or IWORK" 00513 << endreq; break; 00514 case 2: log << MSG::ERROR << " method failed to to converge in maxIter steps." 00515 << endreq; break; 00516 case 3:log << MSG::ERROR << " Error in user input. Check input value of _nXtal,_itol." 00517 << endreq; break; 00518 case 4:log << MSG::ERROR << " User error tolerance set too tight. " 00519 << "Reset to 500.0*D1MACH(3). Iteration proceeded." 00520 << endreq; break; 00521 case 5:log << MSG::ERROR << " Preconditioning matrix, M, is not Positive Definite. " 00522 << endreq; break; 00523 case 6: log << MSG::ERROR << " Matrix A is not Positive Definite." 00524 << endreq; break; 00525 default: log << MSG::ERROR << " unknown flag" << endreq; 00526 } 00527 log << MSG::VERBOSE << " Integer workspace used = " << intWork[9] << endl 00528 << "Double workspace used = " << intWork[10] << endreq; 00529 } 00530 00531 log << MSG::VERBOSE << "------ Calibration fit statistics ------- " << endl 00532 << "maximum number of iterations =" << maxIter << endl 00533 << " number of iterations =" << iters << endl 00534 << "error estimate of error in final solution =" 00535 << errorv << endreq; 00536 00537 if ( 0 != doubleWork) delete [] doubleWork; 00538 if ( 0 != intWork) delete [] intWork; 00539 00540 if ( errorFlag != 0 ) return false; 00541 else return true; 00542 00543 return true; 00544 }
| void EmcBhaCalib::writeOut | ( | ) | [private] |
Definition at line 548 of file EmcBhaCalib.cxx.
References m_fileDir, m_fileExt, m_myCalibData, msgSvc(), deljobs::string, VERBOSE, and EmcBhaCalibData::writeOut().
Referenced by execute().
00548 { 00549 00550 ofstream vectorOut; 00551 std::string vectorFile = m_fileDir; 00552 vectorFile += m_fileExt; 00553 vectorFile += "allCalibVector.dat"; 00554 vectorOut.open(vectorFile.c_str()); 00555 00556 ofstream matrixOut; 00557 std::string matrixFile = m_fileDir; 00558 matrixFile += m_fileExt; 00559 matrixFile += "allCalibMatrix.dat"; 00560 matrixOut.open(matrixFile.c_str()); 00561 00562 MsgStream log(msgSvc(), name()); 00563 00564 log << MSG::VERBOSE << " Write out files " 00565 << vectorFile << " " 00566 << matrixFile 00567 << endreq; 00568 00569 m_myCalibData->writeOut(matrixOut,vectorOut); 00570 00571 vectorOut.close(); 00572 matrixOut.close(); 00573 00574 }
| void EmcBhaCalib::writeOutConst | ( | ) | [private] |
Definition at line 578 of file EmcBhaCalib.cxx.
References genRecEmupikp::i, m_absoluteConstants, m_calibConstUnred, m_dirHitsCut, m_fileDir, m_fileExt, m_myCalibData, m_oldConstants, EmcBhaCalibData::nXtals(), EmcBhaCalibData::nXtalsHit(), deljobs::string, EmcBhaCalibData::xtalHitsDir(), and EmcBhaCalibData::xtalInd().
Referenced by execute().
00578 { 00579 00580 ofstream constOut; 00581 00582 std::string constFile = m_fileDir; 00583 constFile += m_fileExt; 00584 constFile += "CalibConst.dat"; 00585 00586 constOut.open(constFile.c_str()); 00587 00588 constOut << "#crystalIndex relative-constant absolute-constant" << endl; 00589 00590 int chanIndex; 00591 00592 //output constants to file 00593 for ( int i=0; i < m_myCalibData->nXtalsHit(); i++) { 00594 00595 chanIndex=m_myCalibData->xtalInd(i); //xtalInd(i) array of xtal indices 00596 00597 //---- get the new absolute constant ---- 00598 //set it to 0 if we have not enough hits -> we'll keep the old constant 00599 if ( m_myCalibData->xtalHitsDir(i) < m_dirHitsCut ) 00600 m_absoluteConstants[chanIndex] = m_oldConstants[chanIndex]; 00601 else 00602 m_absoluteConstants[chanIndex] = 00603 m_oldConstants[chanIndex] * m_calibConstUnred[chanIndex]; 00604 00605 } 00606 00607 //output constants to file 00608 for ( int i=0; i < m_myCalibData->nXtals(); i++) { 00609 00610 long Xtal_Index = i; 00611 if(m_calibConstUnred[Xtal_Index]>0){ 00612 constOut << Xtal_Index << " " 00613 << m_calibConstUnred[Xtal_Index] << " " 00614 << m_oldConstants[Xtal_Index] << " " 00615 << m_absoluteConstants[Xtal_Index] 00616 << endl; 00617 00618 00619 } 00620 } 00621 constOut.close(); 00622 00623 }
NTuple::Item<float> EmcBhaCalib::calibConst [private] |
double EmcBhaCalib::DigiCalibConst [private] |
Definition at line 190 of file EmcBhaCalib.h.
NTuple::Item<float> EmcBhaCalib::err [private] |
NTuple::Item<float> EmcBhaCalib::gi0 [private] |
NTuple::Item<long> EmcBhaCalib::ixtal [private] |
double* EmcBhaCalib::m_absoluteConstants [private] |
Definition at line 164 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), initCalibConst(), prepareConstants(), writeOutConst(), and ~EmcBhaCalib().
bool EmcBhaCalib::m_askForMatrixInversion [private] |
double* EmcBhaCalib::m_calibConst [private] |
Definition at line 158 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), execute(), initCalibConst(), solveEqua(), and ~EmcBhaCalib().
double* EmcBhaCalib::m_calibConstUnred [private] |
Definition at line 161 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), execute(), initCalibConst(), ntupleOut(), prepareConstants(), writeOutConst(), and ~EmcBhaCalib().
bool EmcBhaCalib::m_calibrationSuccessful [private] |
double EmcBhaCalib::m_convCrit [private] |
int EmcBhaCalib::m_dirHitsCut [private] |
Definition at line 123 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), execute(), prepareConstants(), printInfo(), readInFromFile(), and writeOutConst().
IEmcCalibConstSvc* EmcBhaCalib::m_emcCalibConstSvc [private] |
Definition at line 192 of file EmcBhaCalib.h.
Referenced by digiConstCor(), initCalibConst(), initialize(), and prepareConstants().
std::string EmcBhaCalib::m_equationSolver [private] |
std::string EmcBhaCalib::m_fileDir [private] |
Definition at line 148 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), printInfo(), readInFromFile(), writeOut(), and writeOutConst().
std::string EmcBhaCalib::m_fileExt [private] |
Definition at line 145 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), printInfo(), readInFromFile(), writeOut(), and writeOutConst().
bool EmcBhaCalib::m_fitPolynom [private] |
int EmcBhaCalib::m_MsgFlag [private] |
EmcBhaCalibData* EmcBhaCalib::m_myCalibData [private] |
Definition at line 151 of file EmcBhaCalib.h.
Referenced by digiConstCor(), execute(), initCalibConst(), initialize(), ntupleOut(), prepareConstants(), printInfo(), readInFromFile(), solveEqua(), writeOut(), writeOutConst(), and ~EmcBhaCalib().
long int EmcBhaCalib::m_nrNonZeros [private] |
int EmcBhaCalib::m_nrXtalsEnoughHits [private] |
double* EmcBhaCalib::m_oldConstants [private] |
Definition at line 167 of file EmcBhaCalib.h.
Referenced by digiConstCor(), EmcBhaCalib(), initCalibConst(), ntupleOut(), prepareConstants(), writeOutConst(), and ~EmcBhaCalib().
double EmcBhaCalib::m_peakCor[6240] [private] |
Definition at line 120 of file EmcBhaCalib.h.
bool EmcBhaCalib::m_readDataFromDB [private] |
std::string EmcBhaCalib::m_runNumberFile [private] |
NTuple::Tuple* EmcBhaCalib::m_tuple1 [private] |
Definition at line 183 of file EmcBhaCalib.h.
Referenced by EmcBhaCalib(), initialize(), and ntupleOut().
bool EmcBhaCalib::m_writeToFile [private] |
NTuple::Item<long> EmcBhaCalib::nhitxtal [private] |
1.4.7