NtupleAnalysis: NtupleReader/src/MySystematics.cxx Source File

00001 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00002 // january 2008        Ketevi A. Assamagan
00003 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00004 
00005 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
00006 //:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS MySystematics.h ::
00009 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
00010 
00011 //::::::::::::::::::
00012 //:: HEADER FILES ::
00013 //::::::::::::::::::
00014 
00015 #include "NtupleReader/MySystematics.h"
00016 
00017 //::::::::::::::::::::::::
00018 //:: NAMESPACE SETTINGS ::
00019 //::::::::::::::::::::::::
00020 
00021 using namespace std;
00022 
00023 //*****************************************************************************
00024 
00026 void MySystematics::init( unsigned int  seed ) {
00027 
00028    m_random = new TRandom ( seed );  
00029    m_dataset_info = &MyDatasetInfo::Instance();
00030 
00031 
00032 }
00033 
00034 
00036 void MySystematics::destruct( void ) {
00037 
00038    delete m_random;
00039 
00040 }
00041 
00042 // *****************************************************************************
00043 
00044 //:::::::::::::::::::::::::::::::
00045 //:: METHOD  apply_syst        ::
00046 //::   to apply systematics
00047 //::
00048 //::    vp_syst   = vector of particle to investigate
00049 //::    met_syst  = missing energy has to be recalculated after manipulating particles in vp_syst
00050 //::    f_reso    => 1 = resolution should be smeared
00051 //::    f_eScale_minus  => 1 = energy scale should be manipulated in netagtive direction
00052 //::    f_eScale_plus  => 1 = energy scale should be manipulated in positive direction
00053 //::    f_effi    => 1 = effi should be manipulated ( does not chnage met (execpt for muons) )
00054 //::    f_btag    => 1 = b taging effi should be manipulated (USE bjets matched to TRUTH!!)
00055 //::--------------------------------------------------------
00056 //::NOTE: You should do the systematic study one by one!!!
00057 //::--------------------------------------------------------
00058 //::   outcome: vp_syst:    tlv's are updated
00059 //::            met_syst:   missing energy is recalculated (and should be used in analysis)
00060 //:::::::::::::::::::::::::::::::
00061 
00062 void MySystematics::apply_syst( std::vector<MyParticle*> * vp_syst, MyMissingEt * met_syst, int f_reso, int f_eScale_minus, int f_eScale_plus, int f_effi, int f_btag){
00063   
00064 
00065 //   cout<<" syst: reso="<<f_reso<<" e+="<<f_eScale_minus<<" e-="<<f_eScale_plus<<" effi="<<f_effi<<"\n";
00066 //   cout<<"     number of particles = "<<vp_syst->size()<<"\n";
00067 
00068   double m_gev = m_dataset_info->get_scale();
00069 
00070   double vp_m=0;
00071   double egamma_scale_factor_minus =0.995;
00072   double egamma_scale_factor_plus =1.005;
00073   double electron_eff_error=0.002;
00074   double photon_eff_error=0.002;
00075 
00076   double muon_scale_factor_minus =0.99;
00077   double muon_scale_factor_plus =1.01;
00078   double muon_eff_error=0.01;
00079 
00080   double jet_scale_factor_minus =0.93;
00081   double jet_scale_factor_plus =1.07;
00082   double jet_scale_factor_minus_large_eta =0.85;
00083   double jet_scale_factor_plus_large_eta =1.15;
00084   double jet_eff_error=0.0;                       //100% effi
00085   double jet_eff_error_large_eta=0.0;             //100% effi
00086   double jet_METscale_factor_minus =-0.05;
00087   double jet_METscale_factor_plus =0.05;
00088 
00089   double taujet_scale_factor_minus =0.95;
00090   double taujet_scale_factor_plus =1.05;
00091   double tau_eff_error=0.05;
00092   double taujet_METscale_factor_minus =-0.05;
00093   double taujet_METscale_factor_plus =0.05;
00094 
00095   double bjet_eff_error_scale=0.05;
00096   //double bjet_eff=0.6; // should be in pt/eta bins
00097   // double bjet_eff[4][4]={ };
00098   // double bjet_eta[5]={ };
00099   // double bjet_pt[5]={ };
00100   double bjet_eff_error=0.;
00101 
00102   int nsyst=0;
00103   if (f_reso==1) nsyst++;
00104   if (f_eScale_minus==1) nsyst++;
00105   if (f_eScale_plus==1) nsyst++;
00106   if (f_effi==1) nsyst++;
00107   if (f_btag==1) nsyst++;
00108 
00109   if(nsyst>1){
00110     cout<<" ** ERROR MyTools::apply_syst: You should do the systematic study one by one!!! \n";
00111     cout<<" ** ERROR MyTools::apply_syst:  NO EFFEKT HAS BEEN APPLIED!!! \n";
00112     cerr<<" ** ERROR MyTools::apply_syst: You should do the systematic study one by one!!! \n";
00113     cerr<<" ** ERROR MyTools::apply_syst:  NO EFFEKT HAS BEEN APPLIED!!! \n";
00114     return;
00115   }
00116 
00117   double new_metx=met_syst->met_x();
00118   double new_mety=met_syst->met_y();
00119   double new_met=met_syst->met();
00120   double new_metsum=met_syst->met_sum();
00121 
00122   //for (unsigned int il=0;il<vp_syst->size();il++){    
00123   // loop starting from the last entry in the vector
00124   for (int il=vp_syst->size()-1; il>=0; il--){ //because of erase
00125 
00126     bool syst_remove=false;
00127 
00128     TLorentzVector aod_tlv=(((*vp_syst)[il])->tlv());
00129     int id_pdg=0;
00130     id_pdg=TMath::Abs((*vp_syst)[il]->pdgId());
00131 
00132     if (id_pdg==22) vp_m=0;
00133     if (id_pdg==11) vp_m=m_gev*0.511;
00134     if (id_pdg==13) vp_m=m_gev*105.6;
00135     if (id_pdg==15) vp_m=m_gev*1776.99;
00136     /* jets in scale/reso get vp_m=0 */
00137 
00138     TLorentzVector new_tlv=aod_tlv;
00139 
00140     /*--- to get MeV tlv to GeV for resolution --*/
00141     double mev_adapt=1.;
00142     if (m_gev==1.) mev_adapt=0.001;
00143 
00144     /*---------------------------------------------*/
00145 
00146     if( f_reso==1 ){
00147       /* --------------------------------------
00148          systematics in resolution
00149          -------------------------------------- */
00150       
00151       //electron
00152       if (id_pdg==11){
00153 
00154         /* Note that Et in GeV/MeV???? egamma_resolution (linear: 0.01*et) */
00155         double et_sys=new_tlv.Et();
00156         double new_et=egamma_resolution(et_sys);
00157         double new_e=new_et*cosh(new_tlv.Eta());
00158         double new_p=new_e*new_e-vp_m*vp_m;
00159         new_p=sqrt(new_p);
00160         double new_pt=new_p/cosh(new_tlv.Eta());
00161 
00162         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00163         
00164       }
00165       //photon
00166       if (id_pdg==22){
00167 
00168         /* Note that Et in GeV/MeV???? egamma_resolution (linear: 0.01*et) */
00169         double et_sys=new_tlv.Et();
00170         double new_et=egamma_resolution(et_sys);
00171         double new_e=new_et*cosh(new_tlv.Eta());
00172         double new_p=new_e*new_e-vp_m*vp_m;
00173         new_p=sqrt(new_p);
00174         double new_pt=new_p/cosh(new_tlv.Eta());
00175 
00176         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00177         
00178       }
00179       //muon
00180       if (id_pdg==13){
00181 
00182         /* Note that pT is in GeV in muon_resolution */
00183         double pt_sys=new_tlv.Pt();
00184         double new_pt=muon_resolution(mev_adapt*pt_sys)/mev_adapt;
00185         double new_p=new_pt*cosh(new_tlv.Eta());
00186         double new_e=new_p*new_p+vp_m*vp_m;
00187         new_e=sqrt(new_e);
00188 
00189         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00190         
00191       }
00192       //tau-jet
00193       if (id_pdg==15){
00194 
00195         /* Note that E is in GeV in jet_resolution */
00196         double e_sys=new_tlv.E();
00197         double new_e=taujet_resolution(mev_adapt*e_sys)/mev_adapt;
00198         double new_p=new_e*new_e-vp_m*vp_m;
00199         new_p=sqrt(new_p);
00200         double new_pt=new_p/cosh(new_tlv.Eta());
00201 
00202         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00203         
00204       }
00205       //jet
00206       if (id_pdg!=15&&id_pdg!=13&&id_pdg!=11&&id_pdg!=22){
00207 
00208         //vp_m=0!!
00209 
00210         /* Note that E is in GeV in jet_resolution */
00211         double e_sys=new_tlv.E();
00212         double eta_sys=new_tlv.Eta();
00213         double new_e=jet_resolution(mev_adapt*e_sys,eta_sys)/mev_adapt;
00214         double new_p=new_e*new_e-vp_m*vp_m;
00215         new_p=sqrt(new_p);
00216         double new_pt=new_p/cosh(new_tlv.Eta());
00217 
00218         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00219         
00220       }
00221 
00222     }
00223 
00224     /*---------------------------------------------*/
00225 
00226     if( f_eScale_minus==1 ){
00227       /* --------------------------------------
00228          'negativ' systematics in e scale
00229          -------------------------------------- */
00230       
00231       //electron
00232       if (id_pdg==11){
00233 
00234         double et_sys=new_tlv.Et();
00235         double new_et=energy_scale(et_sys,egamma_scale_factor_minus);
00236         double new_e=new_et*cosh(new_tlv.Eta());
00237         double new_p=new_e*new_e-vp_m*vp_m;
00238         new_p=sqrt(new_p);
00239         double new_pt=new_p/cosh(new_tlv.Eta());
00240 
00241         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00242         
00243       }//photon
00244       if (id_pdg==22){
00245 
00246         double et_sys=new_tlv.Et();
00247         double new_et=energy_scale(et_sys,egamma_scale_factor_minus);
00248         double new_e=new_et*cosh(new_tlv.Eta());
00249         double new_p=new_e*new_e-vp_m*vp_m;
00250         new_p=sqrt(new_p);
00251         double new_pt=new_p/cosh(new_tlv.Eta());
00252 
00253         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00254         
00255       }
00256       //muon
00257       if (id_pdg==13){
00258 
00259         double pt_sys=new_tlv.Pt();
00260         double new_pt=energy_scale(pt_sys,muon_scale_factor_minus);;
00261         double new_p=new_pt*cosh(new_tlv.Eta());
00262         double new_e=new_p*new_p+vp_m*vp_m;
00263         new_e=sqrt(new_e);
00264 
00265         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00266         
00267       }
00268       //tau-jet
00269       if (id_pdg==15){
00270 
00271         double e_sys=new_tlv.E();
00272         double new_e=energy_scale(e_sys,taujet_scale_factor_minus);;
00273         double new_p=new_e*new_e-vp_m*vp_m;
00274         new_p=sqrt(new_p);
00275         double new_pt=new_p/cosh(new_tlv.Eta());
00276 
00277         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00278         
00279       }
00280       //jet
00281       if (id_pdg!=15&&id_pdg!=13&&id_pdg!=11&&id_pdg!=22){
00282 
00283         //vp_m=0!!
00284 
00285         double e_sys=new_tlv.E();
00286         double eta_sys=new_tlv.Eta();
00287         double new_e=energy_scale(e_sys,jet_scale_factor_minus);
00288         if(TMath::Abs(eta_sys)>3.2) new_e=energy_scale(e_sys,jet_scale_factor_minus_large_eta);
00289         double new_p=new_e*new_e-vp_m*vp_m;
00290         new_p=sqrt(new_p);
00291         double new_pt=new_p/cosh(new_tlv.Eta());
00292 
00293         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00294         
00295       }
00296 
00297     }
00298 
00299     /*---------------------------------------------*/
00300 
00301     if( f_eScale_plus==1 ){
00302       /* --------------------------------------
00303          'negativ' systematics in e scale
00304          -------------------------------------- */
00305       
00306       //electron
00307       if (id_pdg==11){
00308 
00309         double et_sys=new_tlv.Et();
00310         double new_et=energy_scale(et_sys,egamma_scale_factor_plus);
00311         double new_e=new_et*cosh(new_tlv.Eta());
00312         double new_p=new_e*new_e-vp_m*vp_m;
00313         new_p=sqrt(new_p);
00314         double new_pt=new_p/cosh(new_tlv.Eta());
00315 
00316         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00317         
00318       }
00319       //photon
00320       if(id_pdg==22){
00321 
00322         double et_sys=new_tlv.Et();
00323         double new_et=energy_scale(et_sys,egamma_scale_factor_plus);
00324         double new_e=new_et*cosh(new_tlv.Eta());
00325         double new_p=new_e*new_e-vp_m*vp_m;
00326         new_p=sqrt(new_p);
00327         double new_pt=new_p/cosh(new_tlv.Eta());
00328 
00329         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00330         
00331       }
00332       //muon
00333       if (id_pdg==13){
00334 
00335         double pt_sys=new_tlv.Pt();
00336         double new_pt=energy_scale(pt_sys,muon_scale_factor_plus);;
00337         double new_p=new_pt*cosh(new_tlv.Eta());
00338         double new_e=new_p*new_p+vp_m*vp_m;
00339         new_e=sqrt(new_e);
00340 
00341         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00342         
00343       }
00344       //tau-jet
00345       if (id_pdg==15){
00346 
00347         double e_sys=new_tlv.E();
00348         double new_e=energy_scale(e_sys,taujet_scale_factor_plus);;
00349         double new_p=new_e*new_e-vp_m*vp_m;
00350         new_p=sqrt(new_p);
00351         double new_pt=new_p/cosh(new_tlv.Eta());
00352 
00353         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00354         
00355       }
00356       //jet
00357       if (id_pdg!=15&&id_pdg!=13&&id_pdg!=11&&id_pdg!=22){
00358 
00359         //vp_m=0!!
00360 
00361         double e_sys=new_tlv.E();
00362         double eta_sys=new_tlv.Eta();
00363         double new_e=energy_scale(e_sys,jet_scale_factor_plus);
00364         if(TMath::Abs(eta_sys)>3.2) new_e=energy_scale(e_sys,jet_scale_factor_plus_large_eta);
00365         double new_p=new_e*new_e-vp_m*vp_m;
00366         new_p=sqrt(new_p);
00367         double new_pt=new_p/cosh(new_tlv.Eta());
00368 
00369         new_tlv.SetPtEtaPhiE(new_pt,new_tlv.Eta(),new_tlv.Phi(),new_e);
00370         
00371       }
00372 
00373     }
00374 
00375     /*---------------------------------------------*/
00376 
00377     if( f_effi==1 ){
00378       /* --------------------------------------
00379          systematics on effi
00380          -------------------------------------- */
00381       
00382       //electron
00383       if (id_pdg==11){
00384         syst_remove=!is_detected(electron_eff_error);
00385       }
00386       //photon
00387       else if (id_pdg==22){
00388         syst_remove=!is_detected(photon_eff_error);
00389       }
00390       //muon
00391       else if (id_pdg==13){
00392         syst_remove=!is_detected(muon_eff_error);
00393       }
00394       //tau jet
00395       else if (id_pdg==15){
00396         syst_remove=!is_detected(tau_eff_error);
00397       }
00398       //jets (b and light)
00399       else {
00400         double eta_sys=new_tlv.Eta();
00401         if (TMath::Abs(eta_sys) > 3.2 ) {
00402           syst_remove=!is_detected(jet_eff_error_large_eta);
00403         } else {
00404           syst_remove=!is_detected(jet_eff_error);
00405         }
00406       }
00407     }
00408     /*---------------------------------------------*/
00409 
00410     if( f_btag==1 && (!syst_remove) ){
00411       /* --------------------------------------
00412          systematics on b tagging
00413              the user himself must make sure, that only b, not light, jets are taken!!
00414          -------------------------------------- */
00415       //b jet
00416       if (id_pdg!=0){ cout<<"** WARNING MyTools::apply_syst: b tagging effi on e, photon, mu or tau \n"; }
00417 
00418       //      bjet_eff_error=bjet_eff_error_scale*bjet_eff; //should be in pt/eta bins
00419       bjet_eff_error=bjet_eff_error_scale; //changed 26.2.08 (K.A. also changed)
00420       
00421       syst_remove=!is_bjet(bjet_eff_error);
00422 
00423     }
00424     
00425     (*vp_syst)[il]->set_tlv(new_tlv);
00426 
00427 //     cout<<"part "<<il<<" id="<<id_pdg<<" syst_remove = "<<syst_remove<<"\n";
00428 //     cout<<"              old pt  = "<<aod_tlv.Pt()<<" new pt = "<<new_tlv.Pt()<<" pt = "<<(*vp_syst)[il]->tlv().Pt()<<"\n";
00429 //     cout<<"              old phi = "<<aod_tlv.Phi()<<" new phi = "<<new_tlv.Phi()<<" ph = "<<(*vp_syst)[il]->tlv().Phi()<<"\n";
00430 //     cout<<"              old eta = "<<aod_tlv.Eta()<<" new eta = "<<new_tlv.Eta()<<" eta = "<<(*vp_syst)[il]->tlv().Eta()<<"\n";
00431 //     cout<<"              old e = "<<aod_tlv.E()<<" new e = "<<new_tlv.E()<<" e = "<<(*vp_syst)[il]->tlv().E()<<"\n";
00432 
00433     if (!syst_remove) {  
00434 
00435       // recalculated etsimate of MET due to scale and reso syst 
00436       // if ONLY f_effi set are aod_px and new_px the same => no change
00437 
00438       double new_px=new_tlv.Px();
00439       double new_py=new_tlv.Py();
00440       double new_pt=new_tlv.Pt();
00441       double aod_px=aod_tlv.Px();
00442       double aod_py=aod_tlv.Py();
00443       double aod_pt=aod_tlv.Pt();
00444 
00445       //the following is ONLY correct, if the effects are investigated seperately!      
00446       if ((f_eScale_plus==1)&&!((id_pdg==11)||(id_pdg==13)||(id_pdg==22))){
00447 
00448         if (id_pdg==15){
00449           new_metx += metCorrection_x(0.0,taujet_METscale_factor_plus*aod_px);
00450           // new_metx +=metScale ( aod_px, 0.0, taujet_METscale_factor_plus) //(K.A. has -=)
00451           new_mety += metCorrection_y(0.0,taujet_METscale_factor_plus*aod_py);
00452           new_metsum += metCorrection_sum(0.0,taujet_METscale_factor_plus*aod_pt);
00453         }else{
00454           new_metx += metCorrection_x(0.0,jet_METscale_factor_plus*aod_px);
00455           new_mety += metCorrection_y(0.0,jet_METscale_factor_plus*aod_py);
00456           new_metsum += metCorrection_sum(0.0,jet_METscale_factor_plus*aod_pt);
00457         }
00458 
00459       }else if((f_eScale_minus==1)&&!((id_pdg==11)||(id_pdg==13)||(id_pdg==22))){
00460 
00461         if (id_pdg==15){
00462           new_metx += metCorrection_x(0.0,taujet_METscale_factor_minus*aod_px);
00463           new_mety += metCorrection_y(0.0,taujet_METscale_factor_minus*aod_py);
00464           new_metsum += metCorrection_sum(0.0,taujet_METscale_factor_minus*aod_pt);
00465         }else{
00466           new_metx += metCorrection_x(0.0,jet_METscale_factor_minus*aod_px);
00467           new_mety += metCorrection_y(0.0,jet_METscale_factor_minus*aod_py);
00468           new_metsum += metCorrection_sum(0.0,jet_METscale_factor_minus*aod_pt);
00469         }
00470 
00471       }else{
00472         
00473         new_metx += metCorrection_x(aod_px,new_px);
00474         new_mety += metCorrection_y(aod_py,new_py);
00475         new_metsum += metCorrection_sum(aod_pt,new_pt);
00476 
00477       }
00478       new_met=sqrt( new_metx*new_metx + new_mety*new_mety );
00479 
00480     } else {
00481 
00482       // delete syst_remove => should not change MET  (exept muons)
00483 
00484       if (id_pdg==13){
00485         double aod_px=aod_tlv.Px();
00486         double aod_py=aod_tlv.Py();
00487         double aod_pt=aod_tlv.Pt();
00488         
00489         new_metx += metCorrection_x(aod_px,0.0);
00490         new_mety += metCorrection_y(aod_py,0.0);
00491         new_metsum += metCorrection_sum(aod_pt,0.0);
00492 
00493         new_met=sqrt( new_metx*new_metx + new_mety*new_mety );
00494       }
00495 
00496      
00497       // cout<<"removing particle at position: "<<il<<endl;
00498       vp_syst->erase(vp_syst->begin()+il);
00499 
00500     }
00501 
00502   }// vp loop
00503 
00504 //  cout<<" MET:  old = "<<met_syst->met()<<" new = "<<new_met<<"\n";
00505 
00506   met_syst->set_met_x(new_metx);
00507   met_syst->set_met_y(new_mety);
00508   met_syst->set_met_sum(new_metsum);
00509 
00510 //  cout<<" MET:  et = "<<met_syst->met()<<"\n";
00511 //  cout<<"after erase:"<<vp_syst->size()<<"\n";
00512   
00513 /*  sort particle vector by pt */
00514   *vp_syst=m_tools.sort_particles(*vp_syst); 
00515 
00516 //   for (int il=0;il<vp_syst->size();il++){
00517 //     cout<<il<<" pt="<<(*vp_syst)[il]->tlv().Pt()<<"\n";
00518 //   }
00519 
00520 }
00521 
00522 //*****************************************************************************