CMSSW Reference Manual

00001 /* ******************************************
00002  * ZIterativeAlgorithmWithFit.cc
00003  *
00004  * 
00005  * Paolo Meridiani 06/07/2005
00006  * Rewritten for CMSSW 04/06/2007
00007  ********************************************/
00008 
00009 #include "Calibration/Tools/interface/ZIterativeAlgorithmWithFit.h"
00010 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
00011 #include "Calibration/Tools/interface/EcalIndexingTools.h"
00012 
00013 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00014 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00015 #include "DataFormats/TrackReco/interface/Track.h"
00016 
00017 #include <TMath.h>
00018 #include <TCanvas.h> 
00019 #include "TH1.h"
00020 #include "TH2.h"
00021 #include "TF1.h"
00022 #include "TH1F.h"
00023 #include "TMinuit.h"
00024 #include "TGraphErrors.h"
00025 #include "THStack.h"
00026 #include "TLegend.h"
00027 
00028 
00029 #include <fstream>
00030 #include <iostream>
00031 #include <vector>
00032 
00033 //#include "Tools.C"
00034 
00035 //Scale and Bins for calibration factor histograms
00036 #define MIN_RESCALE -0.5
00037 #define MAX_RESCALE 0.5
00038 #define NBINS_LOWETA 100
00039 #define NBINS_HIGHETA 50
00040 
00041 const double ZIterativeAlgorithmWithFit::M_Z_=91.187;
00042 
00043 //  #if !defined(__CINT__)
00044 //  ClassImp(Electron)
00045 //  #endif
00046 
00047 ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFit()
00048 {
00049   // cout<< "ZIterativeAlgorithmWithFit::Called Construnctor" << std::endl;
00050   numberOfIterations_=10;
00051   channels_=1;
00052   totalEvents_=0;
00053   currentEvent_=0;
00054   currentIteration_=0;
00055   optimizedCoefficients_.resize(channels_);
00056   optimizedCoefficientsError_.resize(channels_);
00057   optimizedChiSquare_.resize(channels_);
00058   optimizedIterations_.resize(channels_);
00059   calib_fac_.resize(channels_);
00060   weight_sum_.resize(channels_);
00061   electrons_.resize(1);
00062   massReco_.resize(1);
00063 }
00064 
00065 ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFit(const edm::ParameterSet&  ps)
00066   //k, unsigned int events) 
00067 {
00068   //std::cout<< "ZIterativeAlgorithmWithFit::Called Constructor" << std::endl;
00069   numberOfIterations_=ps.getUntrackedParameter<unsigned int>("maxLoops",0);
00070   massMethod = ps.getUntrackedParameter<std::string>("ZCalib_InvMass","SCMass");
00071   calibType_= ps.getUntrackedParameter<std::string>("ZCalib_CalibType","RING"); 
00072 
00073   if (calibType_ == "RING")
00074     channels_ = EcalRingCalibrationTools::N_RING_TOTAL;
00075   else if (calibType_ == "MODULE")
00076 
00077 
00078     channels_ = EcalRingCalibrationTools::N_MODULES_BARREL;  
00079   else if (calibType_ == "ABS_SCALE")
00080     channels_ = 1;  
00081   else if(calibType_ == "ETA_ET_MODE")
00082     channels_ = EcalIndexingTools::getInstance()->getNumberOfChannels();
00083 
00084   std::cout << "[ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFit] channels_ = " << channels_ << std::endl;
00085 
00086   nCrystalCut_=ps.getUntrackedParameter<int>("ZCalib_nCrystalCut",-1);
00087 
00088   //Resetting currentEvent & iteration
00089   currentEvent_=0;
00090   currentIteration_=0;
00091   totalEvents_=0;
00092 
00093   optimizedCoefficients_.resize(channels_);
00094   optimizedCoefficientsError_.resize(channels_);
00095   optimizedChiSquare_.resize(channels_);
00096   optimizedIterations_.resize(channels_);
00097   calib_fac_.resize(channels_);
00098   weight_sum_.resize(channels_);
00099 
00100   //Creating and booking histograms
00101   thePlots_ = new ZIterativeAlgorithmWithFitPlots;
00102   bookHistograms();
00103 
00104   //Setting up rescaling if needed
00105   UseStatWeights_=ps.getUntrackedParameter<bool>("ZCalib_UseStatWeights",false);
00106   if (UseStatWeights_) {
00107     WeightFileName_="weights.txt";
00108     StatWeights_.resize(channels_);
00109     getStatWeights(WeightFileName_);
00110     //    Event_Weight_.resize(events);
00111   }
00112 }
00113 
00114 
00115 void ZIterativeAlgorithmWithFit::bookHistograms()
00116 {
00117   if (!thePlots_)
00118     return;
00119 
00120   for (unsigned int i2 = 0; i2 < numberOfIterations_; i2++) {
00121     for (unsigned int i1 = 0; i1 < channels_; i1++) {
00122       char histoName[200];
00123       char histoTitle[200];
00124 
00125       //WeightedRescaling factor
00126       sprintf(histoName, "WeightedRescaleFactor_channel_%d_Iteration_%d",i1, i2);
00127       sprintf(histoTitle, "WeightedRescaleFactor Channel_%d Iteration %d",i1, i2);
00128       if (i1>15 && i1<155)
00129         thePlots_->weightedRescaleFactor[i2][i1] = new TH1F(histoName, histoTitle, NBINS_LOWETA, MIN_RESCALE, MAX_RESCALE);
00130       else
00131         thePlots_->weightedRescaleFactor[i2][i1] = new TH1F(histoName, histoTitle, NBINS_HIGHETA, MIN_RESCALE, MAX_RESCALE);
00132       thePlots_->weightedRescaleFactor[i2][i1]->GetXaxis()->SetTitle("Rescale factor");
00133       thePlots_->weightedRescaleFactor[i2][i1]->GetYaxis()->SetTitle("a.u.");
00134 
00135       //UnweightedRescaling factor
00136       sprintf(histoName, "UnweightedRescaleFactor_channel_%d_Iteration_%d",i1, i2);
00137       sprintf(histoTitle, "UnweightedRescaleFactor Channel_%d Iteration %d",i1, i2);
00138       if (i1>15 && i1<155)
00139         thePlots_->unweightedRescaleFactor[i2][i1] = new TH1F(histoName, histoTitle, NBINS_LOWETA, MIN_RESCALE, MAX_RESCALE);
00140       else
00141         thePlots_->unweightedRescaleFactor[i2][i1] = new TH1F(histoName, histoTitle, NBINS_HIGHETA, MIN_RESCALE, MAX_RESCALE);
00142       thePlots_->unweightedRescaleFactor[i2][i1]->GetXaxis()->SetTitle("Rescale factor");
00143       thePlots_->unweightedRescaleFactor[i2][i1]->GetYaxis()->SetTitle("a.u.");
00144 
00145       //Weights
00146       sprintf(histoName, "Weight_channel_%d_Iteration_%d",i1, i2);
00147       sprintf(histoTitle, "Weight Channel_%d Iteration %d",i1, i2);
00148       thePlots_->weight[i2][i1] = new TH1F(histoName, histoTitle, 100, 0., 1.);
00149       thePlots_->weight[i2][i1]->GetXaxis()->SetTitle("Weight");
00150       thePlots_->weight[i2][i1]->GetYaxis()->SetTitle("a.u.");
00151     }
00152   }
00153 }
00154 
00155 void ZIterativeAlgorithmWithFit::getStatWeights(char* file) {
00156   ifstream statfile;
00157   statfile.open(file);
00158   if (!statfile) {
00159     std::cout << "ZIterativeAlgorithmWithFit::FATAL: stat weight  file " << file << " not found" << std::endl;
00160     exit(-1);
00161   }
00162   for(unsigned int i=0;i<channels_;i++) {
00163     int imod;
00164     statfile >> imod >> StatWeights_[i];
00165     //std::cout << "Read Stat Weight for module " << imod << ": " <<  StatWeights_[i] << std::endl;
00166   }
00167 }
00168 
00169 
00170 bool ZIterativeAlgorithmWithFit::resetIteration()
00171 {
00172   totalEvents_=0;
00173   currentEvent_=0;
00174   
00175   //Reset correction
00176   massReco_.clear();
00177   for (unsigned int i=0;i<channels_;i++) calib_fac_[i]=0.;
00178   for (unsigned int i=0;i<channels_;i++) weight_sum_[i]=0.;
00179 
00180   return kTRUE;
00181 }    
00182 
00183 
00184 bool ZIterativeAlgorithmWithFit::iterate()
00185 {
00186 
00187   //Found optimized coefficients
00188   for (int i=0;i<(int)channels_;i++) 
00189     { 
00190 
00191       //RP      if (weight_sum_[i]!=0. && calib_fac_[i]!=0.) {
00192       if( (nCrystalCut_ == -1) || ((!(i <=  nCrystalCut_ -1 )) &&
00193                                   !((i > (19-nCrystalCut_)) && (i <= (19+nCrystalCut_))) &&
00194                                   !((i > (39-nCrystalCut_)) && (i <= (39+nCrystalCut_))) &&
00195                                   !((i > (59-nCrystalCut_)) && (i <= (59+nCrystalCut_))) &&
00196                                   !((i > (84-nCrystalCut_)) && (i <= (84+nCrystalCut_))) &&
00197                                   !((i > (109-nCrystalCut_)) && (i <= (109+nCrystalCut_))) &&
00198                                   !((i > (129-nCrystalCut_)) && (i <= (129+nCrystalCut_))) &&
00199                                   !((i > (149-nCrystalCut_)) && (i <= (149+nCrystalCut_))) &&
00200                                   !(i > (169-nCrystalCut_))))
00201         {
00202           if (weight_sum_[i]!=0.) {
00203             //optimizedCoefficients_[i] = calib_fac_[i]/weight_sum_[i];
00204             
00205             double gausFitParameters[3], gausFitParameterErrors[3], gausFitChi2;
00206             int gausFitIterations;
00207             
00208             gausfit( (TH1F*)thePlots_->weightedRescaleFactor[currentIteration_][i], gausFitParameters, gausFitParameterErrors, 2.5, 2.5, &gausFitChi2, &gausFitIterations );
00209             
00210             float peak=gausFitParameters[1];
00211             float peakError=gausFitParameterErrors[1];
00212             float chi2 = gausFitChi2;
00213             
00214             int iters = gausFitIterations;
00215                   
00216             optimizedCoefficientsError_[i] = peakError;
00217             optimizedChiSquare_[i] = chi2;
00218             optimizedIterations_[i] = iters;
00219             
00220             if (peak >=MIN_RESCALE && peak <= MAX_RESCALE)
00221               optimizedCoefficients_[i] = 1 / (1 + peak);
00222             else
00223               optimizedCoefficients_[i] = 1 / (1 + calib_fac_[i]/weight_sum_[i]);
00224           
00225           } else {
00226             optimizedCoefficients_[i]=1.;
00227             optimizedCoefficientsError_[i] = 0.;
00228             optimizedChiSquare_[i] = -1.;
00229             optimizedIterations_[i] = 0;
00230           }
00231 
00232         }
00233 
00234       else
00235         {
00236           optimizedCoefficientsError_[i] = 0.;
00237           optimizedCoefficients_[i]=1.;
00238           optimizedChiSquare_[i] = -1.;
00239           optimizedIterations_[i] = 0;
00240 //        EcalCalibMap::getMap()->setRingCalib(i, optimizedCoefficients_[i]);
00241 //        //      initialCoefficients_[i] *= optimizedCoefficients_[i];
00242         }
00243       
00244       std::cout << "ZIterativeAlgorithmWithFit::run():Energy Rescaling Coefficient for region " 
00245                 << i << " is "  << optimizedCoefficients_[i] << ", with error "<<optimizedCoefficientsError_[i]<<" - number of events: " << weight_sum_[i] << std::endl;
00246     }
00247   
00248   currentIteration_++;
00249   return kTRUE;
00250 }    
00251 
00252 
00253 bool ZIterativeAlgorithmWithFit::addEvent(calib::CalibElectron* ele1, calib::CalibElectron* ele2, float invMassRescFactor)
00254 {
00255   totalEvents_++;
00256   std::pair<calib::CalibElectron*, calib::CalibElectron*> Electrons(ele1, ele2);
00257 #ifdef DEBUG
00258   std::cout  << "In addEvent " ;
00259   std::cout << ele1->getRecoElectron()->superCluster()->rawEnergy() << " " ;
00260   std::cout << ele1->getRecoElectron()->superCluster()->position().eta() << " " ;
00261   std::cout << ele2->getRecoElectron()->superCluster()->rawEnergy() << " " ;
00262   std::cout << ele2->getRecoElectron()->superCluster()->position().eta() << " " ;
00263   std::cout << std::endl;
00264 #endif
00265 
00266   if (massMethod == "SCTRMass" )
00267     {
00268       massReco_.push_back(invMassCalc(ele1->getRecoElectron()->superCluster()->energy(), ele1->getRecoElectron()->eta(), ele1->getRecoElectron()->phi(), ele2->getRecoElectron()->superCluster()->energy(), ele2->getRecoElectron()->eta(), ele2->getRecoElectron()->phi()));
00269     }
00270 
00271   else if (massMethod == "SCMass" )
00272     {
00273       massReco_.push_back(invMassCalc(ele1->getRecoElectron()->superCluster()->energy(), ele1->getRecoElectron()->superCluster()->position().eta(), ele1->getRecoElectron()->superCluster()->position().phi(), ele2->getRecoElectron()->superCluster()->energy(), ele2->getRecoElectron()->superCluster()->position().eta(), ele2->getRecoElectron()->superCluster()->position().phi()));
00274     }  
00275   
00276 #ifdef DEBUG
00277   std::cout << "Mass calculated " << massReco_[currentEvent_] << std::endl;
00278 #endif
00279 
00280   if((ele2->getRecoElectron()->superCluster()->position().eta() > -10.) && (ele2->getRecoElectron()->superCluster()->position().eta() < 10.) && 
00281      (ele2->getRecoElectron()->superCluster()->position().phi() > -10.) && (ele2->getRecoElectron()->superCluster()->position().phi() < 10.)) {
00282     getWeight(currentEvent_, Electrons, invMassRescFactor);
00283   }
00284 
00285   currentEvent_++;
00286   return kTRUE;
00287 }
00288 
00289     
00290 void ZIterativeAlgorithmWithFit::getWeight(unsigned int event_id, std::pair<calib::CalibElectron*,calib::CalibElectron*> elepair, float invMassRescFactor) 
00291 {
00292 
00293   float event_weight;
00294   if (UseStatWeights_) {
00295     event_weight=getEventWeight(event_id);
00296   } else {
00297     event_weight=1/(elepair.first->getRecoElectron()->superCluster()->rawEnergy()+elepair.second->getRecoElectron()->superCluster()->rawEnergy());
00298   }
00299 
00300   getWeight(event_id, elepair.first, invMassRescFactor);
00301   getWeight(event_id, elepair.second, invMassRescFactor);
00302 }
00303 
00304 float ZIterativeAlgorithmWithFit::getEventWeight(unsigned int event_id) {
00305 
00306   return 1.;
00307 }
00308 
00309 void ZIterativeAlgorithmWithFit::getWeight(unsigned int event_id, calib::CalibElectron* ele, float evweight) {
00310   //  std::cout<< "getting weight for module " << module << " in electron " << ele << std::endl;
00311 
00312   std::vector< std::pair<int,float> > modules=(*ele).getCalibModulesWeights(calibType_); 
00313 
00314   for (int imod=0; imod< (int) modules.size(); imod++) {
00315 
00316     int mod = (int) modules[imod].first;
00317     
00318     if (mod< (int) channels_ && mod>=0) {
00319 
00320       if (modules[imod].second >= 0.12 && modules[imod].second < 10000.) 
00321         {
00322           if( (nCrystalCut_ == -1) || ((!(mod <= nCrystalCut_ - 1 )) &&
00323                                      !((mod > (19-nCrystalCut_)) && (mod <= (19+nCrystalCut_))) &&
00324                                      !((mod > (39-nCrystalCut_)) && (mod <= (39+nCrystalCut_))) &&
00325                                      !((mod > (59-nCrystalCut_)) && (mod <= (59+nCrystalCut_))) &&
00326                                      !((mod > (84-nCrystalCut_)) && (mod <= (84+nCrystalCut_))) &&
00327                                      !((mod > (109-nCrystalCut_)) && (mod <= (109+nCrystalCut_))) &&
00328                                      !((mod > (129-nCrystalCut_)) && (mod <= (129+nCrystalCut_))) &&
00329                                      !((mod > (149-nCrystalCut_)) && (mod <= (149+nCrystalCut_))) &&
00330                                      !(mod > (169-nCrystalCut_))))
00331             {
00332 
00333               float weight2 = modules[imod].second / ele->getRecoElectron()->superCluster()->rawEnergy();
00334 #ifdef DEBUG
00335               std::cout << "w2 " << weight2 << std::endl;
00336 #endif
00337               if (weight2>=0. && weight2<=1.)
00338                 {
00339 
00340                   float rescale = (TMath::Power((massReco_[event_id] / evweight), 2.) - 1) / 2.;
00341 #ifdef DEBUG
00342                   std::cout << "rescale " << rescale << std::endl;                
00343 #endif
00344                   if (rescale>= MIN_RESCALE && rescale<=MAX_RESCALE)
00345                     {
00346                       calib_fac_[mod] += weight2 * rescale;
00347                       weight_sum_[mod]+= weight2;
00348 
00349                       thePlots_->weightedRescaleFactor[currentIteration_][mod]->Fill(rescale,weight2);
00350                       thePlots_->unweightedRescaleFactor[currentIteration_][mod]->Fill(rescale,1.);
00351                       thePlots_->weight[currentIteration_][mod]->Fill(weight2,1.);
00352                     }
00353                   else
00354                     {
00355                       std::cout     << "[ZIterativeAlgorithmWithFit]::[getWeight]::rescale out " << rescale << std::endl;
00356                     }
00357                 }
00358 
00359             }
00360         }
00361     } 
00362     else 
00363       {
00364         std::cout << "ZIterativeAlgorithmWithFit::FATAL:found a wrong module_id " << mod << " channels " << channels_ << std::endl;
00365       }
00366   }
00367 }
00368 
00369 
00370 ZIterativeAlgorithmWithFit::~ZIterativeAlgorithmWithFit()
00371 {
00372 
00373 }
00374 
00375 void ZIterativeAlgorithmWithFit::gausfit(TH1F * histoou,double* par,double* errpar,float nsigmalow, float nsigmaup, double* myChi2, int* iterations) {
00376   TF1 *gausa = new TF1("gausa","gaus",histoou->GetMean()-3*histoou->GetRMS(),histoou->GetMean()+3*histoou->GetRMS());
00377   
00378   gausa->SetParameters(histoou->GetMaximum(),histoou->GetMean(),histoou->GetRMS());
00379   
00380   histoou->Fit("gausa","qR0N");
00381   
00382   double p1    = gausa->GetParameter(1);
00383   double sigma = gausa->GetParameter(2);
00384   double nor   = gausa->GetParameter(0);
00385     
00386   double xmi=p1-5*sigma;
00387   double xma=p1+5*sigma;
00388   double chi2=100;
00389   
00390   double xmin_fit=p1-nsigmalow*sigma;
00391   double xmax_fit=p1+nsigmaup*sigma;
00392   
00393   int iter=0;
00394   TF1* fitFunc;
00395 
00396   while ((chi2>1. && iter<5) || iter<2 )
00397     {
00398       xmin_fit=p1-nsigmalow*sigma;
00399       xmax_fit=p1+nsigmaup*sigma;
00400       xmi=p1-5*sigma;
00401       xma=p1+5*sigma;
00402       
00403       char suffix[20];
00404       sprintf (suffix,"_iter_%d",iter); 
00405       fitFunc = new TF1("FitFunc"+TString(suffix),"gaus",xmin_fit,xmax_fit);
00406       fitFunc->SetParameters(nor,p1,sigma);
00407       fitFunc->SetLineColor((int)(iter+1));
00408       fitFunc->SetLineWidth(1.);
00409       //histoou->Fit("FitFunc","lR+","");
00410       histoou->Fit("FitFunc"+TString(suffix),"qR0+","");
00411       
00412       histoou->GetXaxis()->SetRangeUser(xmi,xma);
00413       histoou->GetXaxis()->SetLabelSize(0.055);
00414       
00415       //      cout << fitFunc->GetParameters() << "," << par << endl;
00416       par[0]=(fitFunc->GetParameters())[0];
00417       par[1]=(fitFunc->GetParameters())[1];
00418       par[2]=(fitFunc->GetParameters())[2];
00419       errpar[0]=(fitFunc->GetParErrors())[0];
00420       errpar[1]=(fitFunc->GetParErrors())[1];
00421       errpar[2]=(fitFunc->GetParErrors())[2];
00422       if (fitFunc->GetNDF()!=0)
00423         {
00424           chi2=fitFunc->GetChisquare()/(fitFunc->GetNDF());
00425           *myChi2 = chi2;
00426  
00427         }
00428       else
00429         {
00430           chi2=100.;
00431 //           par[0]=-99;
00432 //           par[1]=-99;
00433 //           par[2]=-99;
00434           std::cout << "WARNING: Not enough NDF" << std::endl;
00435 //           return 0;
00436         }
00437 
00438       // Non visualizzare
00439       //      histoou->Draw();
00440       //      c1->Update();
00441 
00442       //      cout << "iter " << iter << " chi2 " << chi2 << endl;
00443       nor=par[0];
00444       p1=par[1];
00445       sigma=par[2];
00446       iter++;
00447       *iterations = iter;
00448 
00449     }
00450   return;
00451 }
00452 
00453