ScansICHEP2012/code/makeNtuple.C

//-----------------------------------------------
// Takes the text files and dumps the contents 
// in an ntuple.
// It's OK if a text file does not exist.
// However, the text files need to have the 
// same order of gluino-mass LSP-mass etc...
// Meaning that the lines must correspond
// If they do not, the program will barf
//
// It expects that there are files for 9 signal regions.
//
// The txt files are supposed to be "sorted".
// If they are not, you may want to sort them 
// ahead of time, eg:
//
// sort t1tttt_sr0.txt > SR0_sorted.txt
// sort t1tttt_sr1.txt > SR1_sorted.txt
// sort t1tttt_sr2.txt > SR2_sorted.txt
// sort t1tttt_sr3.txt > SR3_sorted.txt
// sort t1tttt_sr4.txt > SR4_sorted.txt
// sort t1tttt_sr5.txt > SR5_sorted.txt
// sort t1tttt_sr6.txt > SR6_sorted.txt
// sort t1tttt_sr7.txt > SR7_sorted.txt
// sort t1tttt_sr8.txt > SR8_sorted.txt
//
// Make sure that the file names, which are 
// hardwired (sorry) are correct.
//
// Also reads in the gluino xsection txt file.  
// Assumes that the gluino mass increases monotonically
// from one line to the next
//
//
//  Usage: 
//  root> .L Functions.cc
//  root> .L makeNtuple()
//  root> makeNtuple()
//------------------------------------------- 

//
//

