CAF: ApplyJES.cpp Source File

00001 #include <iostream>
00002 #include "caf_util/ApplyJES.hpp"
00003 #include "cafe/Config.hpp"
00004 #include "jetcorr/JetDefs.hpp"
00005 #include "jetcorr/CalTool.hpp"
00006 
00007 using namespace cafe;
00008 using namespace std;
00009 
00010 namespace caf_util {
00011 
00012   ApplyJES::ApplyJES(const char *name) : cafe::SelectUserObjects<TMBJet>(name), _event(0),_vars("_nchunks")
00013   {
00014     cafe::Config config(name);
00015     _doJESMU    = config.get("doJESMU",true);
00016     _jetBranch  = config.get("From",""); 
00017     _muonBranch = config.get("muonBranch",""); 
00018     _jes_datatype = config.get("JetEnergyScale_datatype",-1); // default: fill from input event
00019     _jes_jetalgo = config.get("JetEnergyScale_jetalgo",2); // default: JCCB
00020     _addDataErr2MC = config.get("AddDataErrorsToMC", false) ;
00021 
00022     if ((_jetBranch=="JCCA" && _jes_jetalgo == 2) || (_jetBranch=="JCCB" && _jes_jetalgo == 1)) {
00023         cout << " WARNING: in ApplyJES, we are using jet branch " << _jetBranch 
00024              << " as input while asking for JES correction of type " << _jes_jetalgo 
00025              << ". This is not compatible !" << endl;
00026     }
00027 
00028     _muotag_dr       = config.get("MuoTagdR",0.5);
00029     _muotag_ptrel    = config.get("MuoTagPtrel",0.);
00030 
00031     float energy = 50.0;
00032     float etap   = 0.0;
00033     float etad   = 0.0;
00034     float phid   = 0.0;
00035     _pjc = ParticleJetCorr(energy,etap,etad, phid, _jes_jetalgo, 0, 1, 0., 0);
00036 
00037     _debug = config.get("Debug", 0);
00038     if(_muonBranch.empty()) { 
00039         _muonBranch = "Muon";
00040         _useDefaultMuon = true;
00041     }
00042 
00043     _p1803cafPatch = config.get("p1803CAFPatch",true) ; 
00044 
00045   }
00046 
00047   bool ApplyJES::processEvent(cafe::Event &event)
00048   {
00049     if (_jes_datatype == -1) {
00050       _isMC=false;
00051       if (event.getMCEventInfo(_vars)!=NULL) {
00052         if (event.getMCEventInfo(_vars)->nchunks() > 0) _isMC=true;     
00053       }
00054       if (_isMC) _jes_datatype = 2 ; else _jes_datatype = 0;    
00055     } 
00056 
00057     _event = &event ;
00058     SelectUserObjects<TMBJet>::processEvent(event) ;
00059     return true  ;
00060   }
00061 
00062   bool ApplyJES::selectObject(const TMBJet &jet)
00063   {
00064     return true ;
00065   }
00066 
00067   void ApplyJES::after(cafe::Collection<TMBJet>& input, cafe::Collection<TMBJet>& rejected)
00068   {
00069     Collection<TMBMuon> muoncol(_event->getCollection<TMBMuon>(_muonBranch.c_str()));
00070 
00071     Collection<TMBPrimaryVertex> vertexcol(_event->getPrimaryVertices());
00072     int npv = 0;
00073     for (Collection<TMBPrimaryVertex>::const_iterator it_vtx = vertexcol.begin() ; 
00074         it_vtx != vertexcol.end() ; ++ it_vtx) {
00075         if (it_vtx->ntrack() > 2) npv++;
00076     } // for vertex
00077 
00078     int ijet = 0;
00079     if (_debug) cout << " ApplyJES: njets=" << input.size() << " nmuon=" << muoncol.size() << " npv=" << npv << endl;
00080 
00081     for (Collection<TMBJet>::iterator it_jet = input.begin() ;
00082          it_jet != input.end() ; ++ it_jet)  {
00083 
00084       // apply p18.03 patch
00085       if (_p1803cafPatch && it_jet->detEta() == 999 && it_jet->detPhi() == 999
00086           && it_jet->hadec() == 999 && it_jet->hadcc() == 999 && it_jet->ccmg() == 999
00087           && it_jet->Eta() < 100 && it_jet->Phi() < 100) {
00088         //         cerr << "WARNING: Doing something nasty in ApplyJES::before (etaD = 999.) !!! " << endl;
00089         it_jet->Set1 (it_jet->charge(), 64.0, 0.0, it_jet->emf(), it_jet->em1f(),
00090                   it_jet->em2f(),it_jet->em3f(), it_jet->ccmg(), it_jet->icdf(), it_jet->ecmg(), it_jet->icrf(),
00091                   it_jet->fh1f(), it_jet->fh2f(), it_jet->fh3f(), it_jet->chETfraction(), it_jet->emcc(),
00092                   it_jet->hadcc(), it_jet->emec(), it_jet->hadec(), it_jet->hot(), it_jet->mxET(), it_jet->cpsE(),
00093                   it_jet->etaW(), it_jet->phiW(), it_jet->sET(), it_jet->vPT(), it_jet->iPT(), it_jet->seedET(),
00094                   it_jet->pxCH(), it_jet->pyCH(), it_jet->pzCH() );
00095       }
00096     
00097     
00098     if ( it_jet->E() == 0. ) continue;
00099       Int_t actas = it_jet->CurrentlyActAs();
00100       it_jet->ActAsUnCorrected();
00101       
00102       if (_debug) {
00103         TMBJet jet = *it_jet;
00104         double pt = jet.Pt();
00105         jet.ActAsJESCorrected();
00106         double ptcorr = jet.Pt();
00107         jet.ActAsJESMUCorrected();
00108         double ptcorrmu = jet.Pt();
00109                                                                                                                                       
00110         cout << "   ApplyJES: -- old jet: ";
00111         cout << " uncorrected pt=" << pt << " , jes=" << it_jet->JES_C() << " ptcorr=" << ptcorr << " ptcorrmu=" << ptcorrmu << " , eta=" << it_jet->Eta() << " isGood=" << it_jet->isGood() <<  endl;
00112       } // debug
00113 
00114       float etad = CalTool::EtaDToEtaDet(it_jet->detEta());
00115       float phid = CalTool::PhiDToPhiDet(it_jet->detPhi());
00116 
00117       _pjc.reset(it_jet->E(),it_jet->Eta(), etad, phid, _jes_jetalgo, _jes_datatype,npv, it_jet->emf(), it_jet->split_merge_word());
00118       _pjc.calc();
00119       double Enew;
00120       float dummy, C, dC_stat_up, dC_sys_up, dC_stat_down, dC_sys_down, metC, metdC_stat, metdC_sys;
00121       _pjc.forJets(&C, &dummy, &dummy, &Enew);
00122       _pjc.error_up(&dC_stat_up,&dC_sys_up);
00123       _pjc.error_down(&dC_stat_down,&dC_sys_down);
00124       _pjc.forMET(&metC, &metdC_stat, &metdC_sys);
00125 
00126       if (_addDataErr2MC && _jes_datatype == 2) {
00127         _pjc.reset(it_jet->E(),it_jet->Eta(), etad, phid, _jes_jetalgo, 
00128                    0,npv, it_jet->emf(), it_jet->split_merge_word());
00129         _pjc.calc();
00130         double Enew1;
00131         float dummy, C1, dC_stat_up1, dC_sys_up1, dC_stat_down1, dC_sys_down1, metC1, metdC_stat1, metdC_sys1;
00132         _pjc.forJets(&C1, &dummy, &dummy, &Enew1);
00133         _pjc.error_up(&dC_stat_up1,&dC_sys_up1);
00134         _pjc.error_down(&dC_stat_down1,&dC_sys_down1);
00135         _pjc.forMET(&metC1, &metdC_stat1, &metdC_sys1);
00136         
00137         dC_sys_up = C*sqrt(pow(dC_stat_up1/C1,2)  + pow(dC_sys_up1/C1,2) + pow(dC_sys_up/C,2)) ;
00138         dC_sys_down = C*sqrt(pow(dC_stat_down1/C1,2)  + pow(dC_sys_down1/C1,2) + pow(dC_sys_down/C,2)) ;
00139         metdC_sys = metC*sqrt(pow(metdC_stat1/metC1,2)  + pow(metdC_sys1/metC1,2) + pow(metdC_sys/metC,2)) ;
00140       }
00141 
00142 
00143       it_jet->SetCorr(C, dC_stat_up, dC_sys_up, dC_stat_down, dC_sys_down, metC, metdC_stat, metdC_sys);
00144       it_jet->SetCorrMU(C, dC_stat_up, dC_sys_up, dC_stat_down, dC_sys_down, metC, metdC_stat, metdC_sys);
00145 
00146       if (_debug) cout << " APPLYJES: C=" << C << endl;
00147  
00148       // muon matching
00149 
00150       if (!_doJESMU) continue;
00151 
00152       vector<const TMBMuon*>  temp;
00153 
00154       for (Collection<TMBMuon>::const_iterator muptr = muoncol.begin() ;
00155            muptr != muoncol.end() ; ++ muptr)  {
00156         if (_useDefaultMuon) {
00157            if (!muptr->isMedium()) continue;
00158         }
00159 
00160         float muonpx  = muptr->Px();
00161         float muonpy  = muptr->Py();
00162         float muonpz  = muptr->Pz();
00163         float muonphi = muptr->Phi();
00164         float muoneta = muptr->Eta();
00165         float muonpt  = muptr->Pt();
00166 
00167         // calculate dr and store all muons close to the jet
00168         float de = fabs( it_jet->Eta() - muoneta );
00169         float dp = kinem::delta_phi( it_jet->Phi(), muonphi);
00170         float dr = sqrt(dp*dp+de*de);
00171         if (_debug) cout << " in APPLYJES: dR(mu,jet)=" << dr << endl;
00172         if ( dr > _muotag_dr) continue;
00173         temp.push_back(&(*muptr));
00174         if (_debug) cout << " in APPLYJES: muon pt=" << muptr->Pt() << " eta=" << muptr->Eta() << " nseg=" << muptr->nseg() << endl;
00175 
00176      } // for muon
00177 
00178      const TMBMuon* mu;
00179      float pTrel = 0;
00180      if (temp.size()>0) {
00181         bool hasMU = false;
00182         const TMBMuon* mu=0;
00183         float pTrel = 0;
00184         for (int i=0; i<temp.size(); i++) {  
00185            double pTrel_curr = PTrel(temp[i], *it_jet);
00186            // the possibility to cut on ptrel is not implemented into d0correct
00187            if (pTrel_curr >= pTrel && pTrel_curr > _muotag_ptrel) {
00188               mu = temp[i];
00189               pTrel = pTrel_curr;
00190               hasMU = true;
00191            } // ptrel
00192         } // for temp
00193 
00194         if (hasMU) {
00195           it_jet->SetFlag(it_jet->isGood(), it_jet->isL1Conf(), it_jet->isBad(), it_jet->isNoise(), it_jet->isEM(), true, it_jet->isSmeared(), it_jet->isSmearedMU());
00196           _pjc.addMuon(mu->Px(), mu->Py(), mu->Pz(), mu->E(), BJET, pTrel);
00197           _pjc.calc();
00198           double muEnew;
00199           float mudummy, muC, mudC_stat_up, mudC_sys_up, mudC_stat_down, mudC_sys_down, mumetC, mumetdC_stat, mumetdC_sys;
00200           _pjc.forJets(&muC, &mudummy, &mudummy, &muEnew);
00201           _pjc.error_up(&mudC_stat_up,&mudC_sys_up);
00202           _pjc.error_down(&mudC_stat_down,&mudC_sys_down);
00203           _pjc.forMET(&mumetC, &mumetdC_stat, &mumetdC_sys);
00204           if (_debug) cout << " APPLYJES: muC=" << C << endl;
00205           it_jet->SetCorrMU(muC, mudC_stat_up, mudC_sys_up, mudC_stat_down, mudC_sys_down, mumetC, mumetdC_stat, mumetdC_sys);
00206         } // hasMU
00207      } // temp.size()>0
00208 
00209 
00210      if (actas==TMBJet::kUnCorrected) it_jet->ActAsUnCorrected();
00211      if (actas==TMBJet::kJESCorrected) it_jet->ActAsJESCorrected();
00212      if (actas==TMBJet::kJESCorrectedShiftedPlus) it_jet->ActAsJESCorrectedShiftedPlus();
00213      if (actas==TMBJet::kJESCorrectedShiftedMinus) it_jet->ActAsJESCorrectedShiftedMinus();
00214      if (actas==TMBJet::kJESMUCorrected) it_jet->ActAsJESMUCorrected();
00215      if (actas==TMBJet::kJESMUCorrectedShiftedPlus) it_jet->ActAsJESMUCorrectedShiftedPlus();
00216      if (actas==TMBJet::kJESMUCorrectedShiftedMinus) it_jet->ActAsJESMUCorrectedShiftedMinus();
00217      if (actas==TMBJet::kSmeared) it_jet->ActAsSmeared();
00218      if (actas==TMBJet::kSmearedMU) it_jet->ActAsSmearedMU();
00219 
00220       if (_debug) {
00221 
00222         TMBJet jet = *it_jet;
00223         jet.ActAsUnCorrected();   
00224         double pt = jet.Pt();
00225         jet.ActAsJESCorrected();
00226         double ptcorr = jet.Pt();
00227         jet.ActAsJESMUCorrected();
00228         double ptcorrmu = jet.Pt();
00229 
00230         cout << "   ApplyJES: -- new jet: ";
00231         cout << " uncorrected pt=" << pt << " , jes=" << it_jet->JES_C() << " ptcorr=" << ptcorr << " ptcorrmu=" << ptcorrmu << " , eta=" << it_jet->Eta() << " isGood=" << it_jet->isGood() << endl;
00232       } // debug
00233     } // for jet
00234 
00235   }
00236 
00237   /*****************************************************************************/
00238   // same computation as in JetRecoCorr.cpp
00239   float ApplyJES::PTrel(const TMBMuon* muon, const TMBJet& jet)
00240   {
00241     // components of sum momentum of muon+jet
00242     double mujet_x = muon->Px() + jet.Px();
00243     double mujet_y = muon->Py() + jet.Py();
00244     double mujet_z = muon->Pz() + jet.Pz();
00245     double mujet2     = mujet_x*mujet_x       + mujet_y*mujet_y       + mujet_z*mujet_z;
00246     double muon2      = muon->Px()*muon->Px() + muon->Py()*muon->Py() + muon->Pz()*muon->Pz();
00247     double muonXmujet = muon->Px()*mujet_x    + muon->Py()*mujet_y    + muon->Pz()*mujet_z;
00248     double pLrel2 = muonXmujet*muonXmujet/mujet2;
00249     // muon2 = pTrel2 + pLrel2
00250     double pTrel2 = muon2 - pLrel2;
00251     return (pTrel2 > 0)? sqrt(pTrel2): 0;
00252   }
00253 
00254 } // namespace
00255 
00256 ClassImp(caf_util::ApplyJES);
00257