StRoot: StEmcCalibrationMaker/StEmcCalibMaker.cxx Source File

00001 #include "StEmcCalibMaker.h"
00002 
00003 #include <TFile.h>
00004 #include <TROOT.h>
00005 #include <TF1.h>
00006 
00007 #include <StEventTypes.h>
00008 #include <StEvent.h>
00009 #include <StEmcUtil/geometry/StEmcGeom.h>
00010 #include <StEmcUtil/voltageCalib/VoltCalibrator.h>
00011 #include <StDbLib/StDbManager.hh>
00012 #include <StDbLib/StDbTable.h>
00013 #include <StDbLib/StDbConfigNode.hh>
00014 #include <tables/St_emcCalib_Table.h>
00015 
00016 ClassImp(StEmcCalibMaker)
00017 
00018 //_____________________________________________________________________________
00019 StEmcCalibMaker::StEmcCalibMaker(const Char_t *name)
00020     :StMaker(name)
00021     {
00022   mCalib = NULL;
00023   mAccept = NULL;
00024   mSpec = NULL;
00025   mDate = 20330101;
00026   mNMinTracks = 0;
00027   mNMaxTracks = 100000;
00028   mDebug = false;
00029   setDetector(1);
00030   setRange(200);
00031   mZDCMin = -1000000;
00032   mZDCMax =  1000000;
00033   mCTBMin = -1000000;
00034   mCTBMax =  1000000;
00035   mSpecName = "mSpec";
00036   mAcceptName = "mAccept";
00037   setFile("spec.root");
00038 }
00039 //_____________________________________________________________________________
00040 StEmcCalibMaker::~StEmcCalibMaker()
00041 {
00042 }
00043 //_____________________________________________________________________________
00044 Int_t StEmcCalibMaker::Init()
00045 {
00046   mSpec = new TH2F(mSpecName.Data(),"Equal Spectra",mNChannel,+0.5,(float)mNChannel+0.5,(int)mRange,-0.5,mRange-0.5);
00047   mAccept = new TH1F(mAcceptName.Data(),"Accepted events",20,-0.5,19.5);
00048   return StMaker::Init();
00049 }
00050 //_____________________________________________________________________________
00051 void StEmcCalibMaker::Clear(Option_t *option)              
00052 {
00053 }
00054 //_____________________________________________________________________________
00055 Int_t StEmcCalibMaker::Make()
00056 {    
00057   return kStOK;
00058 }
00059 //_____________________________________________________________________________
00060 bool StEmcCalibMaker::accept()
00061 {  
00062   mAccept->Fill(0);
00063         if(!mCalib) mCalib = (StEmcCalibrationMaker*)GetMaker("Calib");
00064   if(!mCalib) return false;
00065   mAccept->Fill(1);
00066         
00067         if(!mCalib->hasDetector(getDetector())) return false;
00068   mAccept->Fill(2);
00069   
00070   if(mCalib->getNTracks()<mNMinTracks) return false;
00071   mAccept->Fill(3);
00072         
00073   if(mCalib->getNTracks()>mNMaxTracks) return false;
00074   mAccept->Fill(4);
00075   
00076   if(GetDate()<mDate) { mDate = GetDate(); mTime = GetTime();}
00077   
00078   if(!mCalib->isTriggerOk()) return false;
00079   mAccept->Fill(5);
00080         
00081         long CTB = mCalib->getCTBSum();
00082         long ZDC = mCalib->getZDCSum();
00083         
00084         if(ZDC<mZDCMin || ZDC > mZDCMax) return false;
00085   mAccept->Fill(6);
00086         
00087         if(CTB<mCTBMin || CTB > mCTBMax) return false;
00088   mAccept->Fill(7);
00089         
00090   
00091   mNEvents++;
00092   if(mNEvents%100==0) saveHist(mFileName.Data());
00093   
00094         //cout <<"Processing event number "<<mNEvents<<" for maker "<<GetName()<<endl;
00095         return true;
00096 }
00097 //_____________________________________________________________________________
00098 Int_t StEmcCalibMaker::Finish()
00099 {
00100   saveHist(mFileName.Data());
00101   return StMaker::Finish();
00102 }
00103 //_____________________________________________________________________________
00104 void StEmcCalibMaker::saveHist(const Char_t *file)
00105 {
00106   cout << "Saving spec to file " << file << endl;
00107   TFile *f = new TFile(file, "RECREATE");
00108   mSpec->Write();
00109   f->Close();
00110   delete f;
00111   return ;
00112 }
00113 //_____________________________________________________________________________
00114 void StEmcCalibMaker::loadHist(const Char_t *file)
00115 {
00116         TFile *f = new TFile(file);
00117   TString N = mSpec->GetName();
00118   N+=";1";
00119   TH2F *h = (TH2F*)f->Get(N.Data());
00120   if(h){ mSpec->Reset(); mSpec->Add(h,1);}
00121   f->Close();
00122   delete f;
00123   return ;
00124 }
00125 //_____________________________________________________________________________
00126 void StEmcCalibMaker::addHist(const Char_t *file)
00127 {
00128   TFile *f = new TFile(file);
00129   TString N = mSpec->GetName();
00130   N+=";1";
00131   TH2F *h = (TH2F*)f->Get(N.Data());
00132   if(h){ mSpec->Add(h,1);}
00133   f->Close();
00134   delete f;
00135   return ;
00136 }
00137 //_____________________________________________________________________________
00138 float StEmcCalibMaker::getTimeInterval(int refDate, int refTime)
00139 {
00140   struct tm t0;
00141   struct tm t1;
00142     
00143   Int_t year  = (Int_t)(mDate/10000);
00144   t0.tm_year=year-1900;
00145   Int_t month = (Int_t)(mDate-year*10000)/100; 
00146   t0.tm_mon=month-1;
00147   Int_t day   = (Int_t)(mDate-year*10000-month*100);
00148   t0.tm_mday=day;
00149   Int_t hour  = (Int_t)(mTime/10000);
00150   t0.tm_hour=hour;
00151   Int_t minute= (Int_t)(mTime-hour*10000)/100;
00152   t0.tm_min=minute;
00153   Int_t second= (Int_t)(mTime-hour*10000-minute*100);
00154   t0.tm_sec=second;
00155   t0.tm_isdst = -1;
00156   
00157   year  = (Int_t)(refDate/10000);
00158   t1.tm_year=year-1900;
00159   month = (Int_t)(refDate-year*10000)/100; 
00160   t1.tm_mon=month-1;
00161   day   = (Int_t)(refDate-year*10000-month*100);
00162   t1.tm_mday=day;
00163   hour  = (Int_t)(refTime/10000);
00164   t1.tm_hour=hour;
00165   minute= (Int_t)(refTime-hour*10000)/100;
00166   t1.tm_min=minute;
00167   second= (Int_t)(refTime-hour*10000-minute*100);
00168   t1.tm_sec=second;
00169   t1.tm_isdst = -1;
00170 
00171   time_t ttime0=mktime(&t0);
00172   time_t ttime1=mktime(&t1);
00173   double diff=difftime(ttime0,ttime1);
00174   
00175   //cout <<"startdate = "<<startDate<<" now = "<<d<<"  dt = "<<diff/60<<endl;
00176   //cout <<diff<<"  "<<diff/60.<<endl;
00177   float ddd=fabs(diff/3600.);
00178   return ddd;
00179 }
00180 //_____________________________________________________________________________
00181 void StEmcCalibMaker::getMeanAndRms(TH1D* tmp,float amin,float amax,
00182                                    float* m,float* r)
00183 {  
00184   float mean=0,rms=0,sum=0;
00185   if(tmp->Integral()>0)
00186   {
00187           Int_t bin0 = tmp->FindBin(amin);
00188           Int_t bin1 = tmp->FindBin(amax);
00189     for(Int_t j=bin0;j<bin1;j++)
00190     {
00191       float w = tmp->GetBinContent(j);
00192                   float x = tmp->GetBinCenter(j);
00193       mean += w*x;
00194       sum  += w;
00195       rms  += x*x*w;
00196     }
00197     if(sum>0)
00198     {
00199       mean /= sum;
00200       rms   = ::sqrt(rms/sum-mean*mean);
00201     }
00202   }
00203   *m    = mean;
00204   *r    = rms;
00205   return;
00206 
00207 }
00208 //_____________________________________________________________________________
00209 void StEmcCalibMaker::getLogMeanAndRms(TH1D* tmp,float amin,float amax,
00210                                       float* m,float* r)
00211 {
00212   float mean=0,rms=0,sum=0;
00213   
00214   if(tmp->Integral()>0)
00215   {
00216           Int_t bin0 = tmp->FindBin(amin);
00217           Int_t bin1 = tmp->FindBin(amax);
00218     for(Int_t j=bin0;j<bin1;j++)
00219     {
00220       float w = tmp->GetBinContent(j);
00221                         if(w==1) w=1.5;
00222                   if(w>0) w = ::log(w);
00223                   float x = tmp->GetBinCenter(j);
00224       mean += w*x;
00225       sum  += w;
00226       rms  += x*x*w;
00227     }
00228     if(sum>0)
00229     {
00230             mean /= sum;
00231       rms   = ::sqrt(rms/sum-mean*mean);
00232     }
00233   }
00234   *m    = mean;
00235   *r    = rms;
00236   return;
00237 
00238 }
00239 StEmcGeom* StEmcCalibMaker::getGeom()
00240 {
00241         return StEmcGeom::instance(getDetector());
00242 }
00243 //_____________________________________________________________________________
00244 void StEmcCalibMaker::calcVoltages(TH1F *gain, char* refFile, char* voltageInput, char* voltageOutput)
00245 {
00246   char GainShift[]="HVGainShift.dat";
00247         char line[300];
00248         
00249         ofstream out(GainShift);
00250         ifstream inp(voltageInput);
00251         //TH1F *h = new TH1F("Gain","Current Gain Distribution",100,0,0.02);
00252         do
00253         {
00254           int a,b,c,d,e,f,id;
00255                 float v;
00256                 inp>>a>>id>>b>>c>>d>>e>>f>>v;           
00257     float gainShift = gain->GetBinContent(id);
00258                 sprintf(line,"%2d %5d  %5d  %5d  %5d  %5d  %5d  %7.5f",a,id,b,c,d,e,f,gainShift);
00259                 out <<line<<endl;
00260                 cout <<line<<endl;
00261                 
00262         } while(!inp.eof());    
00263   out.close();
00264         
00265         cout <<"Calculating voltages\n";
00266         
00267   // Instantiate the voltage calibrator class
00268   VoltCalibrator vc;
00269   // Select the file that contains the Gain vs HV data
00270   vc.setRefFile(refFile);
00271   // Select the file that contains the relative gains to be
00272   // be achieved
00273   vc.setGainFile(GainShift);
00274   // Select the file that containe the current voltages
00275   vc.setVoltInputFile(voltageInput);
00276   // Select the file that will contain the new request voltages
00277   vc.setVoltOutputFile(voltageOutput);
00278   // Perform the calculation
00279   vc.process();
00280 }