Reference Manual

Go to the documentation of this file.
00001 
00013 #include "TopQuarkAnalysis/TopPairBSM/interface/JetCombinatorics.h"
00014 #include "TMath.h"
00015 
00016 #include <iostream>
00017 
00018 std::string itoa(int i) {
00019   char temp[20];
00020   sprintf(temp,"%d",i);
00021   return((std::string)temp);
00022 }
00023 
00024 
00025 
00026 //_______________________________________________________________
00027 JetCombinatorics::JetCombinatorics() {
00028 
00029         this->Clear();
00030                 
00031         minMassLepW_ = -999999.;
00032         maxMassLepW_ = 999999.;
00033         minMassHadW_ = -999999.;
00034         maxMassHadW_ = 999999.;
00035         minMassLepTop_ = -999999.;
00036         maxMassLepTop_ = 999999.;
00037         
00038         minPhi_ = -1.;
00039         removeDuplicates_ = true;
00040         maxNJets_ = 9999;
00041         verbosef = false;
00042         UsebTagging_ = false;
00043         UseMtop_ = true;
00044         SigmasTypef = 0;
00045         UseFlv_ = false;
00046         
00047         Template4jCombos_ = NestedCombinatorics(); // 12 combinations
00048         Template5jCombos_ = Combinatorics(4,5); // 5 combinations of 4 combos
00049         Template6jCombos_ = Combinatorics(4,6); // 15 combinations of 4 combos
00050         Template7jCombos_ = Combinatorics(4,7); // xx combinations of 4 combos
00051         
00052 }
00053 
00054 //_______________________________________________________________
00055 JetCombinatorics::~JetCombinatorics() {
00056   this->Clear();
00057 }
00058 
00059 //_______________________________________________________________
00060 void JetCombinatorics::Clear() {
00061 
00062         allCombos_.clear();
00063         allCombosSumEt_.clear();
00064         Template4jCombos_.clear();
00065         Template5jCombos_.clear();
00066         Template6jCombos_.clear();
00067         Template7jCombos_.clear();
00068         cand1_.clear();
00069         
00070 }
00071         
00072 
00073 //_______________________________________________________________
00074 std::map< int, std::string > JetCombinatorics::Combinatorics(int n, int max) {
00075 
00076         // find a combinatorics template
00077         // This is a simple stupid function to make algebratic combinatorics
00078 
00079         int kcombos = n;
00080         int maxcombos = max;
00081 
00082         std::string list;
00083 
00084         for ( int m=0; m<maxcombos; m++) { list = list + (itoa(m));}
00085 
00086         std::string seed;
00087         for ( int m=0; m<kcombos; m++) { seed = seed + (itoa(m));}
00088 
00089         
00090         std::map< int, std::string > aTemplateCombos;
00091         aTemplateCombos.clear();
00092         
00093         aTemplateCombos[0] = seed;
00094 
00095         int i = 0;
00096         int totalmatches = seed.size();
00097         int totalIte = list.size();
00098 
00099         for ( int ite = 0; ite < ((int)totalIte); ite++) {
00100 
00101                 //cout << "iteration " << ite << endl;
00102                 //i = 0;
00103                 //for ( Itevec = seed.begin(); Itevec != seed.end(); ++Itevec) {
00104                 for ( i=0; i< (int) totalmatches; i++) {
00105 
00106                         std::string newseed = aTemplateCombos[ite];
00107                         std::string newseed2;
00108                         /*
00109                         cout << "newseed size= " << newseed.size() << " : ";
00110                         for (std::vector< std::string>::iterator iite = newseed.begin();
00111                                  iite != newseed.end(); ++iite) {
00112 
00113                                 cout << *iite << " ";
00114                         }
00115                         cout << endl;
00116                         */
00117                         for ( int itemp=0; itemp<(int)newseed.size(); itemp++) {
00118                                 if (itemp!=i) newseed2 = newseed2 + (newseed[itemp]);
00119                         }
00120                         /*
00121                         cout << "newseed2: ";
00122                         for (std::vector< std::string>::iterator iite = newseed2.begin();
00123                                  iite != newseed2.end(); ++iite) {
00124 
00125                                 cout << *iite << " ";
00126                         }
00127                         cout << endl;
00128                         */
00129                         for ( int j=0; j<(int) list.size(); j++) {
00130                                 //cout << " j = " << j << endl;
00131                                 bool Isnewelement = true;
00132                                 std::string newelement = "0";
00133                                 //bool Isnewcombo = true;
00134                                 for (int k=0; k< (int)newseed2.size(); k++) {
00135                                         if ( list[j] == newseed2[k] ) Isnewelement = false;
00136                                 }
00137                                 if (Isnewelement) {
00138 
00139                                         newelement = list[j];
00140                                         //cout << "new element: " << newelement << endl;
00141 
00142                                         std::string candseed = newseed2;
00143                                         candseed = candseed + newelement;
00144 
00145                                         bool IsnewCombo = true;
00146                                         for (int ic=0; ic<(int)aTemplateCombos.size(); ++ic ) {
00147 
00148                                                 int nmatch = 0;
00149                                                 for ( int ij=0; ij<(int)(aTemplateCombos[ic]).size(); ij++) {
00150 
00151                                                         for (int ik=0; ik<(int)candseed.size(); ik++) {
00152                                                                 if ( candseed[ik] == aTemplateCombos[ic][ij] ) nmatch++;
00153                                                         }
00154                                                 }
00155                                                 if (nmatch == (int)totalmatches)
00156                                                         IsnewCombo = false;
00157 
00158                                         }
00159                                         if (IsnewCombo) {
00160                                                 //cout << "new combo"<< " before combo size=" << aTemplateCombos.size() << endl;
00161                                                 aTemplateCombos[(int)aTemplateCombos.size()] = candseed;
00162                                                 //cout << " after size = " << aTemplateCombos.size() << endl;
00163                                         }
00164                                 }
00165 
00166                         }
00167                 }
00168         }//close iterations
00169 
00170         // debug info
00171         
00172         //std::cout << " build templates for total combos = " << aTemplateCombos.size() << std::endl;
00173         //std::cout << " template combos: " << std::endl;
00174         //for (size_t ic=0; ic != aTemplateCombos.size(); ++ic) {
00175 
00176         //std::cout << aTemplateCombos[ic] << std::endl;
00177         //}
00178         
00179         return aTemplateCombos;
00180         
00181         
00182         
00183 }
00184 
00185 
00186 //______________________________________________________________
00187 std::map< int, std::string > JetCombinatorics::NestedCombinatorics() {
00188 
00189         // build by hand 12 combinations for semileptonic top decays
00190 
00191         std::map< int, std::string > aTemplateCombos;
00192         aTemplateCombos.clear();
00193         
00194         aTemplateCombos[0] = "0123";
00195         aTemplateCombos[1] = "0132";
00196         aTemplateCombos[2] = "0213";
00197         aTemplateCombos[3] = "0231";
00198         aTemplateCombos[4] = "0312";
00199         aTemplateCombos[5] = "0321";
00200         aTemplateCombos[6] = "1203";
00201         aTemplateCombos[7] = "1230";
00202         aTemplateCombos[8] = "1302";
00203         aTemplateCombos[9] = "1320";
00204         aTemplateCombos[10] = "2301";
00205         aTemplateCombos[11] = "2310";
00206                 
00207         return aTemplateCombos;
00208 
00209 }
00210 
00211 //______________________________________________________________
00212 void JetCombinatorics::FourJetsCombinations(std::vector<TLorentzVector> jets, std::vector<double> bdiscriminators ) {
00213 
00214 
00215         int n = 0; // total number of combos
00216         std::map< Combo, int, minChi2 > allCombos;
00217         std::map< Combo, int, maxSumEt > allCombosSumEt;
00218         
00219         std::map< int, std::string > aTemplateCombos;
00220         aTemplateCombos.clear();
00221 
00222         if ( jets.size() == 4 ) aTemplateCombos[0] = std::string("0123");
00223         if ( jets.size() == 5 ) aTemplateCombos = Template5jCombos_;
00224         if ( jets.size() == 6 ) aTemplateCombos = Template6jCombos_;
00225         if ( jets.size() == 7 ) aTemplateCombos = Template7jCombos_;    
00226 
00227         // force to use only 4 jets
00228         if ( maxNJets_ == 4 ) aTemplateCombos[0] = std::string("0123");
00229         
00230         if (verbosef) std::cout << "[JetCombinatorics] size of vector of jets = " << jets.size() << std::endl;
00231         
00232         for (size_t ic=0; ic != aTemplateCombos.size(); ++ic) {
00233 
00234                 if (verbosef) std::cout << "[JetCombinatorics] get 4 jets from the list, cluster # " << ic << "/"<< aTemplateCombos.size()-1 << std::endl;
00235                 
00236                 // get a template
00237                 std::string aTemplate = aTemplateCombos[ic];
00238 
00239                 if (verbosef) std::cout << "[JetCombinatorics] template of 4 jets = " << aTemplate << std::endl;
00240                 
00241                 // make a list of 4 jets
00242                 std::vector< TLorentzVector > the4jets;
00243                 std::vector< int > the4Ids;
00244                 std::vector< double > thebdisc;
00245                 std::vector< double > theFlvCorr;
00246                 //the4jets[0] = jets[0];
00247                 
00248                 for (int ij=0; ij<4; ij++) {
00249                         //std::cout << "ij= " << ij << std::endl;
00250                         //std::cout << "atoi = " << atoi((aTemplate.substr(0,1)).c_str()) << std::endl;
00251                         //std::cout << "jets[].Pt = " << jets[ij].Pt() << std::endl;
00252                         int tmpi = atoi((aTemplate.substr(ij,1)).c_str());
00253                         //std::cout << "tmpi= " << tmpi << std::endl;
00254                         the4jets.push_back(jets[tmpi]);
00255                         the4Ids.push_back(tmpi);
00256                         if ( UsebTagging_ ) thebdisc.push_back( bdiscriminators[tmpi] );
00257                         if ( UseFlv_ ) theFlvCorr.push_back( flavorCorrections_[tmpi] );
00258                 }
00259 
00260                 if (verbosef) std::cout<< "[JetCombinatorics] with these 4 jets, make 12 combinations: " <<std::endl;
00261 
00262                 //std::cout << " the4jets[ij].size = " << the4jets.size() << std::endl;
00263                         
00264                 for (size_t itemplate=0; itemplate!= Template4jCombos_.size(); ++itemplate) {
00265                         
00266                         std::string a4template = Template4jCombos_[itemplate];
00267 
00268                         if (verbosef) std::cout << "[JetCombinatorics] ==> combination: " << a4template << " is # " << itemplate << "/"<<  Template4jCombos_.size()-1 << std::endl;
00269                         
00270                         Combo acombo;
00271                         
00272                         acombo.SetWp( the4jets[atoi((a4template.substr(0,1)).c_str())] );
00273                         acombo.SetWq( the4jets[atoi((a4template.substr(1,1)).c_str())] );
00274                         acombo.SetHadb( the4jets[atoi((a4template.substr(2,1)).c_str())] );
00275                         acombo.SetLepb( the4jets[atoi((a4template.substr(3,1)).c_str())] );
00276                         acombo.SetLepW( theLepW_ );
00277 
00278                         acombo.SetIdWp( the4Ids[atoi((a4template.substr(0,1)).c_str())] );
00279                         acombo.SetIdWq( the4Ids[atoi((a4template.substr(1,1)).c_str())] );
00280                         acombo.SetIdHadb( the4Ids[atoi((a4template.substr(2,1)).c_str())] );
00281                         acombo.SetIdLepb( the4Ids[atoi((a4template.substr(3,1)).c_str())] );
00282                         //std::cout << " acombo setup" << std::endl;
00283 
00284                         if ( UseFlv_ ) {
00285                                 acombo.SetFlvCorrWp( theFlvCorr[atoi((a4template.substr(0,1)).c_str())] );
00286                                 acombo.SetFlvCorrWq( theFlvCorr[atoi((a4template.substr(1,1)).c_str())] );
00287                                 acombo.SetFlvCorrHadb( theFlvCorr[atoi((a4template.substr(2,1)).c_str())] );
00288                                 acombo.SetFlvCorrLepb( theFlvCorr[atoi((a4template.substr(3,1)).c_str())] );
00289                                 acombo.ApplyFlavorCorrections();
00290                         }
00291                         if ( UsebTagging_ ) {
00292 
00293                                 acombo.Usebtagging();
00294                                 acombo.SetbDiscPdf(bTagPdffilename_);
00295                                 acombo.SetWp_disc( thebdisc[atoi((a4template.substr(0,1)).c_str())] );
00296                                 acombo.SetWq_disc( thebdisc[atoi((a4template.substr(1,1)).c_str())] );
00297                                 acombo.SetHadb_disc( thebdisc[atoi((a4template.substr(2,1)).c_str())] );
00298                                 acombo.SetLepb_disc( thebdisc[atoi((a4template.substr(3,1)).c_str())] );
00299                                 
00300                         }
00301 
00302                         acombo.UseMtopConstraint(UseMtop_);
00303                         // choose value of sigmas
00304                         acombo.SetSigmas(SigmasTypef);
00305 
00306                         acombo.analyze();
00307 
00308                         if (verbosef) {
00309 
00310                           std::cout << "[JetCombinatorics] ==> combination done:" << std::endl;
00311                           acombo.Print();
00312                         }
00313 
00314                         // invariant mass cuts
00315                         TLorentzVector aHadWP4 = acombo.GetHadW();
00316                         TLorentzVector aLepWP4 = acombo.GetLepW();
00317                         TLorentzVector aLepTopP4=acombo.GetLepTop();
00318                         
00319                         if ( ( aHadWP4.M() > minMassHadW_ && aHadWP4.M() < maxMassHadW_ ) &&
00320                                  ( aLepWP4.M() > minMassLepW_ && aLepWP4.M() < maxMassLepW_ ) &&
00321                                  ( aLepTopP4.M() > minMassLepTop_ && aLepTopP4.M() < maxMassLepTop_) ) {
00322                         
00323                                 allCombos[acombo] = n;
00324                                 allCombosSumEt[acombo] = n;
00325                         
00326                                 n++;
00327                         }
00328                 
00329                 }
00330         }
00331 
00332         allCombos_ = allCombos;
00333         allCombosSumEt_ = allCombosSumEt;
00334        
00335 }
00336 
00337 Combo JetCombinatorics::GetCombination(int n) {
00338 
00339         int j = 0;
00340         Combo a;
00341         for ( std::map<Combo,int,minChi2>::const_iterator ite=allCombos_.begin();
00342                   ite!=allCombos_.end(); ++ite) {
00343                 
00344                 if (j == n) a = ite->first;
00345                 j++;
00346         }
00347 
00348         return a;
00349 
00350         
00351 }
00352 
00353 Combo JetCombinatorics::GetCombinationSumEt(int n) {
00354 
00355         int j = 0;
00356         Combo a;
00357         for ( std::map<Combo,int,maxSumEt>::const_iterator ite=allCombosSumEt_.begin();
00358                   ite!=allCombosSumEt_.end(); ++ite) {
00359                 
00360                 if (j == n) a = ite->first;
00361                 j++;
00362         }
00363 
00364         return a;
00365 
00366         
00367 }
/data/refman/pasoursint/CMSSW_5_3_10_patch2/src/TopQuarkAnalysis/TopPairBSM/src/JetCombinatorics.cc
