Reference Manual

Go to the documentation of this file.
00001 
00002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00003 
00004 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00006 #include "DataFormats/DetId/interface/DetId.h"
00007 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
00008 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
00009 #include "Calibration/Tools/interface/calibXMLwriter.h"
00010 #include "Calibration/Tools/interface/CalibrationCluster.h"
00011 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
00012 #include "Calibration/Tools/interface/MinL3Algorithm.h"
00013 #include "Calibration/EcalCalibAlgos/interface/ElectronCalibration.h"
00014 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00015 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
00016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00017 #include "DataFormats/TrackReco/interface/Track.h"
00018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00019 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
00020 #include "FWCore/Utilities/interface/isFinite.h"
00021 #include "TFile.h"
00022 #include "TH1.h"
00023 #include "TH2.h"
00024 #include "TF1.h"
00025 #include "TRandom.h"
00026 
00027 #include <iostream>
00028 #include <string>
00029 #include <stdexcept>
00030 #include <vector>
00031 
00032 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00033 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00034 
00035 ElectronCalibration::ElectronCalibration(const edm::ParameterSet& iConfig)
00036 {
00037 
00038    rootfile_ = iConfig.getParameter<std::string>("rootfile");
00039    recHitLabel_ = iConfig.getParameter< edm::InputTag > ("ebRecHitsLabel");
00040    electronLabel_ = iConfig.getParameter< edm::InputTag > ("electronLabel");
00041    trackLabel_ = iConfig.getParameter< edm::InputTag > ("trackLabel");
00042    calibAlgo_       = iConfig.getParameter<std::string>("CALIBRATION_ALGO");
00043    std::cout << " The used Algorithm is  " << calibAlgo_ << std::endl;
00044    keventweight_ = iConfig.getParameter<int>("keventweight");
00045    ClusterSize_ = iConfig.getParameter<int>("Clustersize");
00046    ElePt_ = iConfig.getParameter<double>("ElePt");
00047    maxeta_ = iConfig.getParameter<int>("maxeta");
00048    mineta_ = iConfig.getParameter<int>("mineta");
00049    maxphi_ = iConfig.getParameter<int>("maxphi");
00050    minphi_ = iConfig.getParameter<int>("minphi");
00051    cut1_ = iConfig.getParameter<double>("cut1");
00052    cut2_ = iConfig.getParameter<double>("cut2");
00053    cut3_ = iConfig.getParameter<double>("cut3");
00054    elecclass_ = iConfig.getParameter<int>("elecclass");
00055    std::cout << " The electronclass is " << elecclass_ <<std::endl;
00056    numevent_ = iConfig.getParameter<int>("numevent");
00057    miscalibfile_ = iConfig.getParameter<std::string>("miscalibfile");
00058 
00059    cutEPCalo1_ = iConfig.getParameter<double>("cutEPCaloMin");
00060    cutEPCalo2_ = iConfig.getParameter<double>("cutEPCaloMax");
00061    cutEPin1_ = iConfig.getParameter<double>("cutEPinMin");
00062    cutEPin2_ = iConfig.getParameter<double>("cutEPinMax");
00063    cutCalo1_ = iConfig.getParameter<double>("cutCaloMin");
00064    cutCalo2_ = iConfig.getParameter<double>("cutCaloMax");
00065 
00066    cutESeed_ = iConfig.getParameter<double>("cutESeed");
00067 }
00068 
00069 
00070 ElectronCalibration::~ElectronCalibration()
00071 {
00072   
00073   
00074 }
00075 
00076 //========================================================================
00077 void ElectronCalibration::beginJob() {
00078   //========================================================================
00079    f = new TFile(rootfile_.c_str(),"RECREATE");
00080 
00081   // Book histograms 
00082    e9 = new TH1F("e9","E9 energy", 300, 0., 150.);
00083   e25 = new TH1F("e25","E25 energy", 300, 0., 150.);
00084   scE = new TH1F("scE","SC energy", 300, 0., 150.);
00085   trP = new TH1F("trP","Trk momentum", 300, 0., 150.);
00086   EoP = new TH1F("EoP","EoP", 600, 0., 3.);
00087   EoP_all = new TH1F("EoP_all","EoP_all",600, 0., 3.);
00088 
00089   if (elecclass_ ==0 || elecclass_ == -1) { 
00090     calibs = new TH1F("calib","Calibration constants", 4000, 0.5, 2.);
00091   }else{
00092     calibs = new TH1F("calib","Calibration constants", 800, 0.5, 2.);
00093   }
00094 
00095   e25OverScE = new TH1F("e25OverscE","E25 / SC energy", 400, 0., 2.);
00096   E25oP = new TH1F("E25oP","E25 / P", 1200, 0., 1.5);
00097 
00098   Map = new TH2F("Map","Nb Events in Crystal",85,1, 85,70 ,5, 75);
00099   e9Overe25 = new TH1F("e9Overe25","E9 / E25", 400, 0., 2.);
00100   Map3Dcalib = new TH2F("3Dcalib", "3Dcalib",85 ,1 ,85,70, 5, 75 );
00101 
00102   MapCor1 = new TH2F ("MapCor1", "Correlation E25/Pcalo versus E25/Pin",100 ,0. ,5. ,100,0.,5. );
00103   MapCor2 = new TH2F ("MapCor2", "Correlation E25/Pcalo versus E/P",100 ,0. ,5. ,100,0.,5. );
00104   MapCor3 = new TH2F ("MapCor3", "Correlation E25/Pcalo versus Pout/Pin",100 ,0. ,5. ,100,0.,5. );
00105   MapCor4 = new TH2F ("MapCor4", "Correlation E25/Pcalo versus E25/highestP",100 ,0. ,5. ,100,0.,5. );
00106   MapCor5 = new TH2F ("MapCor5", "Correlation E25/Pcalo versus Pcalo/Pout",100 ,0. ,5. ,100,0.,5. );
00107   MapCor6 = new TH2F ("MapCor6", "Correlation Pout/Pin versus E25/Pin",100 ,0. ,5. ,100,0.,5. );
00108   MapCor7 = new TH2F ("MapCor7", "Correlation Pout/Pin versus Pcalo/Pout",100 ,0. ,5. ,100,0.,5. );
00109   MapCor8 = new TH2F ("MapCor8", "Correlation E25/Pin versus Pcalo/Pout",100 ,0. ,5. ,100,0.,5. );
00110   MapCor9 = new TH2F ("MapCor9", "Correlation  E25/Pcalo versus Eseed/Pout",100 ,0. ,5. ,100,0.,5. );
00111   MapCor10 = new TH2F ("MapCor10", "Correlation Eseed/Pout versus Pout/Pin",100 ,0. ,5. ,100,0.,5. );
00112   MapCor11 = new TH2F ("MapCor11", "Correlation Eseed/Pout versus E25/Pin",100 ,0. ,5. ,100,0.,5. );
00113   MapCorCalib = new TH2F ("MapCorCalib", "Correlation Miscalibration versus Calibration constants", 100, 0.5,1.5, 100, 0.5, 1.5);
00114 
00115   PinMinPout = new TH1F("PinMinPout","(Pin - Pout)/Pin",600,-2.0,2.0);
00116 
00117   if(elecclass_ == 0 || elecclass_ == -1) { 
00118     calibinter = new TH1F("calibinter", "internal calibration constants", 2000 , 0.5,2.);
00119     PinOverPout= new TH1F("PinOverPout", "pinOverpout", 600,0., 3.);
00120     eSeedOverPout= new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
00121     MisCalibs = new TH1F("MisCalibs","Miscalibration constants",4000,0.5,2.);
00122     RatioCalibs = new TH1F("RatioCalibs","Ratio in Calibration Constants", 4000, 0.5, 2.0);
00123     DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 4000, -1.0,1.0);
00124   }else { 
00125     calibinter = new TH1F("calibinter", "internal calibration constants",400 , 0.5,2.);
00126     PinOverPout= new TH1F("PinOverPout", "pinOverpout", 600,0., 3.);
00127     eSeedOverPout= new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
00128     MisCalibs = new TH1F("MisCalibs","Miscalibration constants",800,0.5,2.);
00129     RatioCalibs = new TH1F("RatioCalibs","Ratio in Calibration Constants", 800, 0.5, 2.0);
00130     DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 800, -1.0,1.0);
00131   }
00132   Error1 = new TH1F ("Error1","DeltaP/Pin",800 ,-1.0,1.0 );
00133   Error2 = new TH1F ("Error2","DeltaP/Pout",800 ,-1.0,1.0 );
00134   Error3 = new TH1F ("Error3","DeltaP/Pcalo",800 ,-1.0,1.0 );
00135   eSeedOverPout2= new TH1F("eSeedOverPout2", "eSeedOverpout (No Supercluster)", 600, 0., 4.);
00136   hadOverEm= new TH1F("hadOverEm", "Had/EM distribution", 600, -2., 2.);
00137   
00138   // Book histograms  
00139   Map3DcalibNoCuts = new TH2F("3DcalibNoCuts", "3Dcalib (Before Cuts)",85 ,1 ,85,70, 5, 75 );
00140   e9NoCuts = new TH1F("e9NoCuts","E9 energy (Before Cuts)",300, 0., 150.);
00141   e25NoCuts = new TH1F("e25NoCuts","E25 energy (Before Cuts)", 300, 0., 150.);
00142   scENoCuts = new TH1F("scENoCuts","SC energy (Before Cuts)", 300, 0., 150.);
00143   trPNoCuts = new TH1F("trPNoCuts","Trk momentum (Before Cuts)", 300, 0., 150.);
00144   EoPNoCuts = new TH1F("EoPNoCuts","EoP (Before Cuts)", 600, 0., 3.);
00145   if (elecclass_ ==0 || elecclass_ == -1){ 
00146     calibsNoCuts = new TH1F("calibNoCuts","Calibration constants (Before Cuts)", 4000, 0., 2.);
00147   }else{
00148     calibsNoCuts = new TH1F("calibNoCuts","Calibration constants (Before Cuts)", 800, 0., 2.);
00149   }
00150   e25OverScENoCuts = new TH1F("e25OverscENoCuts","E25 / SC energy (Before Cuts)", 400, 0., 2.);
00151   E25oPNoCuts = new TH1F("E25oPNoCuts","E25 / P (Before Cuts)", 1200, 0., 1.5);
00152   MapNoCuts = new TH2F("MapNoCuts","Nb Events in Crystal (Before Cuts)",85,1, 85,70 ,5, 75);
00153   e9Overe25NoCuts = new TH1F("e9Overe25NoCuts","E9 / E25 (Before Cuts)", 400, 0., 2.);
00154   PinOverPoutNoCuts = new TH1F("PinOverPoutNoCuts", "pinOverpout (Before Cuts)", 600,0., 3.);
00155   eSeedOverPoutNoCuts = new TH1F(" eSeedOverPoutNoCuts", "eSeedOverpout (Before Cuts) ", 600, 0., 4.);
00156   PinMinPoutNoCuts = new TH1F("PinMinPoutNoCuts","(Pin - Pout)/Pin (Before Cuts)",600,-2.0,2.0);
00157 
00158   RatioCalibsNoCuts = new TH1F("RatioCalibsNoCuts","Ratio in Calibration Constants (Before Cuts)", 4000, 0.5, 2.0);
00159   DiffCalibsNoCuts = new TH1F("DiffCalibsNoCuts", "Difference in Calibration constants (Before Cuts)", 4000, -1.0,1.0);
00160   calibinterNoCuts = new TH1F("calibinterNoCuts", "internal calibration constants", 2000 , 0.5,2.);
00161  
00162   MapCor1NoCuts = new TH2F ("MapCor1NoCuts", "Correlation E25/PatCalo versus E25/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00163   MapCor2NoCuts = new TH2F ("MapCor2NoCuts", "Correlation E25/PatCalo versus E/P (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00164   MapCor3NoCuts = new TH2F ("MapCor3NoCuts", "Correlation E25/PatCalo versus Pout/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00165   MapCor4NoCuts = new TH2F ("MapCor4NoCuts", "Correlation E25/PatCalo versus E25/highestP (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00166   MapCor5NoCuts = new TH2F ("MapCor5NoCuts", "Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00167   MapCor6NoCuts = new TH2F ("MapCor6NoCuts", "Correlation Pout/Pin versus E25/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00168   MapCor7NoCuts = new TH2F ("MapCor7NoCuts", "Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00169   MapCor8NoCuts = new TH2F ("MapCor8NoCuts", "Correlation E25/Pin versus Pcalo/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00170   MapCor9NoCuts = new TH2F ("MapCor9NoCuts", "Correlation  E25/Pcalo versus Eseed/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00171   MapCor10NoCuts = new TH2F ("MapCor10NoCuts", "Correlation Eseed/Pout versus Pout/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00172   MapCor11NoCuts = new TH2F ("MapCor11NoCuts", "Correlation Eseed/Pout versus E25/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00173   MapCorCalibNoCuts = new TH2F ("MapCorCalibNoCuts", "Correlation Miscalibration versus Calibration constants (Before Cuts)", 100, 0., 3., 100, 0., 3.);
00174 
00175   Error1NoCuts = new TH1F ("Eror1NoCuts","DeltaP/Pin (Before Cuts)",800 ,-1.0,1.0 );
00176   Error2NoCuts = new TH1F ("Error2NoCuts","DeltaP/Pout (Before Cuts)",800 ,-1.0,1.0 );
00177   Error3NoCuts = new TH1F ("Error3NoCuts","DeltaP/Pcalo (Before Cuts)",800 ,-1.0, 1.0);
00178   eSeedOverPout2NoCuts= new TH1F("eSeedOverPout2NoCuts", "eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
00179   hadOverEmNoCuts= new TH1F("hadOverEmNoCuts", "Had/EM distribution (Before Cuts)", 600, -2., 2.);
00180 
00181   //Book histograms after ESeed cut
00182   MapCor1ESeed = new TH2F ("MapCor1ESeed", "Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00183   MapCor2ESeed = new TH2F ("MapCor2ESeed", "Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00184   MapCor3ESeed = new TH2F ("MapCor3ESeed", "Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00185   MapCor4ESeed = new TH2F ("MapCor4ESeed", "Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00186   MapCor5ESeed = new TH2F ("MapCor5ESeed", "Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00187   MapCor6ESeed = new TH2F ("MapCor6ESeed", "Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00188   MapCor7ESeed = new TH2F ("MapCor7ESeed", "Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00189   MapCor8ESeed = new TH2F ("MapCor8ESeed", "Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00190   MapCor9ESeed = new TH2F ("MapCor9ESeed", "Correlation  E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00191   MapCor10ESeed = new TH2F ("MapCor10ESeed", "Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00192   MapCor11ESeed = new TH2F ("MapCor11ESeed", "Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00193  
00194   eSeedOverPout2ESeed= new TH1F("eSeedOverPout2ESeed", "eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
00195 
00196   hadOverEmESeed= new TH1F("hadOverEmESeed", "Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
00197  
00198  //Book histograms without any cut
00199   GeneralMap = new TH2F("GeneralMap","Map without any cuts",85,1,85,70,5,75);
00200 
00201   calibClusterSize=ClusterSize_; 
00202   etaMin = mineta_;
00203   etaMax = maxeta_;
00204   phiMin = minphi_;
00205   phiMax = maxphi_;
00206   if(calibAlgo_=="L3") {
00207     MyL3Algo1 = new MinL3Algorithm(keventweight_,calibClusterSize, etaMin, etaMax, phiMin, phiMax);
00208   }else{ 
00209     if(calibAlgo_=="HH" || calibAlgo_=="HHReg"){
00210       MyHH = new HouseholderDecomposition(calibClusterSize, etaMin,etaMax, phiMin, phiMax); 
00211     }else{ 
00212       std::cout<<" Name of Algorithm is not recognize "<<calibAlgo_<<" Should be either L3, HH or HHReg. Abort! "<<std::endl;
00213     }
00214   }
00215   read_events=0;
00216   
00217   // get Region to be calibrated  
00218   ReducedMap = calibCluster.getMap(etaMin, etaMax, phiMin, phiMax);
00219   
00220   oldCalibs.resize(ReducedMap.size(),0.);
00221 
00222    // table is set to zero
00223   for (int phi=0; phi<360; phi++){for (int eta=0; eta<171; eta++){eventcrystal[eta][phi]=0;}}
00224  
00225 
00226   std::cout<<" Begin JOB "<<std::endl;
00227 }
00228 
00229 
00230 //========================================================================
00231 
00232 void
00233 ElectronCalibration::endJob() {
00234 //========================================================================
00235 
00236 int nIterations =10;
00237 if(calibAlgo_=="L3"){ 
00238   solution = MyL3Algo1->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, nIterations);
00239  }else{
00240   if(calibAlgo_=="HH"){
00241     solution = MyHH->iterate(EventMatrix, MaxCCeta, MaxCCphi,EnergyVector, 1,false);
00242   }else{
00243     if(calibAlgo_=="HHReg"){solution = MyHH->runRegional(EventMatrix, MaxCCeta, MaxCCphi,EnergyVector, 2);
00244     }else{ 
00245       std::cout<<" Calibration not run due to problem in Algo Choice..."<<std::endl;return;
00246     }
00247   }
00248  }
00249  for (int ii=0;ii<(int)solution.size();ii++)
00250    {
00251      std::cout << "solution[" << ii << "] = " << solution[ii] << std::endl;
00252      calibs->Fill(solution[ii]); 
00253    }
00254  
00255  newCalibs.resize(ReducedMap.size(),0.);
00256  
00257  calibXMLwriter write_calibrations;
00258  
00259  FILE* MisCalib;
00260  MisCalib = fopen(miscalibfile_.c_str(),"r");
00261  int fileStatus=1;
00262  int eta=-1;
00263  int phi=-1;
00264  float coeff=-1;
00265  
00266  std::map<EBDetId,float> OldCoeff;
00267  
00268  while(fileStatus != EOF) {
00269    fileStatus = fscanf(MisCalib,"%d %d %f\n",  &eta,&phi,&coeff);
00270    if(eta!=-1&&phi!=-1&& coeff!=-1){
00271      //      std::cout<<" We have read correctly the coefficient " << coeff << " corresponding to eta "<<eta<<" and  phi "<<phi<<std::endl;
00272      OldCoeff.insert(std::make_pair(EBDetId(eta,phi,EBDetId::ETAPHIMODE),coeff )); 
00273    }
00274  } 
00275  
00276  fclose(MisCalib);
00277  
00278  int icry=0;
00279  CalibrationCluster::CalibMap::iterator itmap;
00280  for (itmap=ReducedMap.begin(); itmap != ReducedMap.end();itmap++){
00281    
00282    newCalibs[icry] = solution[icry];
00283    
00284    write_calibrations.writeLine(itmap->first,newCalibs[icry]);
00285    float Compare =1.;   
00286    std::map<EBDetId,float>::iterator iter = OldCoeff.find(itmap->first);
00287    if( iter != OldCoeff.end() )Compare = iter->second;
00288 
00289    if((itmap->first).ieta()>mineta_ && (itmap->first).ieta()<maxeta_ && (itmap->first).iphi()>minphi_ && (itmap->first).iphi()<maxphi_){
00290      Map3Dcalib->Fill((itmap->first).ieta(),(itmap->first).iphi(),newCalibs[icry]*Compare ) ;
00291      MisCalibs->Fill(Compare);
00292 
00293 }
00294    if((itmap->first).ieta()< mineta_+2){icry++; continue;}
00295    if((itmap->first).ieta()> maxeta_-2){icry++; continue;}
00296    if((itmap->first).iphi()< minphi_+2){icry++; continue;} 
00297    if((itmap->first).iphi()> maxphi_-2){icry++; continue;}
00298 
00299    calibinter->Fill(newCalibs[icry]);
00300    DiffCalibs->Fill(newCalibs[icry]-1./Compare);
00301    RatioCalibs->Fill(newCalibs[icry]*Compare);
00302    MapCorCalib->Fill(1./Compare, newCalibs[icry]);
00303    icry++;
00304  }
00305  
00306  if(calibAlgo_=="L3"){
00307    solutionNoCuts = MyL3Algo1->iterate(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts,EnergyVectorNoCuts,nIterations);
00308  }else{
00309    if(calibAlgo_=="HH"){
00310      solutionNoCuts = MyHH->iterate(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts, EnergyVectorNoCuts, 1,false);
00311    }else{ 
00312      if(calibAlgo_=="HHReg"){
00313        solutionNoCuts = MyHH->runRegional(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts,EnergyVectorNoCuts, 2);
00314      }else{
00315        std::cout<<" Calibration not run due to problem in AlgoChoice..."<<std::endl;return;
00316      }
00317    }
00318  }
00319  for (int ii=0;ii<(int)solutionNoCuts.size();ii++){
00320    calibsNoCuts->Fill(solutionNoCuts[ii]); 
00321  }
00322  int icryp=0;
00323  CalibrationCluster::CalibMap::iterator itmapp;
00324  for (itmapp=ReducedMap.begin(); itmapp != ReducedMap.end();itmapp++){
00325    
00326    newCalibs[icryp] = solutionNoCuts[icryp];
00327    float Compare2 =1.;   
00328    std::map<EBDetId,float>::iterator iter2 = OldCoeff.find(itmapp->first);
00329    if( iter2 != OldCoeff.end() )Compare2 = iter2->second;
00330    
00331    if((itmapp->first).ieta()>mineta_ && (itmapp->first).ieta()<maxeta_ && (itmapp->first).iphi()>minphi_ && (itmapp->first).iphi()<maxphi_)Map3DcalibNoCuts->Fill((itmapp->first).ieta(),(itmapp->first).iphi(),newCalibs[icryp]*Compare2) ;
00332    if ((itmapp->first).ieta()< mineta_+2){icryp++; continue;}
00333    if ((itmapp->first).ieta()> maxeta_-2){icryp++; continue;}
00334    if ((itmapp->first).iphi()< minphi_+2){icryp++; continue;} 
00335    if ((itmapp->first).iphi()> maxphi_-2){icryp++; continue;}
00336    calibinterNoCuts->Fill(newCalibs[icryp]);
00337    DiffCalibsNoCuts->Fill(newCalibs[icryp]-1./(Compare2));
00338    RatioCalibsNoCuts->Fill(newCalibs[icryp]*Compare2);
00339    MapCorCalibNoCuts->Fill(1./Compare2 ,newCalibs[icryp]);
00340    icryp++;
00341  }
00342  
00343  
00344  
00346  
00347  std::cout << " " << std::endl;
00348  std::cout << "************* STATISTICS **************" << std::endl;
00349  std::cout << " Events Studied "<<read_events<< std::endl;
00350  
00352    
00353     f->Write();
00354    
00355     f->Close();
00356 }
00357 
00358 
00359 //=================================================================================
00360 EBDetId ElectronCalibration::findMaxHit(edm::Handle<EBRecHitCollection> &  phits) {
00361 //=================================================================================
00362 
00363      EcalRecHitCollection ecrh = *phits;
00364      EcalRecHitCollection::iterator it;
00365      int count=0;
00366      EBDetId save;
00367      float en_save=0;
00368      for (it = ecrh.begin(); it != ecrh.end(); it++)
00369      {
00370        EBDetId p = EBDetId(it->id().rawId());
00371         if(it->energy()> en_save){
00372           en_save=it->energy();
00373           save=p;
00374         }
00375         count++;
00376      }
00377      return save;
00378 
00379 }
00380 
00381 //=================================================================================
00382 EBDetId  ElectronCalibration::findMaxHit2(const std::vector<DetId> & v1,const EBRecHitCollection* hits) {
00383 //=================================================================================
00384 
00385   double currEnergy = 0.;
00386   EBDetId maxHit;
00387   
00388   for( std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
00389     if(idsIt->subdetId()!=1) continue;
00390     EBRecHitCollection::const_iterator itrechit;
00391     itrechit = hits->find(*idsIt);
00392            
00393     if(itrechit == hits->end()){
00394       std::cout << "ElectronCalibration::findMaxHit2: rechit not found! " << std::endl;
00395       continue;
00396     }
00397     if(itrechit->energy() > currEnergy) {
00398       currEnergy=itrechit->energy();
00399       maxHit= *idsIt;
00400     }
00401   }
00402   
00403       return maxHit;
00404 }
00405 
00406 
00407 //=================================================================================
00408 void ElectronCalibration::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){
00409 //=================================================================================
00410    using namespace edm;
00411    
00412    // Get EBRecHits
00413    Handle<EBRecHitCollection> phits;
00414    iEvent.getByLabel( recHitLabel_, phits);
00415    if (!phits.isValid()) {
00416      std::cerr << "Error! can't get the product EBRecHitCollection: " << std::endl;
00417    }
00418 
00419    const EBRecHitCollection* hits = phits.product(); // get a ptr to the product
00420 
00421    // Get pixelElectrons
00422    Handle<reco::GsfElectronCollection> pElectrons;
00423 
00424    iEvent.getByLabel(electronLabel_, pElectrons);
00425    if (!pElectrons.isValid()) {
00426      std::cerr << "Error! can't get the product ElectronCollection: " << std::endl;
00427    }
00428 
00429   const reco::GsfElectronCollection* electronCollection = pElectrons.product();
00430   read_events++;
00431   if(read_events%1000 ==0)std::cout << "read_events = " << read_events << std::endl;
00432   
00433   if(!hits)return;
00434   if(hits->size() == 0)return;
00435   if(!electronCollection)return;
00436   if(electronCollection->size() == 0)return;
00437   
00438   
00440   //                          START HERE....
00442   reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();
00443 
00444   reco::GsfElectron highPtElectron;
00445 
00446   float highestElePt=0.;
00447   bool found=false;
00448   for (eleIt=electronCollection->begin(); eleIt!=electronCollection->end(); eleIt++) {
00449     //Comments
00450     if(fabs(eleIt->eta())>(maxeta_+3) * 0.0175) continue;
00451     if(eleIt->eta()<(mineta_-3) * 0.0175) continue;
00452 
00453      if(eleIt->pt()>highestElePt) {
00454        highestElePt=eleIt->pt();
00455        highPtElectron = *eleIt;
00456        found =true;
00457      }
00458 
00459   }
00460   if(highestElePt<ElePt_)return;
00461       if(!found) return;
00462       const reco::SuperCluster & sc = *(highPtElectron.superCluster()) ;
00463       if(fabs(sc.eta())>(maxeta_+3) * 0.0175){
00464         std::cout<<"++++ Problem with electron, electron eta is "<< highPtElectron.eta()<<" while SC is "<<sc.eta()<<std::endl;return;
00465       }
00466 //      std::cout << "track eta = " << highPtElectron.eta() << std::endl;
00467 //      std::cout << "track phi = " << highPtElectron.phi() << std::endl;
00468     
00469       std::vector<DetId> v1;
00470       //Loop to fill the vector of DetIds
00471 for (std::vector<std::pair<DetId,float> >::const_iterator idsIt = sc.hitsAndFractions().begin();
00472        idsIt != sc.hitsAndFractions().end ();++idsIt)
00473   {v1.push_back(idsIt->first);
00474  }
00475 
00476       //getHitsByDetId(); //Change function name
00477       EBDetId maxHitId;
00478       
00479       maxHitId = findMaxHit2(v1,hits); 
00480       
00481       if(maxHitId.null()){std::cout<<" Null "<<std::endl; return;}
00482       
00483       int maxCC_Eta = maxHitId.ieta();
00484       int maxCC_Phi = maxHitId.iphi();
00485       
00486       if(maxCC_Eta>maxeta_ )return;
00487       if(maxCC_Eta<mineta_ )return;
00488       if(maxCC_Phi>maxphi_ ) return;
00489       if(maxCC_Phi<minphi_ )  return;
00490 
00491       // number of events per crystal is set
00492       if(numevent_>0){
00493         eventcrystal[maxCC_Eta+85][maxCC_Phi-1]+=1;
00494         if (eventcrystal[maxCC_Eta+85][maxCC_Phi-1] > numevent_) return;
00495       }
00496       
00497       std::vector<EBDetId> Xtals5x5 = calibCluster.get5x5Id(maxHitId);
00498       
00499       if((int)Xtals5x5.size()!=ClusterSize_*ClusterSize_)return;
00500  
00501       // fill cluster energy
00502       std::vector<float> energy;
00503       float energy3x3=0.;  
00504       float energy5x5=0.;  
00505       
00506       for (int icry=0;icry<ClusterSize_*ClusterSize_;icry++){
00507         
00508            EBRecHitCollection::const_iterator itrechit;
00509            if(Xtals5x5[icry].subdetId()!=1) continue;
00510            itrechit = hits->find(Xtals5x5[icry]);
00511            if(itrechit==hits->end())
00512              { std::cout << "DetId not is e25" << std::endl;
00513                continue;
00514              }
00515            
00516            if (edm::isNotFinite(itrechit->energy())) return;      
00517            energy.push_back(itrechit->energy());
00518            energy5x5 += energy[icry];
00519            
00520            if ( icry == 6  || icry == 7  || icry == 8 ||
00521                 icry == 11 || icry == 12 || icry ==13 ||
00522                 icry == 16 || icry == 17 || icry ==18   )
00523              {
00524                energy3x3+=energy[icry];
00525              }
00526            
00527       }
00528       if((int)energy.size()!=ClusterSize_*ClusterSize_) return;
00529       //Once we have the matrix 5x5, we have to correct for gaps/cracks/umbrella and maincontainement  
00530       
00531       GeneralMap->Fill(maxCC_Eta,maxCC_Phi);
00532       
00533       EoP_all->Fill(highPtElectron.eSuperClusterOverP()); 
00534       
00535       if(highPtElectron.classification()==elecclass_ || elecclass_== -1 ){
00536         
00537         float Ptrack_in=sqrt( pow(highPtElectron.trackMomentumAtVtx().X(),2) +pow(highPtElectron.trackMomentumAtVtx().Y(),2) + pow(highPtElectron.trackMomentumAtVtx().Z(),2) );
00538         
00539         float UncorrectedPatCalo = sqrt(pow(highPtElectron.trackMomentumAtCalo().X(),2)+pow(highPtElectron.trackMomentumAtCalo().Y(),2)+pow(highPtElectron.trackMomentumAtCalo().Z(),2));
00540         
00541         float Ptrack_out = sqrt( pow(highPtElectron.trackMomentumOut().X(),2)+ pow(highPtElectron.trackMomentumOut().Y(),2)+ pow(highPtElectron.trackMomentumOut().Z(),2) );
00542         
00543         
00544         EventMatrixNoCuts.push_back(energy);
00545         EnergyVectorNoCuts.push_back(UncorrectedPatCalo);
00546         
00547         MaxCCetaNoCuts.push_back(maxCC_Eta);
00548         MaxCCphiNoCuts.push_back(maxCC_Phi);
00549         
00550         WeightVectorNoCuts.push_back(energy5x5/UncorrectedPatCalo);
00551         
00552         //---------------------------------------------------No Cuts-------------------------------------------------------
00553         e9NoCuts->Fill(energy3x3); 
00554         e25NoCuts->Fill(energy5x5); 
00555         e9Overe25NoCuts->Fill(energy3x3/energy5x5);
00556         scENoCuts->Fill(sc.energy()); 
00557         
00558         trPNoCuts->Fill(UncorrectedPatCalo); 
00559         
00560         EoPNoCuts->Fill(highPtElectron.eSuperClusterOverP()); 
00561         e25OverScENoCuts->Fill(energy5x5/sc.energy());
00562         
00563         E25oPNoCuts->Fill(energy5x5/UncorrectedPatCalo);
00564         
00565         MapNoCuts->Fill(maxCC_Eta,maxCC_Phi);
00566         PinOverPoutNoCuts->Fill( sqrt( pow(highPtElectron.trackMomentumAtVtx().X(),2) +pow(highPtElectron.trackMomentumAtVtx().Y(),2) + pow(highPtElectron.trackMomentumAtVtx().Z(),2) )/sqrt( pow(highPtElectron.trackMomentumOut().X(),2)+ pow(highPtElectron.trackMomentumOut().Y(),2)+ pow(highPtElectron.trackMomentumOut().Z(),2) ) );
00567         eSeedOverPoutNoCuts->Fill(highPtElectron.eSuperClusterOverP());
00568         
00569         MapCor1NoCuts->Fill(energy5x5/UncorrectedPatCalo,energy5x5/Ptrack_in);
00570         MapCor2NoCuts->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSuperClusterOverP());
00571         MapCor3NoCuts->Fill(energy5x5/UncorrectedPatCalo,Ptrack_out/Ptrack_in);
00572         MapCor4NoCuts->Fill(energy5x5/UncorrectedPatCalo,energy5x5/highPtElectron.p());
00573         MapCor5NoCuts->Fill(energy5x5/UncorrectedPatCalo,UncorrectedPatCalo/Ptrack_out);
00574         MapCor6NoCuts->Fill(Ptrack_out/Ptrack_in,energy5x5/Ptrack_in);
00575         MapCor7NoCuts->Fill(Ptrack_out/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00576         MapCor8NoCuts->Fill(energy5x5/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00577         MapCor9NoCuts->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSeedClusterOverPout());
00578         MapCor10NoCuts->Fill(highPtElectron.eSeedClusterOverPout(),Ptrack_out/Ptrack_in);
00579         MapCor11NoCuts->Fill(highPtElectron.eSeedClusterOverPout(),energy5x5/Ptrack_in);
00580         
00581         PinMinPoutNoCuts->Fill((Ptrack_in-Ptrack_out)/Ptrack_in);
00582         
00583         Error1NoCuts-> Fill(highPtElectron.trackMomentumError()/Ptrack_in);
00584         Error2NoCuts->Fill(highPtElectron.trackMomentumError()/Ptrack_out);
00585         Error3NoCuts->Fill(highPtElectron.trackMomentumError()/UncorrectedPatCalo);
00586         eSeedOverPout2NoCuts->Fill(highPtElectron.eSeedClusterOverPout());
00587         
00588         hadOverEmNoCuts->Fill(highPtElectron.hadronicOverEm());
00589         
00590         //------------------------------------------------Cuts-----------------------------------------------------
00591         //Cuts!
00592         if((energy3x3/energy5x5)<cut1_)return;
00593         
00594         if((Ptrack_out/Ptrack_in)< cut2_  || (Ptrack_out/Ptrack_in)> cut3_ )return;
00595         if((energy5x5/Ptrack_in)< cutEPin1_  || (energy5x5/Ptrack_in)> cutEPin2_ )return;
00596         
00597         e9->Fill(energy3x3); 
00598         e25->Fill(energy5x5); 
00599         e9Overe25->Fill(energy3x3/energy5x5);
00600         scE->Fill(sc.energy()); 
00601         trP->Fill(UncorrectedPatCalo);
00602         
00603         EoP->Fill(highPtElectron.eSuperClusterOverP()); 
00604         e25OverScE->Fill(energy5x5/sc.energy());
00605         
00606         E25oP->Fill(energy5x5/UncorrectedPatCalo);
00607         
00608         Map->Fill(maxCC_Eta,maxCC_Phi);
00609         PinOverPout->Fill( sqrt( pow(highPtElectron.trackMomentumAtVtx().X(),2) +pow(highPtElectron.trackMomentumAtVtx().Y(),2) + pow(highPtElectron.trackMomentumAtVtx().Z(),2) )/sqrt( pow(highPtElectron.trackMomentumOut().X(),2)+ pow(highPtElectron.trackMomentumOut().Y(),2)+ pow(highPtElectron.trackMomentumOut().Z(),2) ) );
00610         eSeedOverPout->Fill(highPtElectron.eSuperClusterOverP());
00611         
00612         MapCor1->Fill(energy5x5/UncorrectedPatCalo,energy5x5/Ptrack_in);
00613         MapCor2->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSuperClusterOverP());
00614         MapCor3->Fill(energy5x5/UncorrectedPatCalo,Ptrack_out/Ptrack_in);
00615         MapCor4->Fill(energy5x5/UncorrectedPatCalo,energy5x5/highPtElectron.p());
00616         MapCor5->Fill(energy5x5/UncorrectedPatCalo,UncorrectedPatCalo/Ptrack_out);
00617         MapCor6->Fill(Ptrack_out/Ptrack_in,energy5x5/Ptrack_in);
00618         MapCor7->Fill(Ptrack_out/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00619         MapCor8->Fill(energy5x5/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00620         MapCor9->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSeedClusterOverPout());
00621         MapCor10->Fill(highPtElectron.eSeedClusterOverPout(),Ptrack_out/Ptrack_in);
00622         MapCor11->Fill(highPtElectron.eSeedClusterOverPout(),energy5x5/Ptrack_in);
00623         
00624         PinMinPout->Fill((Ptrack_in-Ptrack_out)/Ptrack_in);
00625         
00626         Error1-> Fill(highPtElectron.trackMomentumError()/Ptrack_in);
00627         Error2->Fill(highPtElectron.trackMomentumError()/Ptrack_out);
00628         Error3->Fill(highPtElectron.trackMomentumError()/UncorrectedPatCalo);
00629         
00630         eSeedOverPout2->Fill(highPtElectron.eSeedClusterOverPout());
00631         hadOverEm->Fill(highPtElectron.hadronicOverEm());
00632         
00633         EventMatrix.push_back(energy);
00634         EnergyVector.push_back(UncorrectedPatCalo);
00635         MaxCCeta.push_back(maxCC_Eta);
00636         MaxCCphi.push_back(maxCC_Phi);
00637         
00638         WeightVector.push_back(energy5x5/UncorrectedPatCalo);
00639         
00640         //-------------------------------------------------------Extra Cut-----------------------------------------------------
00641         if(highPtElectron.eSeedClusterOverPout()< cutESeed_ ) return;
00642 
00643         MapCor1ESeed->Fill(energy5x5/UncorrectedPatCalo,energy5x5/Ptrack_in);
00644         MapCor2ESeed->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSuperClusterOverP());
00645         MapCor3ESeed->Fill(energy5x5/UncorrectedPatCalo,Ptrack_out/Ptrack_in);
00646         MapCor4ESeed->Fill(energy5x5/UncorrectedPatCalo,energy5x5/highPtElectron.p());
00647         MapCor5ESeed->Fill(energy5x5/UncorrectedPatCalo,UncorrectedPatCalo/Ptrack_out);
00648         MapCor6ESeed->Fill(Ptrack_out/Ptrack_in,energy5x5/Ptrack_in);
00649         MapCor7ESeed->Fill(Ptrack_out/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00650         MapCor8ESeed->Fill(energy5x5/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00651         MapCor9ESeed->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSeedClusterOverPout());
00652         MapCor10ESeed->Fill(highPtElectron.eSeedClusterOverPout(),Ptrack_out/Ptrack_in);
00653         MapCor11ESeed->Fill(highPtElectron.eSeedClusterOverPout(),energy5x5/Ptrack_in);
00654         
00655         eSeedOverPout2ESeed->Fill(highPtElectron.eSeedClusterOverPout());
00656         
00657         hadOverEmESeed->Fill(highPtElectron.hadronicOverEm());
00658         
00659       }else{return;}
00660 }
00661
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/Calibration/EcalCalibAlgos/src/ElectronCalibration.cc
