CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/CalibMuon/DTCalibration/src/DTMeanTimerFitter.cc

Go to the documentation of this file.
00001 
00002 /*
00003  *  See header file for a description of this class.
00004  *
00005  *  $Date: 2009/02/26 09:46:57 $
00006  *  $Revision: 1.8 $
00007  *  \author S. Bolognesi - INFN Torino
00008  */
00009 
00010 #include "CalibMuon/DTCalibration/interface/DTMeanTimerFitter.h"
00011 //#include "CalibMuon/DTCalibration/plugins/vDriftHistos.h"
00012 #include "CalibMuon/DTCalibration/interface/vDriftHistos.h"
00013 
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 
00016 #include <iostream>
00017 
00018 #include "TFile.h"
00019 #include "TF1.h"
00020 
00021 using namespace std;
00022 
00023 DTMeanTimerFitter::DTMeanTimerFitter(TFile *file) : hInputFile(file), theVerbosityLevel(0) {
00024   //hInputFile = new TFile("DTTMaxHistos.root", "READ");
00025   hDebugFile = new TFile("DTMeanTimerFitter.root", "RECREATE");
00026 }
00027 
00028 DTMeanTimerFitter::~DTMeanTimerFitter() {
00029   hDebugFile->Close();
00030 }
00031 
00032 vector<float> DTMeanTimerFitter::evaluateVDriftAndReso (TString N) {
00033   
00034   // Retrieve histogram sets
00035   hTMaxCell * histos   = new hTMaxCell(N, hInputFile);
00036   vector<float> vDriftAndReso;
00037 
00038   // Check that the histo for this cell exists
00039   if(histos->hTmax123 != 0) {
00040     vector<TH1F*> hTMax;  // histograms for <T_max> calculation
00041     vector <TH1F*> hT0;   // histograms for T0 evaluation
00042     hTMax.push_back(histos->hTmax123); 
00043     hTMax.push_back(histos->hTmax124s72);
00044     hTMax.push_back(histos->hTmax124s78);
00045     hTMax.push_back(histos->hTmax134s72);
00046     hTMax.push_back(histos->hTmax134s78);
00047     hTMax.push_back(histos->hTmax234);
00048 
00049     hT0.push_back(histos->hTmax_3t0);
00050     hT0.push_back(histos->hTmax_2t0);
00051     hT0.push_back(histos->hTmax_t0);
00052     hT0.push_back(histos->hTmax_0);
00053 
00054     vector<Double_t> factor; // factor relating the width of the Tmax distribution 
00055                              // and the cell resolution 
00056     factor.push_back(sqrt(2./3.)); // hTmax123
00057     factor.push_back(sqrt(2./7.)); // hTmax124s72
00058     factor.push_back(sqrt(8./7.)); // hTmax124s78
00059     factor.push_back(sqrt(2./7.)); // hTmax134s72
00060     factor.push_back(sqrt(8./7.)); // hTmax134s78
00061     factor.push_back(sqrt(2./3.)); // hTmax234
00062 
00063 
00064     // Retrieve the gaussian mean and sigma for each TMax histogram    
00065     vector<Double_t> mean;
00066     vector<Double_t> sigma; 
00067     vector<Double_t> count;  //number of entries
00068 
00069     for(vector<TH1F*>::const_iterator ith = hTMax.begin();
00070         ith != hTMax.end(); ith++) {
00071       TF1 *funct = fitTMax(*ith);
00072       if(!funct){
00073         edm::LogError("DTMeanTimerFitter") << "Error when fitting TMax..histogram name" << (*ith)->GetName();
00074         // return empty or -1 filled vector?
00075         vector<float> defvec(6,-1);
00076         return defvec;
00077       }                 
00078       hDebugFile->cd();
00079       (*ith)->Write();
00080 
00081     // Get mean, sigma and number of entries of each histogram
00082       mean.push_back(funct->GetParameter(1));
00083       sigma.push_back(funct->GetParameter(2)); 
00084       count.push_back((*ith)->GetEntries());  
00085     } 
00086           
00087     Double_t tMaxMean=0.;
00088     Double_t wTMaxSum=0.;
00089     Double_t sigmaT=0.;
00090     Double_t wSigmaSum = 0.;
00091   
00092     //calculate total mean and sigma
00093     for(int i=0; i<=5; i++) {
00094       if(count[i]<200) continue;
00095       tMaxMean  += mean[i]*(count[i]/(sigma[i]*sigma[i]));
00096       wTMaxSum  += count[i]/(sigma[i]*sigma[i]);
00097       sigmaT    += count[i]*factor[i]*sigma[i];
00098       wSigmaSum += count[i];
00099       // cout << "TMaxMean "<<i<<": "<< mean[i] << " entries: " << count[i] 
00100       // << " sigma: " << sigma[i] 
00101       // << " weight: " << (count[i]/(sigma[i]*sigma[i])) << endl; 
00102     }
00103     if((!wTMaxSum)||(!wSigmaSum)){
00104       edm::LogError("DTMeanTimerFitter") << "Error zero sum of weights..returning default";
00105       vector<float> defvec(6,-1);
00106       return defvec;    
00107     }
00108 
00109     tMaxMean /= wTMaxSum;
00110     sigmaT /= wSigmaSum;
00111 
00112     //calculate v_drift and resolution
00113     Double_t vDrift = 2.1 / tMaxMean; //2.1 is the half cell length in cm
00114     Double_t reso = vDrift * sigmaT;
00115     vDriftAndReso.push_back(vDrift);
00116     vDriftAndReso.push_back(reso);
00117     //if(theVerbosityLevel >= 1)
00118     edm::LogVerbatim("DTMeanTimerFitter") << " final TMaxMean=" << tMaxMean << " sigma= "  << sigmaT 
00119            << " v_d and reso: " << vDrift << " " << reso << endl;
00120 
00121     // Order t0 histogram by number of entries (choose histograms with higher nr. of entries)
00122     map<Double_t,TH1F*> hEntries;    
00123     for(vector<TH1F*>::const_iterator ith = hT0.begin();
00124         ith != hT0.end(); ith++) {
00125       hEntries[(*ith)->GetEntries()] = (*ith);
00126     } 
00127 
00128     // add at the end of hT0 the two hists with the higher number of entries 
00129     int counter = 0;
00130     for(map<Double_t,TH1F*>::reverse_iterator iter = hEntries.rbegin();
00131         iter != hEntries.rend(); iter++) {
00132       counter++;
00133       if (counter==1) hT0.push_back(iter->second); 
00134       else if (counter==2) {hT0.push_back(iter->second); break;} 
00135     }
00136     
00137     // Retrieve the gaussian mean and sigma of histograms for Delta(t0) evaluation   
00138     vector<Double_t> meanT0;
00139     vector<Double_t> sigmaT0; 
00140     vector<Double_t> countT0;
00141 
00142     for(vector<TH1F*>::const_iterator ith = hT0.begin();
00143         ith != hT0.end(); ith++) {
00144       try{
00145         (*ith)->Fit("gaus");
00146       } catch(...){
00147         edm::LogError("DTMeanTimerFitter") << "Exception when fitting T0..histogram " << (*ith)->GetName();    
00148         // return empty or -1 filled vector?
00149         vector<float> defvec(6,-1);
00150         return defvec;
00151       } 
00152       TF1 *f1 = (*ith)->GetFunction("gaus");
00153       // Get mean, sigma and number of entries of the  histograms
00154       meanT0.push_back(f1->GetParameter(1));
00155       sigmaT0.push_back(f1->GetParameter(2));
00156       countT0.push_back((*ith)->GetEntries());
00157     }
00158     //calculate Delta(t0)
00159     if(hT0.size() != 6) { // check if you have all the t0 hists
00160       edm::LogVerbatim("DTMeanTimerFitter") << "t0 histograms = " << hT0.size();
00161       for(int i=1; i<=4;i++) {
00162         vDriftAndReso.push_back(-1);
00163       }
00164       return vDriftAndReso;
00165     }
00166     
00167     for(int it0=0; it0<=2; it0++) {      
00168       if((countT0[it0] > 200) && (countT0[it0+1] > 200)) {
00169         Double_t deltaT0 = meanT0[it0] - meanT0[it0+1]; 
00170         vDriftAndReso.push_back(deltaT0);
00171       }  
00172       else
00173         vDriftAndReso.push_back(999.);
00174     }
00175     //deltat0 using hists with max nr. of entries
00176     if((countT0[4] > 200) && (countT0[5] > 200)) {
00177       Double_t t0Diff = histos->GetT0Factor(hT0[4]) - histos->GetT0Factor(hT0[5]);
00178       Double_t deltaT0MaxEntries =  (meanT0[4] - meanT0[5])/ t0Diff;
00179       vDriftAndReso.push_back(deltaT0MaxEntries);
00180     }
00181     else
00182       vDriftAndReso.push_back(999.);
00183   }
00184   else {
00185     for(int i=1; i<=6; i++) { 
00186       // 0=vdrift, 1=reso,  2=(3deltat0-2deltat0), 3=(2deltat0-1deltat0), 
00187       // 4=(1deltat0-0deltat0), 5=deltat0 from hists with max entries,
00188       vDriftAndReso.push_back(-1);
00189     }
00190   }
00191   return vDriftAndReso;
00192 }
00193 
00194 
00195 
00196 TF1* DTMeanTimerFitter::fitTMax(TH1F* histo){
00197  // Find distribution peak and fit range
00198       Double_t peak = (((((histo->GetXaxis())->GetXmax())-((histo->GetXaxis())->GetXmin()))/histo->GetNbinsX())*
00199                        (histo->GetMaximumBin()))+((histo->GetXaxis())->GetXmin());
00200       //if(theVerbosityLevel >= 1)
00201       LogDebug("DTMeanTimerFitter") <<"Peak "<<peak<<" : "<<"xmax "<<((histo->GetXaxis())->GetXmax())
00202             <<"            xmin "<<((histo->GetXaxis())->GetXmin())
00203             <<"            nbin "<<histo->GetNbinsX()
00204             <<"            bin with max "<<(histo->GetMaximumBin());
00205       Double_t range = 2.*histo->GetRMS(); 
00206 
00207       // Fit each Tmax histogram with a Gaussian in a restricted interval
00208       TF1 *rGaus = new TF1("rGaus","gaus",peak-range,peak+range);
00209       rGaus->SetMarkerSize();// just silence gcc complainining about unused var
00210       try{      
00211         histo->Fit("rGaus","R");
00212       } catch(...){
00213         edm::LogError("DTMeanTimerFitter") << "Exception when fitting TMax..histogram " << histo->GetName()
00214              << "   setting return function pointer to zero";   
00215         return 0;
00216       } 
00217       TF1 *f1 = histo->GetFunction("rGaus");
00218       return f1;
00219  }