ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Angles/src/fillBdtTree.cc
(Generate patch)

Comparing UserCode/MitHzz4l/Angles/src/fillBdtTree.cc (file contents):
Revision 1.3 by dkralph, Mon Dec 17 12:35:15 2012 UTC vs.
Revision 1.4 by andrey, Tue Dec 18 00:25:34 2012 UTC

# Line 1 | Line 1
1   #include <iostream>
2   #include <assert.h>
3 < #include <stdio.h>    
4 < #include <stdlib.h>    
3 > #include <stdio.h>
4 > #include <stdlib.h>
5   #include <getopt.h>
6   #include <string.h>
7   #include <utility>
# Line 18 | Line 18
18   #include "ZZMatrixElement/MELA/interface/Mela.h"
19   #include "ZZMatrixElement/MELA/interface/PseudoMELA.h"
20  
21 + #include "ZZMatrixElement/MELA/src/computeAngles.h"
22 +
23 +
24   #include "PwhgWrapper.h"
25  
26   #include "TMVA/Reader.h"
# Line 43 | Line 46 | using namespace std;
46  
47   pwhegwrapper *wrap;
48  
49 < class CtrlFlags
50 < {
49 > class CtrlFlags
50 > {
51   public :
52    CtrlFlags():inputfile(""),
53                inputtree("zznt"),
# Line 85 | Line 88 | public :
88   };
89   void CtrlFlags::dump()
90   {
91 <  cout << "mH           : " << mH << endl;
92 <  cout << "inputfile    : " << inputfile << endl;
93 <  cout << "inputtree    : " << inputtree << endl;
94 <  cout << "weightfile   : " << weightfile << endl;
95 <  cout << "output       : " << output << endl;
96 <  cout << "mH_lo        : " << mH_lo << endl;
97 <  cout << "mH_hi        : " << mH_hi << endl;
98 <  cout << "wL           : " << wL << endl;
99 <  cout << "wH           : " << wH << endl;
100 <  cout << "varstring    : " << varstring << endl;
101 <  cout << "makeMela     : " << makeMela << endl;
102 <  cout << "multiclass   : " << multiclass << endl;
103 <  cout << "flattree     : " << flattree << endl;
104 <  cout << "multisigs    : " << multisigs << endl;
105 <  cout << "withPtY      : " << withPtY << endl;
106 <  cout << "hwidth       : " << hwidth << endl;
91 >  cout << "mH           : " << mH << endl;
92 >  cout << "inputfile    : " << inputfile << endl;
93 >  cout << "inputtree    : " << inputtree << endl;
94 >  cout << "weightfile   : " << weightfile << endl;
95 >  cout << "output       : " << output << endl;
96 >  cout << "mH_lo        : " << mH_lo << endl;
97 >  cout << "mH_hi        : " << mH_hi << endl;
98 >  cout << "wL           : " << wL << endl;
99 >  cout << "wH           : " << wH << endl;
100 >  cout << "varstring    : " << varstring << endl;
101 >  cout << "makeMela     : " << makeMela << endl;
102 >  cout << "multiclass   : " << multiclass << endl;
103 >  cout << "flattree     : " << flattree << endl;
104 >  cout << "multisigs    : " << multisigs << endl;
105 >  cout << "withPtY      : " << withPtY << endl;
106 >  cout << "hwidth       : " << hwidth << endl;
107    assert(era==2011 || era==2012);
108 <  cout << "era          : " << era << endl;
108 >  cout << "era          : " << era << endl;
109   }
110   PseudoMELA *psMela;
111   void parse_args( int, char**, CtrlFlags & );
112   void setFracs(float mH, float mH_lo, float mH_hi, float &wgt_lo, float &wgt_hi);
113   void setZzntVals(CtrlFlags &flags, varkeepter *flatoutput, filestuff *fs, TString cls, TMVA::Reader *reader=0, Mela *mela=0, TMVA::Reader *spinReader=0);
114 +
115 + //**AP: added this:
116 + void setInterferenceWeight(Bool_t isHiggs, varkeepter *flatoutput, filestuff *fs,  Mela *mela=0);
117 + Bool_t isHiggs = kTRUE;
118 + // ***
119 +
120   double getWeight(TH1F *hist, float val);
121   void setPtWgtHists(TString cls, TString fname, Long_t mH, TFile *wgtFile);
122   //----------------------------------------------------------------------------------------
123   map<TString,TH1F*> wgtHists;
124  
125   //------------------------------------------------------------
126 < int main(int argc, char ** argv)  
126 > int main(int argc, char ** argv)
127   //------------------------------------------------------------
128   {
129    CtrlFlags flags;
# Line 125 | Line 134 | int main(int argc, char ** argv)
134  
135    TFile *pt4lWgts = TFile::Open("data/pt4l/pt4l_weights.root");
136    wrap = new pwhegwrapper;
137 <  
137 >
138    TString cmsswpath(CMSSW_BASE + TString("/src"));
139    string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
140    SimpleTable xstab(xspath.c_str());
# Line 153 | Line 162 | int main(int argc, char ** argv)
162        TString dummy;
163        // if(TString(line).Contains("^ntupledir")) ss >> dummy >> ntupledir;
164        if(TString(line).Contains("^spinWeightFile")) { ss >> dummy >> spinWeightFile; cout << "spinWeightFile: " << spinWeightFile << endl; }
165 <      continue;
165 >      continue;
166      }
167  
168      TString cls,fname,type;
# Line 172 | Line 181 | int main(int argc, char ** argv)
181    }
182    ifs.close();
183  
184 <  // mass cut based on weightfile name ...
184 >  // mass cut based on weightfile name ...
185    float mH = flags.mH;
186    // boost::regex expression2(".*_mH(.*)_.*");
187    // boost::cmatch match2;
# Line 182 | Line 191 | int main(int argc, char ** argv)
191    // assert(flags.wL>=0 && flags.wL < 1000 && flags.wH >= 0 && flags.wH < 1000);
192    char cutstr[512];
193    sprintf(cutstr, "m4l>%.3f&&m4l<%.3f", flags.wL, flags.wH);
194 <  cout << "mH: "   << setw(8) << mH
195 <       << "  wL: " << setw(8) << flags.wL
196 <       << "  wH: " << setw(8) << flags.wH
197 <       << "  mH_lo: " << setw(8) << flags.mH_lo
198 <       << "  mH_hi: " << setw(8) << flags.mH_hi
194 >  cout << "mH: "   << setw(8) << mH
195 >       << "  wL: " << setw(8) << flags.wL
196 >       << "  wH: " << setw(8) << flags.wH
197 >       << "  mH_lo: " << setw(8) << flags.mH_lo
198 >       << "  mH_hi: " << setw(8) << flags.mH_hi
199         << endl;
200  
201    float frac_lo=1,frac_hi=1;
# Line 258 | Line 267 | int main(int argc, char ** argv)
267    // setup output
268    //
269    filestuff fs(tmpfname,tmpfname,"",classes[0]=="data");
270 <  float f_bdt;
270 >  float f_bdt;
271    TBranch *br_bdt=0;
272    varkeepter *flatoutput=0;
273    TString varfname("zzntvars");
# Line 341 | Line 350 | int main(int argc, char ** argv)
350    int ifile=-1;
351    int nent=fs.get_zz_chain()->GetEntries();
352    cout << "beginning tree from mass " << mH << "(" << flags.wL << "," << flags.wH << ") with " << nent << "  entries" << endl;
353 <  for( int i=0; i<nent; i++ ) {
353 >  for( int i=0; i<nent; i++ ) {
354      if( !(i%5000) ) cerr << "processed:  " << i << "/" << nent << endl;
355      fs.getentry(i,"","zznt");
356  
# Line 386 | Line 395 | int main(int argc, char ** argv)
395        assert(0);
396        fs.outtuple->Fill();
397      } else {
398 +
399 +      //AP: Here, right?:
400 +      setInterferenceWeights(isHiggs, flatoutput, &fs, mela);
401 +      // **
402 +
403        setZzntVals(flags, flatoutput, &fs, classes[ifile], reader, mela, spinReader);
404        flatoutput->fill();
405      }
# Line 400 | Line 414 | int main(int argc, char ** argv)
414    }
415   }
416   //----------------------------------------------------------
417 < void parse_args( int argc, char** argv, CtrlFlags &flags )
417 > void parse_args( int argc, char** argv, CtrlFlags &flags )
418   //----------------------------------------------------------
419   {
420    int c;
# Line 429 | Line 443 | void parse_args( int argc, char** argv,
443        {"hwidth", 1, 0, 'q'},
444        {0, 0, 0, 0}
445      };
446 <    
446 >
447      c = getopt_long (argc, argv, "a:e",
448                       long_options, &option_index);
449      if (c == -1)
450        break;
451 <    
451 >
452      switch (c) {
453      case 0:
454        printf ("option %s", long_options[option_index].name);
# Line 442 | Line 456 | void parse_args( int argc, char** argv,
456          printf (" with arg %s", optarg);
457        printf ("\n");
458        break;
459 <      
459 >
460      case 'a':
461        flags.inputfile = std::string(optarg);
462        break;
# Line 501 | Line 515 | void parse_args( int argc, char** argv,
515        digit_optind = this_option_optind;
516        printf ("option %c\n", c);
517        break;
518 <      
518 >
519      case '?':
520        break;
521 <      
521 >
522      default:
523        printf ("?? getopt returned character code 0%o ??\n", c);
524      }
525    }
526 <  
526 >
527    if( flags.inputfile.empty() ) {
528      cerr << "please specify an inputfile with --inputfile " << endl;
529      exit(1);
# Line 601 | Line 615 | void setZzntVals(CtrlFlags &flags, varke
615      float costhetastar,costheta1,costheta2,phi,phistar1,kd,psig,pbkg,spinKd;
616      mela->computeKD(Z1_lept1, Z1_lept1Id, Z1_lept2, Z1_lept2Id, Z2_lept1, Z2_lept1Id, Z2_lept2, Z2_lept2Id,
617                      costhetastar, costheta1, costheta2, phi, phistar1, kd, psig, pbkg, flags.withPtY, flags.withPtY);
618 <  
618 >
619      flatoutput->setval("melaLD",                (Float_t)kd);
620  
621      if(flags.multisigs) {
# Line 700 | Line 714 | void setPtWgtHists(TString cls, TString
714      findHist(cls, mH, fname, wgtFile, "hres",    "Hi");
715    }
716   }
717 +
718 +
719 +
720 +
721 +
722 + /* AP ///  Adding interference weights using Mela's functionality.
723 + I don't know if it's actually work in here, couldn't compile it.
724 + Note: in order for this to work, one needs:
725 + 1) Mela's compteAngles.h and Mela.h
726 + 2) mcfm library linked
727 + 3) There maight be a problem with name 'mela'.
728 +  -- computeAngles() is defined in a namespace 'mela::'
729 +  -- The object of class Mela is also called 'mela'. If there is a problem - I'd rename it to 'mela1'
730 + */
731 +
732 + void setInterferenceWeight(Bool_t isHiggs, varkeepter *flatoutput, filestuff *fs,  Mela *mela1)
733 + {
734 +
735 +  float   w_interference = 1;
736 +
737 +  if (!isHiggs)
738 +    {
739 +      // If the sample we're running is not higgs - set this weight to 1 and exit
740 +      flatoutput->setval("w_interference",(Double_t)w_interference);
741 +      return kTRUE;
742 +    }
743 +
744 +  //Otherwise continue.
745 +
746 +
747 +  
748 +  TLorentzVector lep1,lep2,lep3,lep4;
749 +  lep1.SetPtEtaPhiM(fs.gen->pt_1_a, fs.gen->eta_1_a, fs.gen->phi_1_a, 0.0); //Note: masses are set to 0, must now me important
750 +  lep2.SetPtEtaPhiM(fs.gen->pt_2_a, fs.gen->eta_2_a, fs.gen->phi_2_a, 0.0);
751 +  lep3.SetPtEtaPhiM(fs.gen->pt_1_b, fs.gen->eta_1_b, fs.gen->phi_1_b, 0.0);
752 +  lep4.SetPtEtaPhiM(fs.gen->pt_2_b, fs.gen->eta_2_b, fs.gen->phi_2_b, 0.0);
753 +
754 +
755 +  Float_t gen_costhetastar, gen_costheta1,gen_costheta2, gen_Phi, gen_Phi1;
756 +  Float_t gen_Z1Mass, gen_Z2Mass, gen_ZZMass;
757 +
758 +  TLorentzVector vec_zz, vec_a, vec_b;
759 +  vec_a.SetPtEtaPhiM(fs.gen->vpt_a, fs.gen->veta_a, fs.gen->vphi_a, fs.gen->vmass_a);
760 +  vec_b.SetPtEtaPhiM(fs.gen->vpt_b, fs.gen->veta_b, fs.gen->vphi_b, fs.gen->vmass_b);
761 +  vec_zz = vec_a+vec_b;
762 +
763 +  gen_Z1Mass = fs.gen->vmass_a;
764 +  gen_Z2Mass = fs.gen->vmass_b;
765 +  gen_ZZMass = vec_zz.M();
766 +
767 +  if (abs(fs.gen->id_1_a) == abs(fs.gen->id_1_b)) //Only make the weights for 4mu and 4e channels, not for 2e2mu
768 +    {
769 +      if (abs(fs.gen->id_2_a) != abs(fs.gen->id_1_a))
770 +        cout<<" *** Something is wrong - IDs don't match! "<<endl;
771 +
772 +
773 +      //First need to compute the angles in Mela's Framework using gen particles:
774 +      mela::computeAngles(lep1, fs.gen->id_1_a, lep2, fs.gen->id_2_a, lep3, fs.gen->id_1_b, lep4, fs.gen->id_2_b, gen_costhetastar, gen_costheta1, gen_costheta2, gen_Phi, gen_Phi1);
775 +      
776 +      //cout<<"gen mass zz = "<<vec_zz.M()<<"   gen mass a = "<<fs.gen->vmass_a<<"     gen mass b = "<<fs.gen->vmass_b<<endl;
777 +
778 +      //Now we can determine the weight:
779 +      mela1.computeWeight(gen_ZZMass, fs.gen->vmass_a, fs.gen->vmass_b, gen_costhetastar,gen_costheta1,gen_costheta2,gen_Phi,gen_Phi1, w_interference);
780 +    }
781 +
782 +  //Fill the weight into the ntuple:
783 +  flatoutput->setval("w_interference",(Double_t)w_interference);
784 +
785 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines