Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:30 

0001 /* ******************************************
0002  * ZIterativeAlgorithmWithFit.cc
0003  *
0004  * 
0005  * Paolo Meridiani 06/07/2005
0006  * Rewritten for CMSSW 04/06/2007
0007  ********************************************/
0008 
0009 #include "Calibration/Tools/interface/ZIterativeAlgorithmWithFit.h"
0010 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
0011 #include "Calibration/Tools/interface/EcalIndexingTools.h"
0012 
0013 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0014 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016 
0017 #include <TMath.h>
0018 #include <TCanvas.h>
0019 #include "TH1.h"
0020 #include "TH2.h"
0021 #include "TF1.h"
0022 #include "TH1F.h"
0023 #include "TMinuit.h"
0024 #include "TGraphErrors.h"
0025 #include "THStack.h"
0026 #include "TLegend.h"
0027 
0028 #include <fstream>
0029 #include <iostream>
0030 #include <vector>
0031 
0032 //#include "Tools.C"
0033 
0034 //Scale and Bins for calibration factor histograms
0035 #define MIN_RESCALE -0.5
0036 #define MAX_RESCALE 0.5
0037 #define NBINS_LOWETA 100
0038 #define NBINS_HIGHETA 50
0039 
0040 const double ZIterativeAlgorithmWithFit::M_Z_ = 91.187;
0041 
0042 //  #if !defined(__CINT__)
0043 //  ClassImp(Electron)
0044 //  #endif
0045 
0046 ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFit() {
0047   // std::cout<< "ZIterativeAlgorithmWithFit::Called Construnctor" << std::endl;
0048   numberOfIterations_ = 10;
0049   channels_ = 1;
0050   totalEvents_ = 0;
0051   currentEvent_ = 0;
0052   currentIteration_ = 0;
0053   optimizedCoefficients_.resize(channels_);
0054   optimizedCoefficientsError_.resize(channels_);
0055   optimizedChiSquare_.resize(channels_);
0056   optimizedIterations_.resize(channels_);
0057   calib_fac_.resize(channels_);
0058   weight_sum_.resize(channels_);
0059   electrons_.resize(1);
0060   massReco_.resize(1);
0061 }
0062 
0063 ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFit(const edm::ParameterSet& ps)
0064 //k, unsigned int events)
0065 {
0066   //std::cout<< "ZIterativeAlgorithmWithFit::Called Constructor" << std::endl;
0067   numberOfIterations_ = ps.getUntrackedParameter<unsigned int>("maxLoops", 0);
0068   massMethod = ps.getUntrackedParameter<std::string>("ZCalib_InvMass", "SCMass");
0069   calibType_ = ps.getUntrackedParameter<std::string>("ZCalib_CalibType", "RING");
0070 
0071   if (calibType_ == "RING")
0072     channels_ = EcalRingCalibrationTools::N_RING_TOTAL;
0073   else if (calibType_ == "MODULE")
0074 
0075     channels_ = EcalRingCalibrationTools::N_MODULES_BARREL;
0076   else if (calibType_ == "ABS_SCALE")
0077     channels_ = 1;
0078   else if (calibType_ == "ETA_ET_MODE")
0079     channels_ = EcalIndexingTools::getInstance()->getNumberOfChannels();
0080 
0081   std::cout << "[ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFit] channels_ = " << channels_ << std::endl;
0082 
0083   nCrystalCut_ = ps.getUntrackedParameter<int>("ZCalib_nCrystalCut", -1);
0084 
0085   //Resetting currentEvent & iteration
0086   currentEvent_ = 0;
0087   currentIteration_ = 0;
0088   totalEvents_ = 0;
0089 
0090   optimizedCoefficients_.resize(channels_);
0091   optimizedCoefficientsError_.resize(channels_);
0092   optimizedChiSquare_.resize(channels_);
0093   optimizedIterations_.resize(channels_);
0094   calib_fac_.resize(channels_);
0095   weight_sum_.resize(channels_);
0096 
0097   //Creating and booking histograms
0098   thePlots_ = new ZIterativeAlgorithmWithFitPlots;
0099   bookHistograms();
0100 
0101   //Setting up rescaling if needed
0102   UseStatWeights_ = ps.getUntrackedParameter<bool>("ZCalib_UseStatWeights", false);
0103   if (UseStatWeights_) {
0104     WeightFileName_ = "weights.txt";
0105     StatWeights_.resize(channels_);
0106     getStatWeights(WeightFileName_);
0107     //    Event_Weight_.resize(events);
0108   }
0109 }
0110 
0111 void ZIterativeAlgorithmWithFit::bookHistograms() {
0112   if (!thePlots_)
0113     return;
0114 
0115   for (unsigned int i2 = 0; i2 < numberOfIterations_; i2++) {
0116     for (unsigned int i1 = 0; i1 < channels_; i1++) {
0117       char histoName[200];
0118       char histoTitle[200];
0119 
0120       //WeightedRescaling factor
0121       sprintf(histoName, "WeightedRescaleFactor_channel_%d_Iteration_%d", i1, i2);
0122       sprintf(histoTitle, "WeightedRescaleFactor Channel_%d Iteration %d", i1, i2);
0123       if (i1 > 15 && i1 < 155)
0124         thePlots_->weightedRescaleFactor[i2][i1] =
0125             new TH1F(histoName, histoTitle, NBINS_LOWETA, MIN_RESCALE, MAX_RESCALE);
0126       else
0127         thePlots_->weightedRescaleFactor[i2][i1] =
0128             new TH1F(histoName, histoTitle, NBINS_HIGHETA, MIN_RESCALE, MAX_RESCALE);
0129       thePlots_->weightedRescaleFactor[i2][i1]->GetXaxis()->SetTitle("Rescale factor");
0130       thePlots_->weightedRescaleFactor[i2][i1]->GetYaxis()->SetTitle("a.u.");
0131 
0132       //UnweightedRescaling factor
0133       sprintf(histoName, "UnweightedRescaleFactor_channel_%d_Iteration_%d", i1, i2);
0134       sprintf(histoTitle, "UnweightedRescaleFactor Channel_%d Iteration %d", i1, i2);
0135       if (i1 > 15 && i1 < 155)
0136         thePlots_->unweightedRescaleFactor[i2][i1] =
0137             new TH1F(histoName, histoTitle, NBINS_LOWETA, MIN_RESCALE, MAX_RESCALE);
0138       else
0139         thePlots_->unweightedRescaleFactor[i2][i1] =
0140             new TH1F(histoName, histoTitle, NBINS_HIGHETA, MIN_RESCALE, MAX_RESCALE);
0141       thePlots_->unweightedRescaleFactor[i2][i1]->GetXaxis()->SetTitle("Rescale factor");
0142       thePlots_->unweightedRescaleFactor[i2][i1]->GetYaxis()->SetTitle("a.u.");
0143 
0144       //Weights
0145       sprintf(histoName, "Weight_channel_%d_Iteration_%d", i1, i2);
0146       sprintf(histoTitle, "Weight Channel_%d Iteration %d", i1, i2);
0147       thePlots_->weight[i2][i1] = new TH1F(histoName, histoTitle, 100, 0., 1.);
0148       thePlots_->weight[i2][i1]->GetXaxis()->SetTitle("Weight");
0149       thePlots_->weight[i2][i1]->GetYaxis()->SetTitle("a.u.");
0150     }
0151   }
0152 }
0153 
0154 void ZIterativeAlgorithmWithFit::getStatWeights(const std::string& file) {
0155   std::ifstream statfile;
0156   statfile.open(file.c_str());
0157   if (!statfile) {
0158     std::cout << "ZIterativeAlgorithmWithFit::FATAL: stat weight  file " << file << " not found" << std::endl;
0159     exit(-1);
0160   }
0161   for (unsigned int i = 0; i < channels_; i++) {
0162     int imod;
0163     statfile >> imod >> StatWeights_[i];
0164     //std::cout << "Read Stat Weight for module " << imod << ": " <<  StatWeights_[i] << std::endl;
0165   }
0166 }
0167 
0168 bool ZIterativeAlgorithmWithFit::resetIteration() {
0169   totalEvents_ = 0;
0170   currentEvent_ = 0;
0171 
0172   //Reset correction
0173   massReco_.clear();
0174   for (unsigned int i = 0; i < channels_; i++)
0175     calib_fac_[i] = 0.;
0176   for (unsigned int i = 0; i < channels_; i++)
0177     weight_sum_[i] = 0.;
0178 
0179   return kTRUE;
0180 }
0181 
0182 bool ZIterativeAlgorithmWithFit::iterate() {
0183   //Found optimized coefficients
0184   for (int i = 0; i < (int)channels_; i++) {
0185     //RP      if (weight_sum_[i]!=0. && calib_fac_[i]!=0.) {
0186     if ((nCrystalCut_ == -1) ||
0187         ((!(i <= nCrystalCut_ - 1)) && !((i > (19 - nCrystalCut_)) && (i <= (19 + nCrystalCut_))) &&
0188          !((i > (39 - nCrystalCut_)) && (i <= (39 + nCrystalCut_))) &&
0189          !((i > (59 - nCrystalCut_)) && (i <= (59 + nCrystalCut_))) &&
0190          !((i > (84 - nCrystalCut_)) && (i <= (84 + nCrystalCut_))) &&
0191          !((i > (109 - nCrystalCut_)) && (i <= (109 + nCrystalCut_))) &&
0192          !((i > (129 - nCrystalCut_)) && (i <= (129 + nCrystalCut_))) &&
0193          !((i > (149 - nCrystalCut_)) && (i <= (149 + nCrystalCut_))) && !(i > (169 - nCrystalCut_)))) {
0194       if (weight_sum_[i] != 0.) {
0195         //optimizedCoefficients_[i] = calib_fac_[i]/weight_sum_[i];
0196 
0197         double gausFitParameters[3], gausFitParameterErrors[3], gausFitChi2;
0198         int gausFitIterations;
0199 
0200         gausfit((TH1F*)thePlots_->weightedRescaleFactor[currentIteration_][i],
0201                 gausFitParameters,
0202                 gausFitParameterErrors,
0203                 2.5,
0204                 2.5,
0205                 &gausFitChi2,
0206                 &gausFitIterations);
0207 
0208         float peak = gausFitParameters[1];
0209         float peakError = gausFitParameterErrors[1];
0210         float chi2 = gausFitChi2;
0211 
0212         int iters = gausFitIterations;
0213 
0214         optimizedCoefficientsError_[i] = peakError;
0215         optimizedChiSquare_[i] = chi2;
0216         optimizedIterations_[i] = iters;
0217 
0218         if (peak >= MIN_RESCALE && peak <= MAX_RESCALE)
0219           optimizedCoefficients_[i] = 1 / (1 + peak);
0220         else
0221           optimizedCoefficients_[i] = 1 / (1 + calib_fac_[i] / weight_sum_[i]);
0222 
0223       } else {
0224         optimizedCoefficients_[i] = 1.;
0225         optimizedCoefficientsError_[i] = 0.;
0226         optimizedChiSquare_[i] = -1.;
0227         optimizedIterations_[i] = 0;
0228       }
0229 
0230     }
0231 
0232     else {
0233       optimizedCoefficientsError_[i] = 0.;
0234       optimizedCoefficients_[i] = 1.;
0235       optimizedChiSquare_[i] = -1.;
0236       optimizedIterations_[i] = 0;
0237       //      EcalCalibMap::getMap()->setRingCalib(i, optimizedCoefficients_[i]);
0238       //      //      initialCoefficients_[i] *= optimizedCoefficients_[i];
0239     }
0240 
0241     std::cout << "ZIterativeAlgorithmWithFit::run():Energy Rescaling Coefficient for region " << i << " is "
0242               << optimizedCoefficients_[i] << ", with error " << optimizedCoefficientsError_[i]
0243               << " - number of events: " << weight_sum_[i] << std::endl;
0244   }
0245 
0246   currentIteration_++;
0247   return kTRUE;
0248 }
0249 
0250 bool ZIterativeAlgorithmWithFit::addEvent(calib::CalibElectron* ele1,
0251                                           calib::CalibElectron* ele2,
0252                                           float invMassRescFactor) {
0253   totalEvents_++;
0254   std::pair<calib::CalibElectron*, calib::CalibElectron*> Electrons(ele1, ele2);
0255 #ifdef DEBUG
0256   std::cout << "In addEvent ";
0257   std::cout << ele1->getRecoElectron()->superCluster()->rawEnergy() << " ";
0258   std::cout << ele1->getRecoElectron()->superCluster()->position().eta() << " ";
0259   std::cout << ele2->getRecoElectron()->superCluster()->rawEnergy() << " ";
0260   std::cout << ele2->getRecoElectron()->superCluster()->position().eta() << " ";
0261   std::cout << std::endl;
0262 #endif
0263 
0264   if (massMethod == "SCTRMass") {
0265     massReco_.push_back(invMassCalc(ele1->getRecoElectron()->superCluster()->energy(),
0266                                     ele1->getRecoElectron()->eta(),
0267                                     ele1->getRecoElectron()->phi(),
0268                                     ele2->getRecoElectron()->superCluster()->energy(),
0269                                     ele2->getRecoElectron()->eta(),
0270                                     ele2->getRecoElectron()->phi()));
0271   }
0272 
0273   else if (massMethod == "SCMass") {
0274     massReco_.push_back(invMassCalc(ele1->getRecoElectron()->superCluster()->energy(),
0275                                     ele1->getRecoElectron()->superCluster()->position().eta(),
0276                                     ele1->getRecoElectron()->superCluster()->position().phi(),
0277                                     ele2->getRecoElectron()->superCluster()->energy(),
0278                                     ele2->getRecoElectron()->superCluster()->position().eta(),
0279                                     ele2->getRecoElectron()->superCluster()->position().phi()));
0280   }
0281 
0282 #ifdef DEBUG
0283   std::cout << "Mass calculated " << massReco_[currentEvent_] << std::endl;
0284 #endif
0285 
0286   if ((ele2->getRecoElectron()->superCluster()->position().eta() > -10.) &&
0287       (ele2->getRecoElectron()->superCluster()->position().eta() < 10.) &&
0288       (ele2->getRecoElectron()->superCluster()->position().phi() > -10.) &&
0289       (ele2->getRecoElectron()->superCluster()->position().phi() < 10.)) {
0290     getWeight(currentEvent_, Electrons, invMassRescFactor);
0291   }
0292 
0293   currentEvent_++;
0294   return kTRUE;
0295 }
0296 
0297 void ZIterativeAlgorithmWithFit::getWeight(unsigned int event_id,
0298                                            std::pair<calib::CalibElectron*, calib::CalibElectron*> elepair,
0299                                            float invMassRescFactor) {
0300   getWeight(event_id, elepair.first, invMassRescFactor);
0301   getWeight(event_id, elepair.second, invMassRescFactor);
0302 }
0303 
0304 float ZIterativeAlgorithmWithFit::getEventWeight(unsigned int event_id) { return 1.; }
0305 
0306 void ZIterativeAlgorithmWithFit::getWeight(unsigned int event_id, calib::CalibElectron* ele, float evweight) {
0307   //  std::cout<< "getting weight for module " << module << " in electron " << ele << std::endl;
0308 
0309   std::vector<std::pair<int, float> > modules = (*ele).getCalibModulesWeights(calibType_);
0310 
0311   for (int imod = 0; imod < (int)modules.size(); imod++) {
0312     int mod = (int)modules[imod].first;
0313 
0314     if (mod < (int)channels_ && mod >= 0) {
0315       if (modules[imod].second >= 0.12 && modules[imod].second < 10000.) {
0316         if ((nCrystalCut_ == -1) ||
0317             ((!(mod <= nCrystalCut_ - 1)) && !((mod > (19 - nCrystalCut_)) && (mod <= (19 + nCrystalCut_))) &&
0318              !((mod > (39 - nCrystalCut_)) && (mod <= (39 + nCrystalCut_))) &&
0319              !((mod > (59 - nCrystalCut_)) && (mod <= (59 + nCrystalCut_))) &&
0320              !((mod > (84 - nCrystalCut_)) && (mod <= (84 + nCrystalCut_))) &&
0321              !((mod > (109 - nCrystalCut_)) && (mod <= (109 + nCrystalCut_))) &&
0322              !((mod > (129 - nCrystalCut_)) && (mod <= (129 + nCrystalCut_))) &&
0323              !((mod > (149 - nCrystalCut_)) && (mod <= (149 + nCrystalCut_))) && !(mod > (169 - nCrystalCut_)))) {
0324           float weight2 = modules[imod].second / ele->getRecoElectron()->superCluster()->rawEnergy();
0325 #ifdef DEBUG
0326           std::cout << "w2 " << weight2 << std::endl;
0327 #endif
0328           if (weight2 >= 0. && weight2 <= 1.) {
0329             float rescale = (TMath::Power((massReco_[event_id] / evweight), 2.) - 1) / 2.;
0330 #ifdef DEBUG
0331             std::cout << "rescale " << rescale << std::endl;
0332 #endif
0333             if (rescale >= MIN_RESCALE && rescale <= MAX_RESCALE) {
0334               calib_fac_[mod] += weight2 * rescale;
0335               weight_sum_[mod] += weight2;
0336 
0337               thePlots_->weightedRescaleFactor[currentIteration_][mod]->Fill(rescale, weight2);
0338               thePlots_->unweightedRescaleFactor[currentIteration_][mod]->Fill(rescale, 1.);
0339               thePlots_->weight[currentIteration_][mod]->Fill(weight2, 1.);
0340             } else {
0341               std::cout << "[ZIterativeAlgorithmWithFit]::[getWeight]::rescale out " << rescale << std::endl;
0342             }
0343           }
0344         }
0345       }
0346     } else {
0347       std::cout << "ZIterativeAlgorithmWithFit::FATAL:found a wrong module_id " << mod << " channels " << channels_
0348                 << std::endl;
0349     }
0350   }
0351 }
0352 
0353 ZIterativeAlgorithmWithFit::~ZIterativeAlgorithmWithFit() {}
0354 
0355 void ZIterativeAlgorithmWithFit::gausfit(
0356     TH1F* histoou, double* par, double* errpar, float nsigmalow, float nsigmaup, double* myChi2, int* iterations) {
0357   auto gausa = std::make_unique<TF1>(
0358       "gausa", "gaus", histoou->GetMean() - 3 * histoou->GetRMS(), histoou->GetMean() + 3 * histoou->GetRMS());
0359 
0360   gausa->SetParameters(histoou->GetMaximum(), histoou->GetMean(), histoou->GetRMS());
0361 
0362   histoou->Fit(gausa.get(), "qR0N");
0363 
0364   double p1 = gausa->GetParameter(1);
0365   double sigma = gausa->GetParameter(2);
0366   double nor = gausa->GetParameter(0);
0367 
0368   double xmi, xma, xmin_fit, xmax_fit;
0369   double chi2 = 100;
0370   int iter = 0;
0371 
0372   while ((chi2 > 1. && iter < 5) || iter < 2) {
0373     xmin_fit = p1 - nsigmalow * sigma;
0374     xmax_fit = p1 + nsigmaup * sigma;
0375     xmi = p1 - 5 * sigma;
0376     xma = p1 + 5 * sigma;
0377 
0378     char suffix[20];
0379     sprintf(suffix, "_iter_%d", iter);
0380     auto fitFunc = std::make_unique<TF1>("FitFunc" + TString(suffix), "gaus", xmin_fit, xmax_fit);
0381     fitFunc->SetParameters(nor, p1, sigma);
0382     fitFunc->SetLineColor((int)(iter + 1));
0383     fitFunc->SetLineWidth(1);
0384     //histoou->Fit("FitFunc","lR+","");
0385     histoou->Fit(fitFunc.get(), "qR0+", "");
0386 
0387     histoou->GetXaxis()->SetRangeUser(xmi, xma);
0388     histoou->GetXaxis()->SetLabelSize(0.055);
0389 
0390     //      std::cout << fitFunc->GetParameters() << "," << par << std::endl;
0391     par[0] = (fitFunc->GetParameters())[0];
0392     par[1] = (fitFunc->GetParameters())[1];
0393     par[2] = (fitFunc->GetParameters())[2];
0394     errpar[0] = (fitFunc->GetParErrors())[0];
0395     errpar[1] = (fitFunc->GetParErrors())[1];
0396     errpar[2] = (fitFunc->GetParErrors())[2];
0397     if (fitFunc->GetNDF() != 0) {
0398       chi2 = fitFunc->GetChisquare() / (fitFunc->GetNDF());
0399       *myChi2 = chi2;
0400 
0401     } else {
0402       chi2 = 100.;
0403       //           par[0]=-99;
0404       //           par[1]=-99;
0405       //           par[2]=-99;
0406       std::cout << "WARNING: Not enough NDF" << std::endl;
0407       //           return 0;
0408     }
0409 
0410     // Non visualizzare
0411     //      histoou->Draw();
0412     //      c1->Update();
0413 
0414     //      std::cout << "iter " << iter << " chi2 " << chi2 << std::endl;
0415     nor = par[0];
0416     p1 = par[1];
0417     sigma = par[2];
0418     iter++;
0419     *iterations = iter;
0420   }
0421   return;
0422 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!