void makeNtuple(int iflag=0){

  //------------------------------
  // Check the flags
  //------------------------------
  if (iflag < 1 || iflag > 6) {
    cout << "Usage:" << endl;
    cout << "makeNtuple(1)  for T1tttt" << endl;
    cout << "makeNtuple(2)  for glstop lsp50"  << endl;
    cout << "makeNtuple(3)  for glstop lsp250" << endl;
    cout << "makeNtuple(4)  for glsbottom chi150" << endl;
    cout << "makeNtuple(5)  for glsbottom chi300" << endl;
    cout << "makeNtuple(6)  for glsbottom sbpair" << endl;
    return;
  } 
  //---------------------------------------------------
  // This the model name, used for the root file name
  //----------------------------------------------------
  if (iflag == 1) modelName = "T1tttt";
  if (iflag == 2) modelName = "glstop50";
  if (iflag == 3) modelName = "glstop250";
  if (iflag == 4) modelName = "glsb150";
  if (iflag == 5) modelName = "glsb300";
  if (iflag == 6) modelName = "sbpair";

  //-----------------------------------------------
  // The name of the file with the output ntuple
  //-----------------------------------------------
  char* ntFile = Form("ntuple_%s.root",modelName);
  cout << "The output ntuple will be:   " << ntFile << endl;

  //-----------------------------------------------
  // Read in the gluino or sbottom pair production xsection  
  //-----------------------------------------------
  float xsec_mass[2000];
  float xsec[2000];
  float xsecup[2000];
  float xsecdwn[2000];
  int nxsec=0;
  char* xsecFile;
  if (iflag == 6) {
    xsecFile = "../xsec/xsec_sbottom.txt";
  } else {
    xsecFile = "../xsec/xsec_gluino.txt";
  }
  ifstream fromfile(Form("%s",xsecFile));
  if ( !fromfile.good() ) {
    cout << "gluino xsec file does not exist" << endl ;
    return;
  }
  cout << "Reading xsection from " << xsecFile << endl;
  //  cout << "ASSUMING THAT THE XSEC IS IN fb...TURN IT INTO pb" << endl;
  int xs=0;
  while (!fromfile.eof()) {
    TString d1, d2, d3, d4;
    fromfile >> d1 >> xsec_mass[xs] >> d2 >> xsec[xs] >> d3 
             >> xsecup[xs] >> d4 >> xsecdwn[xs];
    // xsec[xs]    = xsec[xs]/1000.;
    // xsecup[xs]  = xsecup[xs]/1000.;
    // xsecdwn[xs] = xsecdwn[xs]/1000.;
    xs++;
    nxsec++;
  }
  fromfile.close();
  nxsec = nxsec-1;  //last line is junk


  //for (int i=0; i<nxsec; i++) {
  // cout<<xsec_mass[i]<<" "<<xsec[i]<<" "<<xsecup[i]<<" "<< xsecdwn[i]<< endl;
  //}

  //--------------------------------------------------------
  // Check that the gluino mass in the xsection file increases monotonically
  //---------------------------------------------------------
  float previous = xsec_mass[0];
  for (int i=1; i<nxsec; i++) {
    if (xsec_mass[i] < previous) {
      cout << "Mass in xsection file does not increase monotonically..Abort"<<endl;
      return;
    }
    previous = xsec_mass[i];
  }

  //-----------------------------------------------
  // Define the variables in the txt files
  //-----------------------------------------------
  float glmass[9][1000];
  float lspmass[9][1000];
  float eff[9][1000];
  float err[9][1000];
  float obslim[9][1000];
  float explim[9][1000];
  float explimplus[9][1000];
  float explimminus[9][1000];
  bool  usedflag[9][1000];  // a flag to be used later
  int nentries[9];          // number of entries in the file
  
  //-----------------------------------------------
  // The file names, ordered by signal region...
  // Again, have to sort out the file names
  //-----------------------------------------------
  char* fileWhere;
  if (iflag == 1) fileWhere = "../T1tttt/L2_T1tttt";
  if (iflag == 2) fileWhere = "../glstop/lsp50/L2_lsp50";
  if (iflag == 3) fileWhere = "../glstop/lsp250/L2_lsp250";
  if (iflag == 4) fileWhere = "../glsbottom/chi150/L2_chi150";
  if (iflag == 5) fileWhere = "../glsbottom/chi300/L2_chi300";
  if (iflag == 6) fileWhere = "../sbottompair/L2_sbottompair";
  vector<TString> filenames;
  for (int kk=0; kk<9; kk++) {
    filenames.push_back(Form("%s_SR%d_CC.txt", fileWhere,kk));
    cout << "Will read from file " << filenames[kk] << endl;
  }
  //filenames.push_back("../sbottompair/L2_sbottompair_SR0_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR1_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR2_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR3_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR4_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR5_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR6_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR7_CC.txt");
  //filenames.push_back("../sbottompair/L2_sbottompair_SR8_CC.txt");


  //-----------------------------------------------
  // Loop over signal regions and load the arrays from the txt files
  //-----------------------------------------------
  for (int ireg=0; ireg<9; ireg++) {
    nentries[ireg] = 0;
    ifstream fromfile(filenames.at(ireg));

    // if the file does not exist, go to the next one
    if ( !fromfile.good() ) {
      cout << "file " << filenames.at(ireg) << " does not exist" << endl ;
      continue;
    }
    // the file exist, read from it
    cout << "Reading from " << filenames.at(ireg) << endl;
    int j=0;
    while (!fromfile.eof()) {
      fromfile >> glmass[ireg][j]
               >> lspmass[ireg][j]
               >> eff[ireg][j]
               >> err[ireg][j];
      j++;
      nentries[ireg]++;
    }
    fromfile.close();
    // the last entry is junk  
    nentries[ireg] = nentries[ireg]-1;
  }
  //-------------------------------------------
  // Now check that the line order is OK
  // And that the number of lines is consistent
  //--------------------------------------------
  bool bad=false;
  int numbLines = 0;
  for (int ii=0; ii<9; ii++) {
    for (int jj=ii+1; jj<9; jj++) {
      if (nentries[ii]*nentries[jj] !=0 ) {
        numbLines = nentries[ii];
        if (nentries[ii] != nentries[jj]) {
          cout << "Mismatched lines: file " << filenames.at(ii);
          cout << " has " << nentries[ii] << " lines" << endl;
          cout << "                  file " << filenames.at(jj);
          cout << " has " << nentries[jj] << " lines" << endl;
          bad=true;
        }
        for (int line=0; line<nentries[ii]; line++) {
          if (glmass[ii][line] != glmass[jj][line]) {
            cout << "Gluino mass on line " << line  
                  << " of files " << filenames.at(ii) << " and " 
                 << filenames.at(jj) << " do not match" << endl;
            bad=true;
          }
          if (lspmass[ii][line] != lspmass[jj][line]) {
            cout << "LSP mass on line " << line 
                 << " of files " << filenames.at(ii) << " and " 
                 << filenames.at(jj) << " do not match" << endl;
            bad=true;
          }
        }
      }
    }
  }
  if (bad) return;
  cout << " The file formats match...Now fill ntuple" <<endl;

  //-------------------------------------
  //Calculate the total uncertainty
  //--------------------------------------
  float lumiErr   = 0.045;
  float trgErr    = 0.030;
  float isoErr    = 0.100;  // 5% per lepton
  float fixedErr2 = lumiErr*lumiErr + trgErr*trgErr + isoErr*isoErr;
  for (int il=0; il<numbLines; il++) {
    for (int sr=0; sr<9; sr++) {
      err[sr][il] = sqrt(err[sr][il]*err[sr][il] + fixedErr2);
    }
  }
  //-----------------------------------------------------
  // And now calculate the limits on the fly
  //----------------------------------------------------
    for (int sr=0; sr<9; sr++) {
      cout << "Calculating limits for SR" << sr << endl;
      for (int il=0; il<numbLines; il++) {
      float  thisErr      = err[sr][il];
      double blahObs      = observedLimitOnNumberOfEvents(sr+1, 100.*thisErr, true  );
      double blahExp      = expectedLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
      double blahExpplus  = expectedPlusOneSigmaLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
      double blahExpminus = expectedMinusOneSigmaLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
      obslim[sr][il]      = (float) blahObs;
      explim[sr][il]      = (float) blahExp;      
      explimplus[sr][il]  = (float) blahExpplus;      
      explimminus[sr][il] = (float) blahExpminus;      
    }
  }
  //-----------------------------------------------------
  // Now the content of the files is loaded in arrays
  // We are going to stick everything in an ntuple
  //-----------------------------------------------------
  TFile* babyFile_ = TFile::Open(Form("%s", ntFile), "RECREATE");
  babyFile_->cd();
  babyTree_ = new TTree("tree", "A Baby Ntuple");
  babyTree_->SetDirectory(0);

  // define the variables
  float glmass_;       // gluino mass
  float lspmass_;      // lsp mass
  float xsec_;         // xsec
  float xsecup_;       // xsec, one sigma high
  float xsecdwn_;      // xsec, one sigma down
  int   bestsr_ = 0;   // best signal region

  // Here come the efficiencies for the various SR.  
  // The last one is the "best", based on the best expected limit 
  float effsr0_ = -1.;
  float effsr1_ = -1.;
  float effsr2_ = -1.;
  float effsr3_ = -1.;
  float effsr4_ = -1.;
  float effsr5_ = -1.;
  float effsr6_ = -1.;
  float effsr7_ = -1.;
  float effsr8_ = -1.;
  float effsrb_ = -1.;    // b = best SR

  // same as above but for the efficiency error
  float errsr0_ = -1.;
  float errsr1_ = -1.;
  float errsr2_ = -1.;
  float errsr3_ = -1.;
  float errsr4_ = -1.;
  float errsr5_ = -1.;
  float errsr6_ = -1.;
  float errsr7_ = -1.;
  float errsr8_ = -1.;
  float errsrb_ = -1.;   // b = best SR

  // same as above but for observed limit
  float obslimsr0_ = -1.;
  float obslimsr1_ = -1.;
  float obslimsr2_ = -1.;
  float obslimsr3_ = -1.;
  float obslimsr4_ = -1.;
  float obslimsr5_ = -1.;
  float obslimsr6_ = -1.;
  float obslimsr7_ = -1.;
  float obslimsr8_ = -1.;
  float obslimsrb_ = -1.;   // b = best SR

  // same as above but for expected limit
  float explimsr0_ = -1.;
  float explimsr1_ = -1.;
  float explimsr2_ = -1.;
  float explimsr3_ = -1.;
  float explimsr4_ = -1.;
  float explimsr5_ = -1.;
  float explimsr6_ = -1.;
  float explimsr7_ = -1.;
  float explimsr8_ = -1.;
  float explimsrb_ = -1.;   // b = best SR

  float explimplussr0_ = -1.;
  float explimplussr1_ = -1.;
  float explimplussr2_ = -1.;  
  float explimplussr3_ = -1.;
  float explimplussr4_ = -1.;
  float explimplussr5_ = -1.;
  float explimplussr6_ = -1.;
  float explimplussr7_ = -1.;
  float explimplussr8_ = -1.;
  float explimplussrb_ = -1.;   // b = best SR

  float explimminussr0_ = -1.;
  float explimminussr1_ = -1.;
  float explimminussr2_ = -1.;  
  float explimminussr3_ = -1.;
  float explimminussr4_ = -1.;
  float explimminussr5_ = -1.;
  float explimminussr6_ = -1.;
  float explimminussr7_ = -1.;
  float explimminussr8_ = -1.;
  float explimminussrb_ = -1.;   // b = best SR


  // put the branches in the tree
  babyTree_->Branch("glmass",  &glmass_);
  babyTree_->Branch("lspmass", &lspmass_);
  babyTree_->Branch("xsec",    &xsec_);
  babyTree_->Branch("xsecup",  &xsecup_);
  babyTree_->Branch("xsecdwn", &xsecdwn_);
  babyTree_->Branch("bestsr",  &bestsr_);


  babyTree_->Branch("effsr0", &effsr0_);
  babyTree_->Branch("effsr1", &effsr1_);
  babyTree_->Branch("effsr2", &effsr2_);
  babyTree_->Branch("effsr3", &effsr3_);
  babyTree_->Branch("effsr4", &effsr4_);
  babyTree_->Branch("effsr5", &effsr5_);
  babyTree_->Branch("effsr6", &effsr6_);
  babyTree_->Branch("effsr7", &effsr7_);
  babyTree_->Branch("effsr8", &effsr8_);
  babyTree_->Branch("effsrb", &effsrb_);   // b = best SR

  babyTree_->Branch("errsr0", &errsr0_);
  babyTree_->Branch("errsr1", &errsr1_);
  babyTree_->Branch("errsr2", &errsr2_);
  babyTree_->Branch("errsr3", &errsr3_);
  babyTree_->Branch("errsr4", &errsr4_);
  babyTree_->Branch("errsr5", &errsr5_);
  babyTree_->Branch("errsr6", &errsr6_);
  babyTree_->Branch("errsr7", &errsr7_);
  babyTree_->Branch("errsr8", &errsr8_);
  babyTree_->Branch("errsrb", &errsrb_);   // b = best SR

  babyTree_->Branch("obslimsr0", &obslimsr0_);
  babyTree_->Branch("obslimsr1", &obslimsr1_);
  babyTree_->Branch("obslimsr2", &obslimsr2_);
  babyTree_->Branch("obslimsr3", &obslimsr3_);
  babyTree_->Branch("obslimsr4", &obslimsr4_);
  babyTree_->Branch("obslimsr5", &obslimsr5_);
  babyTree_->Branch("obslimsr6", &obslimsr6_);
  babyTree_->Branch("obslimsr7", &obslimsr7_);
  babyTree_->Branch("obslimsr8", &obslimsr8_);
  babyTree_->Branch("obslimsrb", &obslimsrb_);   // b = best SR

  babyTree_->Branch("explimsr0", &explimsr0_);
  babyTree_->Branch("explimsr1", &explimsr1_);
  babyTree_->Branch("explimsr2", &explimsr2_);
  babyTree_->Branch("explimsr3", &explimsr3_);
  babyTree_->Branch("explimsr4", &explimsr4_);
  babyTree_->Branch("explimsr5", &explimsr5_);
  babyTree_->Branch("explimsr6", &explimsr6_);
  babyTree_->Branch("explimsr7", &explimsr7_);
  babyTree_->Branch("explimsr8", &explimsr8_);
  babyTree_->Branch("explimsrb", &explimsrb_);  // b = best SR

  babyTree_->Branch("explimplussr0", &explimplussr0_);
  babyTree_->Branch("explimplussr1", &explimplussr1_);
  babyTree_->Branch("explimplussr2", &explimplussr2_);
  babyTree_->Branch("explimplussr3", &explimplussr3_);
  babyTree_->Branch("explimplussr4", &explimplussr4_);
  babyTree_->Branch("explimplussr5", &explimplussr5_);
  babyTree_->Branch("explimplussr6", &explimplussr6_);
  babyTree_->Branch("explimplussr7", &explimplussr7_);
  babyTree_->Branch("explimplussr8", &explimplussr8_);
  babyTree_->Branch("explimplussrb", &explimplussrb_);  // b = best SR

  babyTree_->Branch("explimminussr0", &explimminussr0_);
  babyTree_->Branch("explimminussr1", &explimminussr1_);
  babyTree_->Branch("explimminussr2", &explimminussr2_);
  babyTree_->Branch("explimminussr3", &explimminussr3_);
  babyTree_->Branch("explimminussr4", &explimminussr4_);
  babyTree_->Branch("explimminussr5", &explimminussr5_);
  babyTree_->Branch("explimminussr6", &explimminussr6_);
  babyTree_->Branch("explimminussr7", &explimminussr7_);
  babyTree_->Branch("explimminussr8", &explimminussr8_);
  babyTree_->Branch("explimminussrb", &explimminussrb_);  // b = best SR


  // fill the variables.  
  // Note: the best region is the one with the smallest
  // ratio of expected limit and efficiency
  cout << "Filling an ntuple with " << numbLines << " entries" << endl;
  cout << "SR0 is not allowed to be the best region" << endl;
  for (int il=0; il<numbLines; il++) {
    float best = 999999.;
    int ibest = -1;

    if (nentries[0] != 0) {
      glmass_  = glmass[0][il];
      lspmass_ = lspmass[0][il];
      effsr0_ = eff[0][il];
      errsr0_ = eff[0][il];
      obslimsr0_ = obslim[0][il];
      explimsr0_ = explim[0][il];
      explimplussr0_  = explimplus[0][il];
      explimminussr0_ = explimminus[0][il];
      float temp = explimsr0_/effsr0_;
      //  if (temp < best) {
      //        best = temp;
      //        ibest = 0;
      // }
    }
    if (nentries[1] != 0) {
      glmass_  = glmass[1][il];
      lspmass_ = lspmass[1][il];
      effsr1_ = eff[1][il];
      errsr1_ = eff[1][il];
      obslimsr1_ = obslim[1][il];
      explimsr1_ = explim[1][il];
      explimplussr1_  = explimplus[1][il];
      explimminussr1_ = explimminus[1][il];
      float temp = explimsr1_/effsr1_;
      if (temp < best) {
        best = temp;
        ibest = 1;
      }
    }
    if (nentries[2] != 0) {
      glmass_  = glmass[2][il];
      lspmass_ = lspmass[2][il];
      effsr2_ = eff[2][il];
      errsr2_ = eff[2][il];
      obslimsr2_ = obslim[2][il];
      explimsr2_ = explim[2][il];
      explimplussr2_  = explimplus[2][il];
      explimminussr2_ = explimminus[2][il];
      float temp = explimsr2_/effsr2_;
      if (temp < best) {
        best = temp;
        ibest = 2;
      }
    }
    if (nentries[3] != 0) {
      glmass_  = glmass[3][il];
      lspmass_ = lspmass[3][il];
      effsr3_ = eff[3][il];
      errsr3_ = eff[3][il];
      obslimsr3_ = obslim[3][il];
      explimsr3_ = explim[3][il];
      explimplussr3_  = explimplus[3][il];
      explimminussr3_ = explimminus[3][il];
      float temp = explimsr3_/effsr3_;
      if (temp < best) {
        best = temp;
        ibest = 3;
      }
    }
    if (nentries[4] != 0) {
      glmass_  = glmass[4][il];
      lspmass_ = lspmass[4][il];
      effsr4_ = eff[4][il];
      errsr4_ = eff[4][il];
      obslimsr4_ = obslim[4][il];
      explimsr4_ = explim[4][il];
      explimplussr4_  = explimplus[4][il];
      explimminussr4_ = explimminus[4][il];
      float temp = explimsr4_/effsr4_;
      if (temp < best) {
        best = temp;
        ibest = 4;
      }
    }
    if (nentries[5] != 0) {
      glmass_  = glmass[5][il];
      lspmass_ = lspmass[5][il];
      effsr5_ = eff[5][il];
      errsr5_ = eff[5][il];
      obslimsr5_ = obslim[5][il];
      explimsr5_ = explim[5][il];
      explimplussr5_  = explimplus[5][il];
      explimminussr5_ = explimminus[5][il];
      float temp = explimsr5_/effsr5_;
      if (temp < best) {
        best = temp;
        ibest = 5;
      }
    }
    if (nentries[6] != 0) {
      glmass_  = glmass[6][il];
      lspmass_ = lspmass[6][il];
      effsr6_ = eff[6][il];
      errsr6_ = eff[6][il];
      obslimsr6_ = obslim[6][il];
      explimsr6_ = explim[6][il];
      explimplussr6_  = explimplus[6][il];
      explimminussr6_ = explimminus[6][il];
      float temp = explimsr6_/effsr6_;
      if (temp < best) {
        best = temp;
        ibest = 6;
      }
    }
    if (nentries[7] != 0) {
      glmass_  = glmass[7][il];
      lspmass_ = lspmass[7][il];
      effsr7_ = eff[7][il];
      errsr7_ = eff[7][il];
      obslimsr7_ = obslim[7][il];
      explimsr7_ = explim[7][il];
      explimplussr7_  = explimplus[7][il];
      explimminussr7_ = explimminus[7][il];
      float temp = explimsr7_/effsr7_;
      if (temp < best) {
        best = temp;
        ibest = 7;
      }
    }
    if (nentries[8] != 0) {
      glmass_  = glmass[8][il];
      lspmass_ = lspmass[8][il];
      effsr8_ = eff[8][il];
      errsr8_ = eff[8][il];
      obslimsr8_ = obslim[8][il];
      explimsr8_ = explim[8][il];
      explimplussr8_  = explimplus[8][il];
      explimminussr8_ = explimminus[8][il];
      float temp = explimsr8_/effsr8_;
      if (temp < best) {
        best = temp;
        ibest = 8;
      }
    }

    bestsr_ = ibest;
    effsrb_ = eff[ibest][il];
    errsrb_ = eff[ibest][il];
    obslimsrb_ = obslim[ibest][il];
    explimsrb_ = explim[ibest][il];
    explimplussrb_ = explimplus[ibest][il];
    explimminussrb_ = explimminus[ibest][il];

    // Fill the cross-section
    bool done = false;
    for (int xs=0; xs<nxsec-1; xs++) {
      if ( (glmass_ >= xsec_mass[xs] && glmass_ < xsec_mass[xs+1]) ||
           (glmass_ <  xsec_mass[xs] && xs==0) ) {

        float dd1 = glmass_ - xsec_mass[xs];
        float dd2 = xsec_mass[xs+1] - xsec_mass[xs];
        xsec_     = xsec[xs]    + (xsec[xs+1]    - xsec[xs])   *(dd1/dd2);
        xsecup_   = xsecup[xs]  + (xsecup[xs+1]  - xsecup[xs]) *(dd1/dd2);
        xsecdwn_  = xsecdwn[xs] + (xsecdwn[xs+1] - xsecdwn[xs])*(dd1/dd2);
        done = true;
        break;
      }
    }
    if (!done) {
      int xs = nxsec - 1; 
      float dd1 = glmass_ - xsec_mass[xs];
      float dd2 = xsec_mass[xs] - xsec_mass[xs-1];
      xsec_     = xsec[xs]    + (xsec[xs]    - xsec[xs-1])   *(dd1/dd2);
      xsecup_   = xsecup[xs]  + (xsecup[xs]  - xsecup[xs-1]) *(dd1/dd2);
      xsecdwn_  = xsecdwn[xs] + (xsecdwn[xs] - xsecdwn[xs-1])*(dd1/dd2);
    }


    // Now fill the ntuple
    babyTree_->Fill(); 

  }  
  // Now close the file

    babyFile_->cd();
    babyTree_->Write();
    babyFile_->Close();


}

Revision:	1.4
Committed:	Mon Jul 2 04:42:05 2012 UTC (12 years, 10 months ago) by claudioc
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+76 -5 lines
Error occurred while calculating annotation data.
Log Message:	lipstick
#	Content
1	//-----------------------------------------------
2	// Takes the text files and dumps the contents
3	// in an ntuple.
4	// It's OK if a text file does not exist.
5	// However, the text files need to have the
6	// same order of gluino-mass LSP-mass etc...
7	// Meaning that the lines must correspond
8	// If they do not, the program will barf
9	//
10	// It expects that there are files for 9 signal regions.
11	//
12	// The txt files are supposed to be "sorted".
13	// If they are not, you may want to sort them
14	// ahead of time, eg:
15	//
16	// sort t1tttt_sr0.txt > SR0_sorted.txt
17	// sort t1tttt_sr1.txt > SR1_sorted.txt
18	// sort t1tttt_sr2.txt > SR2_sorted.txt
19	// sort t1tttt_sr3.txt > SR3_sorted.txt
20	// sort t1tttt_sr4.txt > SR4_sorted.txt
21	// sort t1tttt_sr5.txt > SR5_sorted.txt
22	// sort t1tttt_sr6.txt > SR6_sorted.txt
23	// sort t1tttt_sr7.txt > SR7_sorted.txt
24	// sort t1tttt_sr8.txt > SR8_sorted.txt
25	//
26	// Make sure that the file names, which are
27	// hardwired (sorry) are correct.
28	//
29	// Also reads in the gluino xsection txt file.
30	// Assumes that the gluino mass increases monotonically
31	// from one line to the next
32	//
33	//
34	// Usage:
35	// root> .L Functions.cc
36	// root> .L makeNtuple()
37	// root> makeNtuple()
38	//-------------------------------------------
39
40	//
41	//
42
43	void makeNtuple(int iflag=0){
44
45	//------------------------------
46	// Check the flags
47	//------------------------------
48	if (iflag < 1 \|\| iflag > 6) {
49	cout << "Usage:" << endl;
50	cout << "makeNtuple(1) for T1tttt" << endl;
51	cout << "makeNtuple(2) for glstop lsp50" << endl;
52	cout << "makeNtuple(3) for glstop lsp250" << endl;
53	cout << "makeNtuple(4) for glsbottom chi150" << endl;
54	cout << "makeNtuple(5) for glsbottom chi300" << endl;
55	cout << "makeNtuple(6) for glsbottom sbpair" << endl;
56	return;
57	}
58	//---------------------------------------------------
59	// This the model name, used for the root file name
60	//----------------------------------------------------
61	if (iflag == 1) modelName = "T1tttt";
62	if (iflag == 2) modelName = "glstop50";
63	if (iflag == 3) modelName = "glstop250";
64	if (iflag == 4) modelName = "glsb150";
65	if (iflag == 5) modelName = "glsb300";
66	if (iflag == 6) modelName = "sbpair";
67
68	//-----------------------------------------------
69	// The name of the file with the output ntuple
70	//-----------------------------------------------
71	char* ntFile = Form("ntuple_%s.root",modelName);
72	cout << "The output ntuple will be: " << ntFile << endl;
73
74	//-----------------------------------------------
75	// Read in the gluino or sbottom pair production xsection
76	//-----------------------------------------------
77	float xsec_mass[2000];
78	float xsec[2000];
79	float xsecup[2000];
80	float xsecdwn[2000];
81	int nxsec=0;
82	char* xsecFile;
83	if (iflag == 6) {
84	xsecFile = "../xsec/xsec_sbottom.txt";
85	} else {
86	xsecFile = "../xsec/xsec_gluino.txt";
87	}
88	ifstream fromfile(Form("%s",xsecFile));
89	if ( !fromfile.good() ) {
90	cout << "gluino xsec file does not exist" << endl ;
91	return;
92	}
93	cout << "Reading xsection from " << xsecFile << endl;
94	// cout << "ASSUMING THAT THE XSEC IS IN fb...TURN IT INTO pb" << endl;
95	int xs=0;
96	while (!fromfile.eof()) {
97	TString d1, d2, d3, d4;
98	fromfile >> d1 >> xsec_mass[xs] >> d2 >> xsec[xs] >> d3
99	>> xsecup[xs] >> d4 >> xsecdwn[xs];
100	// xsec[xs] = xsec[xs]/1000.;
101	// xsecup[xs] = xsecup[xs]/1000.;
102	// xsecdwn[xs] = xsecdwn[xs]/1000.;
103	xs++;
104	nxsec++;
105	}
106	fromfile.close();
107	nxsec = nxsec-1; //last line is junk
108
109
110	//for (int i=0; i<nxsec; i++) {
111	// cout<<xsec_mass[i]<<" "<<xsec[i]<<" "<<xsecup[i]<<" "<< xsecdwn[i]<< endl;
112	//}
113
114	//--------------------------------------------------------
115	// Check that the gluino mass in the xsection file increases monotonically
116	//---------------------------------------------------------
117	float previous = xsec_mass[0];
118	for (int i=1; i<nxsec; i++) {
119	if (xsec_mass[i] < previous) {
120	cout << "Mass in xsection file does not increase monotonically..Abort"<<endl;
121	return;
122	}
123	previous = xsec_mass[i];
124	}
125
126	//-----------------------------------------------
127	// Define the variables in the txt files
128	//-----------------------------------------------
129	float glmass[9][1000];
130	float lspmass[9][1000];
131	float eff[9][1000];
132	float err[9][1000];
133	float obslim[9][1000];
134	float explim[9][1000];
135	float explimplus[9][1000];
136	float explimminus[9][1000];
137	bool usedflag[9][1000]; // a flag to be used later
138	int nentries[9]; // number of entries in the file
139
140	//-----------------------------------------------
141	// The file names, ordered by signal region...
142	// Again, have to sort out the file names
143	//-----------------------------------------------
144	char* fileWhere;
145	if (iflag == 1) fileWhere = "../T1tttt/L2_T1tttt";
146	if (iflag == 2) fileWhere = "../glstop/lsp50/L2_lsp50";
147	if (iflag == 3) fileWhere = "../glstop/lsp250/L2_lsp250";
148	if (iflag == 4) fileWhere = "../glsbottom/chi150/L2_chi150";
149	if (iflag == 5) fileWhere = "../glsbottom/chi300/L2_chi300";
150	if (iflag == 6) fileWhere = "../sbottompair/L2_sbottompair";
151	vector<TString> filenames;
152	for (int kk=0; kk<9; kk++) {
153	filenames.push_back(Form("%s_SR%d_CC.txt", fileWhere,kk));
154	cout << "Will read from file " << filenames[kk] << endl;
155	}
156	//filenames.push_back("../sbottompair/L2_sbottompair_SR0_CC.txt");
157	//filenames.push_back("../sbottompair/L2_sbottompair_SR1_CC.txt");
158	//filenames.push_back("../sbottompair/L2_sbottompair_SR2_CC.txt");
159	//filenames.push_back("../sbottompair/L2_sbottompair_SR3_CC.txt");
160	//filenames.push_back("../sbottompair/L2_sbottompair_SR4_CC.txt");
161	//filenames.push_back("../sbottompair/L2_sbottompair_SR5_CC.txt");
162	//filenames.push_back("../sbottompair/L2_sbottompair_SR6_CC.txt");
163	//filenames.push_back("../sbottompair/L2_sbottompair_SR7_CC.txt");
164	//filenames.push_back("../sbottompair/L2_sbottompair_SR8_CC.txt");
165
166
167	//-----------------------------------------------
168	// Loop over signal regions and load the arrays from the txt files
169	//-----------------------------------------------
170	for (int ireg=0; ireg<9; ireg++) {
171	nentries[ireg] = 0;
172	ifstream fromfile(filenames.at(ireg));
173
174	// if the file does not exist, go to the next one
175	if ( !fromfile.good() ) {
176	cout << "file " << filenames.at(ireg) << " does not exist" << endl ;
177	continue;
178	}
179	// the file exist, read from it
180	cout << "Reading from " << filenames.at(ireg) << endl;
181	int j=0;
182	while (!fromfile.eof()) {
183	fromfile >> glmass[ireg][j]
184	>> lspmass[ireg][j]
185	>> eff[ireg][j]
186	>> err[ireg][j];
187	j++;
188	nentries[ireg]++;
189	}
190	fromfile.close();
191	// the last entry is junk
192	nentries[ireg] = nentries[ireg]-1;
193	}
194	//-------------------------------------------
195	// Now check that the line order is OK
196	// And that the number of lines is consistent
197	//--------------------------------------------
198	bool bad=false;
199	int numbLines = 0;
200	for (int ii=0; ii<9; ii++) {
201	for (int jj=ii+1; jj<9; jj++) {
202	if (nentries[ii]*nentries[jj] !=0 ) {
203	numbLines = nentries[ii];
204	if (nentries[ii] != nentries[jj]) {
205	cout << "Mismatched lines: file " << filenames.at(ii);
206	cout << " has " << nentries[ii] << " lines" << endl;
207	cout << " file " << filenames.at(jj);
208	cout << " has " << nentries[jj] << " lines" << endl;
209	bad=true;
210	}
211	for (int line=0; line<nentries[ii]; line++) {
212	if (glmass[ii][line] != glmass[jj][line]) {
213	cout << "Gluino mass on line " << line
214	<< " of files " << filenames.at(ii) << " and "
215	<< filenames.at(jj) << " do not match" << endl;
216	bad=true;
217	}
218	if (lspmass[ii][line] != lspmass[jj][line]) {
219	cout << "LSP mass on line " << line
220	<< " of files " << filenames.at(ii) << " and "
221	<< filenames.at(jj) << " do not match" << endl;
222	bad=true;
223	}
224	}
225	}
226	}
227	}
228	if (bad) return;
229	cout << " The file formats match...Now fill ntuple" <<endl;
230
231	//-------------------------------------
232	//Calculate the total uncertainty
233	//--------------------------------------
234	float lumiErr = 0.045;
235	float trgErr = 0.030;
236	float isoErr = 0.100; // 5% per lepton
237	float fixedErr2 = lumiErrlumiErr + trgErrtrgErr + isoErr*isoErr;
238	for (int il=0; il<numbLines; il++) {
239	for (int sr=0; sr<9; sr++) {
240	err[sr][il] = sqrt(err[sr][il]*err[sr][il] + fixedErr2);
241	}
242	}
243	//-----------------------------------------------------
244	// And now calculate the limits on the fly
245	//----------------------------------------------------
246	for (int sr=0; sr<9; sr++) {
247	cout << "Calculating limits for SR" << sr << endl;
248	for (int il=0; il<numbLines; il++) {
249	float thisErr = err[sr][il];
250	double blahObs = observedLimitOnNumberOfEvents(sr+1, 100.*thisErr, true );
251	double blahExp = expectedLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
252	double blahExpplus = expectedPlusOneSigmaLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
253	double blahExpminus = expectedMinusOneSigmaLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
254	obslim[sr][il] = (float) blahObs;
255	explim[sr][il] = (float) blahExp;
256	explimplus[sr][il] = (float) blahExpplus;
257	explimminus[sr][il] = (float) blahExpminus;
258	}
259	}
260	//-----------------------------------------------------
261	// Now the content of the files is loaded in arrays
262	// We are going to stick everything in an ntuple
263	//-----------------------------------------------------
264	TFile* babyFile_ = TFile::Open(Form("%s", ntFile), "RECREATE");
265	babyFile_->cd();
266	babyTree_ = new TTree("tree", "A Baby Ntuple");
267	babyTree_->SetDirectory(0);
268
269	// define the variables
270	float glmass_; // gluino mass
271	float lspmass_; // lsp mass
272	float xsec_; // xsec
273	float xsecup_; // xsec, one sigma high
274	float xsecdwn_; // xsec, one sigma down
275	int bestsr_ = 0; // best signal region
276
277	// Here come the efficiencies for the various SR.
278	// The last one is the "best", based on the best expected limit
279	float effsr0_ = -1.;
280	float effsr1_ = -1.;
281	float effsr2_ = -1.;
282	float effsr3_ = -1.;
283	float effsr4_ = -1.;
284	float effsr5_ = -1.;
285	float effsr6_ = -1.;
286	float effsr7_ = -1.;
287	float effsr8_ = -1.;
288	float effsrb_ = -1.; // b = best SR
289
290	// same as above but for the efficiency error
291	float errsr0_ = -1.;
292	float errsr1_ = -1.;
293	float errsr2_ = -1.;
294	float errsr3_ = -1.;
295	float errsr4_ = -1.;
296	float errsr5_ = -1.;
297	float errsr6_ = -1.;
298	float errsr7_ = -1.;
299	float errsr8_ = -1.;
300	float errsrb_ = -1.; // b = best SR
301
302	// same as above but for observed limit
303	float obslimsr0_ = -1.;
304	float obslimsr1_ = -1.;
305	float obslimsr2_ = -1.;
306	float obslimsr3_ = -1.;
307	float obslimsr4_ = -1.;
308	float obslimsr5_ = -1.;
309	float obslimsr6_ = -1.;
310	float obslimsr7_ = -1.;
311	float obslimsr8_ = -1.;
312	float obslimsrb_ = -1.; // b = best SR
313
314	// same as above but for expected limit
315	float explimsr0_ = -1.;
316	float explimsr1_ = -1.;
317	float explimsr2_ = -1.;
318	float explimsr3_ = -1.;
319	float explimsr4_ = -1.;
320	float explimsr5_ = -1.;
321	float explimsr6_ = -1.;
322	float explimsr7_ = -1.;
323	float explimsr8_ = -1.;
324	float explimsrb_ = -1.; // b = best SR
325
326	float explimplussr0_ = -1.;
327	float explimplussr1_ = -1.;
328	float explimplussr2_ = -1.;
329	float explimplussr3_ = -1.;
330	float explimplussr4_ = -1.;
331	float explimplussr5_ = -1.;
332	float explimplussr6_ = -1.;
333	float explimplussr7_ = -1.;
334	float explimplussr8_ = -1.;
335	float explimplussrb_ = -1.; // b = best SR
336
337	float explimminussr0_ = -1.;
338	float explimminussr1_ = -1.;
339	float explimminussr2_ = -1.;
340	float explimminussr3_ = -1.;
341	float explimminussr4_ = -1.;
342	float explimminussr5_ = -1.;
343	float explimminussr6_ = -1.;
344	float explimminussr7_ = -1.;
345	float explimminussr8_ = -1.;
346	float explimminussrb_ = -1.; // b = best SR
347
348
349	// put the branches in the tree
350	babyTree_->Branch("glmass", &glmass_);
351	babyTree_->Branch("lspmass", &lspmass_);
352	babyTree_->Branch("xsec", &xsec_);
353	babyTree_->Branch("xsecup", &xsecup_);
354	babyTree_->Branch("xsecdwn", &xsecdwn_);
355	babyTree_->Branch("bestsr", &bestsr_);
356
357
358	babyTree_->Branch("effsr0", &effsr0_);
359	babyTree_->Branch("effsr1", &effsr1_);
360	babyTree_->Branch("effsr2", &effsr2_);
361	babyTree_->Branch("effsr3", &effsr3_);
362	babyTree_->Branch("effsr4", &effsr4_);
363	babyTree_->Branch("effsr5", &effsr5_);
364	babyTree_->Branch("effsr6", &effsr6_);
365	babyTree_->Branch("effsr7", &effsr7_);
366	babyTree_->Branch("effsr8", &effsr8_);
367	babyTree_->Branch("effsrb", &effsrb_); // b = best SR
368
369	babyTree_->Branch("errsr0", &errsr0_);
370	babyTree_->Branch("errsr1", &errsr1_);
371	babyTree_->Branch("errsr2", &errsr2_);
372	babyTree_->Branch("errsr3", &errsr3_);
373	babyTree_->Branch("errsr4", &errsr4_);
374	babyTree_->Branch("errsr5", &errsr5_);
375	babyTree_->Branch("errsr6", &errsr6_);
376	babyTree_->Branch("errsr7", &errsr7_);
377	babyTree_->Branch("errsr8", &errsr8_);
378	babyTree_->Branch("errsrb", &errsrb_); // b = best SR
379
380	babyTree_->Branch("obslimsr0", &obslimsr0_);
381	babyTree_->Branch("obslimsr1", &obslimsr1_);
382	babyTree_->Branch("obslimsr2", &obslimsr2_);
383	babyTree_->Branch("obslimsr3", &obslimsr3_);
384	babyTree_->Branch("obslimsr4", &obslimsr4_);
385	babyTree_->Branch("obslimsr5", &obslimsr5_);
386	babyTree_->Branch("obslimsr6", &obslimsr6_);
387	babyTree_->Branch("obslimsr7", &obslimsr7_);
388	babyTree_->Branch("obslimsr8", &obslimsr8_);
389	babyTree_->Branch("obslimsrb", &obslimsrb_); // b = best SR
390
391	babyTree_->Branch("explimsr0", &explimsr0_);
392	babyTree_->Branch("explimsr1", &explimsr1_);
393	babyTree_->Branch("explimsr2", &explimsr2_);
394	babyTree_->Branch("explimsr3", &explimsr3_);
395	babyTree_->Branch("explimsr4", &explimsr4_);
396	babyTree_->Branch("explimsr5", &explimsr5_);
397	babyTree_->Branch("explimsr6", &explimsr6_);
398	babyTree_->Branch("explimsr7", &explimsr7_);
399	babyTree_->Branch("explimsr8", &explimsr8_);
400	babyTree_->Branch("explimsrb", &explimsrb_); // b = best SR
401
402	babyTree_->Branch("explimplussr0", &explimplussr0_);
403	babyTree_->Branch("explimplussr1", &explimplussr1_);
404	babyTree_->Branch("explimplussr2", &explimplussr2_);
405	babyTree_->Branch("explimplussr3", &explimplussr3_);
406	babyTree_->Branch("explimplussr4", &explimplussr4_);
407	babyTree_->Branch("explimplussr5", &explimplussr5_);
408	babyTree_->Branch("explimplussr6", &explimplussr6_);
409	babyTree_->Branch("explimplussr7", &explimplussr7_);
410	babyTree_->Branch("explimplussr8", &explimplussr8_);
411	babyTree_->Branch("explimplussrb", &explimplussrb_); // b = best SR
412
413	babyTree_->Branch("explimminussr0", &explimminussr0_);
414	babyTree_->Branch("explimminussr1", &explimminussr1_);
415	babyTree_->Branch("explimminussr2", &explimminussr2_);
416	babyTree_->Branch("explimminussr3", &explimminussr3_);
417	babyTree_->Branch("explimminussr4", &explimminussr4_);
418	babyTree_->Branch("explimminussr5", &explimminussr5_);
419	babyTree_->Branch("explimminussr6", &explimminussr6_);
420	babyTree_->Branch("explimminussr7", &explimminussr7_);
421	babyTree_->Branch("explimminussr8", &explimminussr8_);
422	babyTree_->Branch("explimminussrb", &explimminussrb_); // b = best SR
423
424
425	// fill the variables.
426	// Note: the best region is the one with the smallest
427	// ratio of expected limit and efficiency
428	cout << "Filling an ntuple with " << numbLines << " entries" << endl;
429	cout << "SR0 is not allowed to be the best region" << endl;
430	for (int il=0; il<numbLines; il++) {
431	float best = 999999.;
432	int ibest = -1;
433
434	if (nentries[0] != 0) {
435	glmass_ = glmass[0][il];
436	lspmass_ = lspmass[0][il];
437	effsr0_ = eff[0][il];
438	errsr0_ = eff[0][il];
439	obslimsr0_ = obslim[0][il];
440	explimsr0_ = explim[0][il];
441	explimplussr0_ = explimplus[0][il];
442	explimminussr0_ = explimminus[0][il];
443	float temp = explimsr0_/effsr0_;
444	// if (temp < best) {
445	// best = temp;
446	// ibest = 0;
447	// }
448	}
449	if (nentries[1] != 0) {
450	glmass_ = glmass[1][il];
451	lspmass_ = lspmass[1][il];
452	effsr1_ = eff[1][il];
453	errsr1_ = eff[1][il];
454	obslimsr1_ = obslim[1][il];
455	explimsr1_ = explim[1][il];
456	explimplussr1_ = explimplus[1][il];
457	explimminussr1_ = explimminus[1][il];
458	float temp = explimsr1_/effsr1_;
459	if (temp < best) {
460	best = temp;
461	ibest = 1;
462	}
463	}
464	if (nentries[2] != 0) {
465	glmass_ = glmass[2][il];
466	lspmass_ = lspmass[2][il];
467	effsr2_ = eff[2][il];
468	errsr2_ = eff[2][il];
469	obslimsr2_ = obslim[2][il];
470	explimsr2_ = explim[2][il];
471	explimplussr2_ = explimplus[2][il];
472	explimminussr2_ = explimminus[2][il];
473	float temp = explimsr2_/effsr2_;
474	if (temp < best) {
475	best = temp;
476	ibest = 2;
477	}
478	}
479	if (nentries[3] != 0) {
480	glmass_ = glmass[3][il];
481	lspmass_ = lspmass[3][il];
482	effsr3_ = eff[3][il];
483	errsr3_ = eff[3][il];
484	obslimsr3_ = obslim[3][il];
485	explimsr3_ = explim[3][il];
486	explimplussr3_ = explimplus[3][il];
487	explimminussr3_ = explimminus[3][il];
488	float temp = explimsr3_/effsr3_;
489	if (temp < best) {
490	best = temp;
491	ibest = 3;
492	}
493	}
494	if (nentries[4] != 0) {
495	glmass_ = glmass[4][il];
496	lspmass_ = lspmass[4][il];
497	effsr4_ = eff[4][il];
498	errsr4_ = eff[4][il];
499	obslimsr4_ = obslim[4][il];
500	explimsr4_ = explim[4][il];
501	explimplussr4_ = explimplus[4][il];
502	explimminussr4_ = explimminus[4][il];
503	float temp = explimsr4_/effsr4_;
504	if (temp < best) {
505	best = temp;
506	ibest = 4;
507	}
508	}
509	if (nentries[5] != 0) {
510	glmass_ = glmass[5][il];
511	lspmass_ = lspmass[5][il];
512	effsr5_ = eff[5][il];
513	errsr5_ = eff[5][il];
514	obslimsr5_ = obslim[5][il];
515	explimsr5_ = explim[5][il];
516	explimplussr5_ = explimplus[5][il];
517	explimminussr5_ = explimminus[5][il];
518	float temp = explimsr5_/effsr5_;
519	if (temp < best) {
520	best = temp;
521	ibest = 5;
522	}
523	}
524	if (nentries[6] != 0) {
525	glmass_ = glmass[6][il];
526	lspmass_ = lspmass[6][il];
527	effsr6_ = eff[6][il];
528	errsr6_ = eff[6][il];
529	obslimsr6_ = obslim[6][il];
530	explimsr6_ = explim[6][il];
531	explimplussr6_ = explimplus[6][il];
532	explimminussr6_ = explimminus[6][il];
533	float temp = explimsr6_/effsr6_;
534	if (temp < best) {
535	best = temp;
536	ibest = 6;
537	}
538	}
539	if (nentries[7] != 0) {
540	glmass_ = glmass[7][il];
541	lspmass_ = lspmass[7][il];
542	effsr7_ = eff[7][il];
543	errsr7_ = eff[7][il];
544	obslimsr7_ = obslim[7][il];
545	explimsr7_ = explim[7][il];
546	explimplussr7_ = explimplus[7][il];
547	explimminussr7_ = explimminus[7][il];
548	float temp = explimsr7_/effsr7_;
549	if (temp < best) {
550	best = temp;
551	ibest = 7;
552	}
553	}
554	if (nentries[8] != 0) {
555	glmass_ = glmass[8][il];
556	lspmass_ = lspmass[8][il];
557	effsr8_ = eff[8][il];
558	errsr8_ = eff[8][il];
559	obslimsr8_ = obslim[8][il];
560	explimsr8_ = explim[8][il];
561	explimplussr8_ = explimplus[8][il];
562	explimminussr8_ = explimminus[8][il];
563	float temp = explimsr8_/effsr8_;
564	if (temp < best) {
565	best = temp;
566	ibest = 8;
567	}
568	}
569
570	bestsr_ = ibest;
571	effsrb_ = eff[ibest][il];
572	errsrb_ = eff[ibest][il];
573	obslimsrb_ = obslim[ibest][il];
574	explimsrb_ = explim[ibest][il];
575	explimplussrb_ = explimplus[ibest][il];
576	explimminussrb_ = explimminus[ibest][il];
577
578	// Fill the cross-section
579	bool done = false;
580	for (int xs=0; xs<nxsec-1; xs++) {
581	if ( (glmass_ >= xsec_mass[xs] && glmass_ < xsec_mass[xs+1]) \|\|
582	(glmass_ < xsec_mass[xs] && xs==0) ) {
583
584	float dd1 = glmass_ - xsec_mass[xs];
585	float dd2 = xsec_mass[xs+1] - xsec_mass[xs];
586	xsec_ = xsec[xs] + (xsec[xs+1] - xsec[xs]) *(dd1/dd2);
587	xsecup_ = xsecup[xs] + (xsecup[xs+1] - xsecup[xs]) *(dd1/dd2);
588	xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs+1] - xsecdwn[xs])*(dd1/dd2);
589	done = true;
590	break;
591	}
592	}
593	if (!done) {
594	int xs = nxsec - 1;
595	float dd1 = glmass_ - xsec_mass[xs];
596	float dd2 = xsec_mass[xs] - xsec_mass[xs-1];
597	xsec_ = xsec[xs] + (xsec[xs] - xsec[xs-1]) *(dd1/dd2);
598	xsecup_ = xsecup[xs] + (xsecup[xs] - xsecup[xs-1]) *(dd1/dd2);
599	xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs] - xsecdwn[xs-1])*(dd1/dd2);
600	}
601
602
603	// Now fill the ntuple
604	babyTree_->Fill();
605
606	}
607	// Now close the file
608
609	babyFile_->cd();
610	babyTree_->Write();
611	babyFile_->Close();
612
613
614	}
615