Plotting/Modules/GeneralToolBox.C

#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <vector>
#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <limits>

#include <TFile.h>
#include <TTree.h>
#include <TCut.h>
#include <TLegend.h>
#include <TLatex.h>
#include <TText.h>
#include <TGraph.h>
#include <TH1.h>
#include <TF1.h>
#include <TMath.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TError.h>
#include <TVirtualPad.h>
#include <TGraphAsymmErrors.h>
#include <TPaveText.h>
#include <TRandom.h>
#ifndef Verbosity
#define Verbosity 0
#endif

/*
#ifndef SampleClassLoaded
#include "SampleClass.C"
#endif
*/
#define GeneralToolBoxLoaded

using namespace std;

bool dopng=false;
bool doC=false;
bool doeps=false;
bool dopdf=false;
string basedirectory="";

TLegend* make_legend(string title,float posx,float posy);
TText* write_title(string title);
TText* write_title_low(string title);

TText* write_text(float xpos,float ypos,string title);
float computeRatioError(float a, float da, float b, float db);
float computeProductError(float a, float da, float b, float db);
TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning);
void setlumi(float l);
void CompleteSave(TCanvas *can, string filename, bool feedback);
void CompleteSave(TVirtualPad *can, string filename, bool feedback);
void write_warning(string funcname, string text);
void write_error(string funcname, string text);
void write_info(string funcname, string text);
void DrawPrelim();
string get_directory();
//-------------------------------------------------------------------------------------
float lumi;

template<typename U>
inline bool isanyinf(U value)
{
return !(value >= std::numeric_limits<U>::min() && value <=
std::numeric_limits<U>::max());
}

template<class A>
string any2string(const A& a){
    ostringstream out;
    out << a;
    return out.str();
}

void do_png(bool s) { dopng=s;}
void do_eps(bool s) { doeps=s;}
void do_C(bool s)   { doC=s;}
void do_pdf(bool s)   { dopdf=s;}

string topdir(string child) {
    string tempdirectory=child;
    if(tempdirectory.substr(tempdirectory.length()-1,1)=="/") tempdirectory=tempdirectory.substr(0,tempdirectory.length());
    //we now have a directory without the trailing slash so we can just look for the last non-slash character :-) 
    for(int ichar=tempdirectory.length()-1;ichar>=0;ichar--) {
      if(tempdirectory.substr(ichar,1)=="/") {
        return tempdirectory.substr(0,ichar);
      }
    }
}

template < typename CHAR_TYPE,
           typename TRAITS_TYPE = std::char_traits<CHAR_TYPE> >

struct basic_teebuf : public std::basic_streambuf< CHAR_TYPE, TRAITS_TYPE >
{
    typedef std::basic_streambuf< CHAR_TYPE, TRAITS_TYPE > streambuf_type ;
    typedef typename TRAITS_TYPE::int_type int_type ;

    basic_teebuf( streambuf_type* buff_a, streambuf_type* buff_b )
            : first(buff_a), second(buff_b) {}

    protected:
        virtual int_type overflow( int_type c )
        {
            const int_type eof = TRAITS_TYPE::eof() ;
            if( TRAITS_TYPE::eq_int_type( c, eof ) )
                return TRAITS_TYPE::not_eof(c) ;
            else
            {
                const CHAR_TYPE ch = TRAITS_TYPE::to_char_type(c) ;
                if( TRAITS_TYPE::eq_int_type( first->sputc(ch), eof ) ||
                    TRAITS_TYPE::eq_int_type( second->sputc(ch), eof ) )
                        return eof ;
                else return c ;
            }
        }

        virtual int sync()
        { return !first->pubsync() && !second->pubsync() ? 0 : -1 ; }

    private:
        streambuf_type* first ;
        streambuf_type* second ;
};

template < typename CHAR_TYPE,
           typename TRAITS_TYPE = std::char_traits<CHAR_TYPE> >
struct basic_teestream : public std::basic_ostream< CHAR_TYPE, TRAITS_TYPE >
{
    typedef std::basic_ostream< CHAR_TYPE, TRAITS_TYPE > stream_type ;
    typedef basic_teebuf< CHAR_TYPE, TRAITS_TYPE > streambuff_type ;

    basic_teestream( stream_type& first, stream_type& second )
         : stream_type( &stmbuf), stmbuf( first.rdbuf(), second.rdbuf() ) {}

    ~basic_teestream() { stmbuf.pubsync() ; }

    private: streambuff_type stmbuf ;
};

typedef basic_teebuf<char> teebuf ;
typedef basic_teestream<char> teestream ;

std::ofstream file("LOG.txt",ios::app) ;
std::ofstream efile("LOGerr.txt",ios::app) ;
teestream dout( file, std::cout ) ; // double out
teestream eout( efile, std::cout ) ; // double out (errors)

void ensure_directory_exists(string thisdirectory) {
  struct stat st;
  if(stat(thisdirectory.c_str(),&st) == 0) {
    if(Verbosity>0) dout << "Directory " << thisdirectory << " exists!" << endl;
  }
  else {
    if(Verbosity>0) dout << "Directory " << thisdirectory << " does not exist. Need to create it!" << endl;
    ensure_directory_exists(topdir(thisdirectory));
    if (mkdir(thisdirectory.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH))
            if(Verbosity>0) dout << "Created the directory " << thisdirectory << endl;
  }
}

void initialize_log() {
  dout << "____________________________________________________________" << endl;
  dout << endl;
  dout << "                                                            " << endl;
  dout << "          JJJJJJJJJJJZZZZZZZZZZZZZZZZZZZBBBBBBBBBBBBBBBBB   " << endl;
  dout << "          J:::::::::JZ:::::::::::::::::ZB::::::::::::::::B  " << endl;
  dout << "          J:::::::::JZ:::::::::::::::::ZB::::::BBBBBB:::::B " << endl;
  dout << "          JJ:::::::JJZ:::ZZZZZZZZ:::::Z BB:::::B     B:::::B" << endl;
  dout << "            J:::::J  ZZZZZ     Z:::::Z    B::::B     B:::::B" << endl;
  dout << "            J:::::J          Z:::::Z      B::::B     B:::::B" << endl;
  dout << "            J:::::J         Z:::::Z       B::::BBBBBB:::::B " << endl;
  dout << "            J:::::j        Z:::::Z        B:::::::::::::BB  " << endl;
  dout << "            J:::::J       Z:::::Z         B::::BBBBBB:::::B " << endl;
  dout << "JJJJJJJ     J:::::J      Z:::::Z          B::::B     B:::::B" << endl;
  dout << "J:::::J     J:::::J     Z:::::Z           B::::B     B:::::B" << endl;
  dout << "J::::::J   J::::::J  ZZZ:::::Z     ZZZZZ  B::::B     B:::::B" << endl;
  dout << "J:::::::JJJ:::::::J  Z::::::ZZZZZZZZ:::ZBB:::::BBBBBB::::::B" << endl;
  dout << " JJ:::::::::::::JJ   Z:::::::::::::::::ZB:::::::::::::::::B " << endl;
  dout << "   JJ:::::::::JJ     Z:::::::::::::::::ZB::::::::::::::::B  " << endl;
  dout << "     JJJJJJJJJ       ZZZZZZZZZZZZZZZZZZZBBBBBBBBBBBBBBBBB   " << endl;
  dout << "                                                            " << endl;
  dout << endl << endl;
  dout << "____________________________________________________________" << endl;
  time_t rawtime;
  struct tm * timeinfo;
  time ( &rawtime );
  dout << "       Analysis run on " << asctime (localtime ( &rawtime ));
  dout << "____________________________________________________________" << endl;
  dout << " Results saved in : " << get_directory() << endl << endl;
}

void set_directory(string basedir="") {
  if(basedir.substr(0,1)=="/") basedir=basedir.substr(1,basedir.length()-1);
  if(basedir.substr(basedir.length()-1,1)!="/") basedir+="/";
  char currentpath[1024];
  getcwd(currentpath,1024);
  basedirectory=(string)currentpath+"/Plots/"+basedir;
  ensure_directory_exists(basedirectory);
  initialize_log();
}

string get_directory() {
  return basedirectory;
}

string extract_directory(string savethis) {
  bool foundslash=false;
  int position=savethis.length();
  while(!foundslash&&position>0) {
    position--;
    if(savethis.substr(position,1)=="/") foundslash=true;
  }
  if(position>0) return savethis.substr(0,position+1);
  else return "";
}

void CompleteSave(TCanvas *can, string filename, bool feedback=false, bool redraw=true) {
  //any change you make here should also be done below in the CompleteSave function for virtual pads
  Int_t currlevel=gErrorIgnoreLevel;
  if(!feedback) gErrorIgnoreLevel=1001;
  if(redraw) can->RedrawAxis();
  ensure_directory_exists(extract_directory(basedirectory+filename));
  if(dopng) can->SaveAs((basedirectory+filename+".png").c_str());
  if(doeps) can->SaveAs((basedirectory+filename+".eps").c_str());
  if(dopdf) can->SaveAs((basedirectory+filename+".pdf").c_str());
  if(doC) can->SaveAs((basedirectory+filename+".C").c_str());
  gErrorIgnoreLevel=currlevel;
  dout << "Saved " << filename << " in all requested formats" << endl;
}

void CompleteSave(TVirtualPad *can, string filename, bool feedback=false, bool redraw=true) {
  Int_t currlevel=gErrorIgnoreLevel;
  if(!feedback) gErrorIgnoreLevel=1001;
  if(redraw) can->RedrawAxis();
  ensure_directory_exists(extract_directory(basedirectory+filename));
  if(dopng) can->SaveAs((basedirectory+filename+".png").c_str());
  if(doeps) can->SaveAs((basedirectory+filename+".eps").c_str());
  if(doC) can->SaveAs((basedirectory+filename+".C").c_str());
  gErrorIgnoreLevel=currlevel;
  dout << "Saved " << filename << " in all requested formats" << endl;
}


void setlumi(float l) {
lumi=l;
}

int write_first_line(vector<vector<string> > &entries) {
  if(entries.size()>0) {
    vector<string> firstline = entries[0];
    int ncolumns=firstline.size();
    int ndividers=ncolumns+1;
    int cellwidth=(int)(((float)(60-ndividers))/(ncolumns));
    dout << " |";
    for(int idiv=0;idiv<ncolumns;idiv++) {
      for(int isig=0;isig<cellwidth;isig++) dout << "-";
      dout << "|";
    }
    dout << endl;
    return ncolumns;
  } else {
    return 0;
  }
}

void write_entry(string entry,int width,int iline=0,int ientry=0) {
  int currwidth=entry.size();
  while(currwidth<width) {
    entry=" "+entry;
    if(entry.size()<width) entry=entry+" ";
    currwidth=entry.size();
  }
  bool do_special=false;
  if(iline==1&&ientry==1) {  dout << "\033[1;32m" << entry << "\033[0m|";do_special=true;}//observed
  if(iline==1&&ientry==2) {  dout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
  if(iline==2&&ientry==1) {  dout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
  if(iline==2&&ientry==2) {  dout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
  if(!do_special) dout << entry << "|";
}

void make_nice_table(vector<vector <string> > &entries) {
  int ncolumns=write_first_line(entries);
  int cellwidth=(int)(((float)(60-(ncolumns+1)))/(ncolumns));
  for(int iline=0;iline<entries.size();iline++) {
    vector<string> currline = entries[iline];
    dout << " |";
    for(int ientry=0;ientry<currline.size();ientry++) {
      write_entry(currline[ientry],cellwidth);
    }
    dout << endl;
    if(iline==0)   write_first_line(entries);
  }
  write_first_line(entries);
}

void make_nice_jzb_table(vector<vector <string> > &entries) {
  int ncolumns=write_first_line(entries);
  int cellwidth=(int)(((float)(60-(ncolumns+1)))/(ncolumns));
  for(int iline=0;iline<entries.size();iline++) {
    vector<string> currline = entries[iline];
    dout << " |";
    for(int ientry=0;ientry<currline.size();ientry++) {
      write_entry(currline[ientry],cellwidth,iline,ientry);
    }
    dout << endl;
    if(iline==0)   write_first_line(entries);
  }
  write_first_line(entries);
}


void write_warning(string funcname, string text) {
  eout << endl << endl;  
  eout << "\033[1;33m" << "                          _             " << endl;
  eout << "\033[1;33m" << "                         (_)            " << endl;
  eout << "\033[1;33m" << "__      ____ _ _ __ _ __  _ _ __   __ _ " << endl;
  eout << "\033[1;33m" << "\\ \\ /\\ / / _` | '__| '_ \\| | '_ \\ / _` |" << endl;
  eout << "\033[1;33m" << " \\ V  V / (_| | |  | | | | | | | | (_| |" << endl;
  eout << "\033[1;33m" << "  \\_/\\_/ \\__,_|_|  |_| |_|_|_| |_|\\__, |" << endl;
  eout << "\033[1;33m" << "                                   __/ |" << endl;
  eout << "\033[1;33m" << "                                  |___/ " << endl;
  eout << endl;
  eout << "\033[1;33m [" << funcname << "] " << text << " \033[0m" << endl;
  eout << endl << endl;
}
void write_error(string funcname, string text) {
  eout << endl << endl;
  eout << "\033[1;31m  ___ _ __ _ __ ___  _ __ " << endl;
  eout << "\033[1;31m / _ \\  __|  __/ _ \\| '__|" << endl;
  eout << "\033[1;31m|  __/ |  | | | (_) | |   " << endl;
  eout << "\033[1;31m \\___|_|  |_|  \\___/|_|   " << endl;
  eout << endl;
  eout << "\033[1;31m [" << funcname << "] " << text << " \033[0m" << endl;
  eout << endl << endl;
}

void write_info(string funcname, string text) {
  dout << endl << endl;
  dout << "\033[1;34m   _____        __        " << endl;
  dout << "\033[1;34m  |_   _|      / _|      " << endl;
  dout << "\033[1;34m    | |  _ __ | |_ ___   " << endl;
  dout << "\033[1;34m    | | | '_ \\|  _/ _ \\  " << endl;
  dout << "\033[1;34m   _| |_| | | | || (_) | " << endl;
  dout << "\033[1;34m  |_____|_| |_|_| \\___/  " << endl;
  dout << endl;
  dout << "\033[1;34m [" << funcname << "] " << text << " \033[0m" << endl;
  dout << endl << endl;
}

TText* write_text(float xpos,float ypos,string title)
{
   TLatex* titlebox = new TLatex (xpos,ypos,title.c_str());
   titlebox->SetNDC(true);
   titlebox->SetTextFont(42);
   titlebox->SetTextSize(0.04);
   titlebox->SetTextAlign(21);
   return titlebox;
}

TText* write_title(string title)
{
   TText* titlebox = write_text(0.5,0.945,title);
   return titlebox;
}

TText* write_cut_on_canvas(string cut) {
//  TLatex *normbox = new TLatex(0.96,0.5,cut.c_str());
  TLatex *normbox = new TLatex(0.96,0.5,"");//currently deactivated
  normbox->SetNDC(true);
  normbox->SetTextFont(42);
  normbox->SetTextSize(0.01);
  normbox->SetTextAlign(21);
  normbox->SetTextAngle(270);
  return normbox;
}

TText* write_title_low(string title)
{
   TText* titlebox = write_text(0.5,0.94,title);
   return titlebox;
}

TLegend* make_legend(string title="", float posx=0.6, float posy=0.55)
{
//  TLegend *leg = new TLegend(0.65,0.65,0.89,0.89);
  gStyle->SetTextFont(42);
  TLegend *leg = new TLegend(posx,posy,0.89,0.89);
  if(title!="") leg->SetHeader(title.c_str());
  leg->SetTextFont(42);
  leg->SetTextSize(0.04);
  leg->SetFillColor(kWhite);
  leg->SetBorderSize(0);
  leg->SetLineColor(kWhite);
  DrawPrelim();
/*  TText *writeline1 = write_text(0.77,0.87,"CMS Preliminary 2011");
  stringstream lumitext;
  lumitext<<"#sqrt{s}=7, L = "<<lumi<<" pb^{-1}";
  TText *writeline2 = write_text(0.77,0.83,lumitext.str().c_str());
  writeline1->SetTextSize(0.03);
  writeline2->SetTextSize(0.03);
  writeline1->Draw();
  writeline2->Draw();
  */
  return leg;
}

TGraph* make_nice_ratio(int nbins,float binning[],TH1F* histo)
{
  float errorsquared[nbins];
  float errors[nbins];
  float bincontent[nbins];
  for (int i=0;i<nbins;i++) {
    errorsquared[i]=0;
    bincontent[i]=0;
    errors[i]=0;
  }
  float currlimit=binning[0];
  int currtoplim=1;
  for(int ibin=1;ibin<=histo->GetNbinsX();ibin++)
  {
    if(binning[currtoplim]<histo->GetBinCenter(ibin)) currtoplim++;
    dout << "Bin i=" << ibin << " with bin center " << histo->GetBinCenter(ibin) << " contains " << histo->GetBinContent(ibin) << " is within " << binning[currtoplim-1] << " and " << binning[currtoplim] << endl;
    
  }
    
  return 0;
}
  
float statErrorN(float x){return x - 0.5*TMath::ChisquareQuantile(0.3173/2,2*x);}
float statErrorP(float x){return 0.5*TMath::ChisquareQuantile(1-0.3173/2,2*(x+1))-x;}
float lowLimit(float a, float x){return 0.5*TMath::ChisquareQuantile(a,2*x);}
float highLimit(float a,float x){return 0.5*TMath::ChisquareQuantile(1-a,2*(x+1));}

float computeRatioError(float a, float da, float b, float db)
{
    float val=0.;
    float errorSquare = (a/b)*(a/b)*( (da/a)*(da/a) + (db/b)*(db/b));
    val = TMath::Sqrt(errorSquare);
    return val;

}
float computeProductError(float a, float da, float b, float db)
{
    float val=0.;
    float errorSquare = (a*b)*(a*b)*( (da/a)*(da/a) + (db/b)*(db/b));
    val = TMath::Sqrt(errorSquare);
    return val;
}

TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning)
{
        int absJZBbinsNumber = binning.size()-1;
        TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
        
        for(unsigned int i=0;i<absJZBbinsNumber;i++)
        {
                float xCenter=h1->GetBinCenter(i+1);
                float xWidth=(h1->GetBinWidth(i+1))*0.5;
                float nominatorError = h1->GetBinError(i+1);
                float nominator=h1->GetBinContent(i+1);
                float denominatorError=h2->GetBinError(i+1);
                float denominator=h2->GetBinContent(i+1);
                float errorN = 0;
                float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
                if(id==1) // (is data)
                {
                        errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
                        errorN = errorP; // symmetrize using statErrorP
                } else {
                        errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
                        errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
                }
                if(denominator!=0) {
                        graph->SetPoint(i, xCenter, nominator/denominator);
                        graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
                }
                else {
                        graph->SetPoint(i, xCenter, -999);
                        graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
                }
        }       
        return graph;
}

string print_range(float cent, float down, float up) {//note that up&down can be flipped, we don't care, but the central value needs to come 1st!
 float uperr=0,downerr=0;
 if(down>up&&down>cent) uperr=down-cent;
 if(up>down&&up>cent) uperr=up-cent;
 if(down<cent&&down<up) downerr=cent-down;
 if(up<cent&&up<down) downerr=cent-up;
 if(cent>up&&cent>down&&(up!=0&&down!=0)) write_error("print_range"," WATCH OUT: THE CENTRAL VALUE SEEMS TO BE LARGER THAN BOTH UP&DOWN!");
 if(cent<up&&cent<down&&(up!=0&&down!=0)) write_error("print_range"," WATCH OUT: THE CENTRAL VALUE SEEMS TO BE SMALLER THAN BOTH UP&DOWN!");
 stringstream result;
 result << cent << " + " << uperr << " - " << downerr;
 return result.str();
}

void bubbleSort ( int arr [ ], int size, int order [ ]) // nice way to sort an array (called arr) which is currently in a random order (indices in (order")
{
   int last = size - 2;
   int isChanged = 1;

   while ( last >= 0 && isChanged )
   {
           isChanged = 0;
           for ( int k = 0; k <= last; k++ )
               if ( arr[k] > arr[k+1] )
               {
                   swap ( arr[k], arr[k+1] );
                   isChanged = 1;
                   int bkp=order[k];
                   order[k]=order[k+1];
                   order[k+1]=bkp;
               }
           last--;
   }
}

void swapvec(vector<float> &vec,int j, int k) {
float bkp=vec[j];
vec[j]=vec[k];
vec[k]=bkp;
}

void bubbleSort ( vector<float> &arr , vector<int> &order) // nice way to sort an array (called arr) which is currently in a random order (indices in (order")
{
   int last = arr.size() - 2;
   int isChanged = 1;

   while ( last >= 0 && isChanged )
   {
           isChanged = 0;
           for ( int k = 0; k <= last; k++ )
               if ( arr[k] > arr[k+1] )
               {
                   swapvec (arr,k,k+1);
                   isChanged = 1;
                   int bkp=order[k];
                   order[k]=order[k+1];
                   order[k+1]=bkp;
               }
           last--;
   }
}

int numerichistoname=0;
bool givingnumber=false;
string GetNumericHistoName() {
  while(givingnumber) sleep(1);
  givingnumber=true;
  stringstream b;
  b << "h_" << numerichistoname;
  numerichistoname++;
  givingnumber=false;
  return b.str();
}

//********************** BELOW : CUT INTERPRETATION **************************//
void splitupcut(string incut, vector<string> &partvector)
{
    //idea: go thru the string called incut; if a parantheses is opened, then the cut cannot be split up until the parantheses is closed.
    //ok anyway screw the parantheses.
    int paranthesis_open=0;
    int substr_start=0;
    string currchar="";
    for (int ichar=0;ichar<incut.length();ichar++)
    {
        currchar=incut.substr(ichar,1);
//        if(currchar=="(") paranthesis_open++;
//        if(currchar==")") paranthesis_open--;
        if(currchar=="&"&&incut.substr(ichar+1,1)=="&"&&paranthesis_open==0) {
            partvector.push_back(incut.substr(substr_start,ichar-substr_start));
            substr_start=ichar+2;
        }
    }
    partvector.push_back(incut.substr(substr_start,incut.length()-substr_start));
    if(Verbosity>1) {
        dout << "[ splitupcut() ] : The cut vector now contains the following elements: "<< endl;
        for (int ipart=0;ipart<partvector.size();ipart++)
        {
            dout << "     - " << partvector[ipart] << endl;
        }
    }
}

int atleastvalue(string expression, int &morethanlessthan) // takes in an expression such as ">2" or ">=3" and returns e.g. 3 (in both examples) 
{
  int retval=0;
  if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
//    dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
    morethanlessthan=1;
    return atoi(expression.substr(2,1).c_str());
  }
  if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
//    dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
    morethanlessthan=0;
    return atoi(expression.substr(1,1).c_str());
  }
  if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
//    dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
    morethanlessthan=-1;
    return 1+atoi(expression.substr(1,1).c_str());
  }
  if(expression.substr(0,1)==">") {
//    dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
    morethanlessthan=1;
    return 1+atoi(expression.substr(2,1).c_str());
  }
  if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
//    dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
    morethanlessthan=-1;
    return 1+atoi(expression.substr(2,1).c_str());
  }
}

int do_jet_cut(string incut, int *nJets) {
  string expression=(incut.substr(12,incut.size()-12));
  dout << "Going to analyze the jet cut : " << expression << " with 0,1 being " << expression.substr(0,1) << " and 1,1 being " << expression.substr(1,1) << endl;
  if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
    int nJet=atoi(expression.substr(2,1).c_str());
    for(int i=nJet+1;i<20;i++) nJets[i]=0;
    dout << "Is of type <=" << endl;
    return 0;
  }
  if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
    int nJet=atoi(expression.substr(2,1).c_str());
    for(int i=0;i<20&&i!=nJet;i++) nJets[i]=0;
    dout << "Is of type ==" << endl;
    return 0;
  }
  if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
    int nJet=atoi(expression.substr(2,1).c_str());
    for(int i=0;i<nJet&&i!=nJet;i++) nJets[i]=0;
    dout << "Is of type >=" << endl;
    return 0;
  }
  if(expression.substr(0,1)=="<") {
    int nJet=atoi(expression.substr(1,1).c_str());
    for(int i=nJet;i<20;i++) nJets[i]=0;
    dout << "Is of type <" << endl;
    return 0;
  }
  if(expression.substr(0,1)==">") {
    int nJet=atoi(expression.substr(1,1).c_str());
    for(int i=0;i<nJet+1&&i!=nJet;i++) nJets[i]=0;
    dout << "Is of type >" << endl;
    return 0;
  }
}

string interpret_cut(string incut, bool &isJetCut, int *permittednJets)
{
//    isJetCut=false;nJets=-1;
    if(incut=="()") return "";
    while(incut.substr(0,1)=="(") incut=incut.substr(1,incut.length()-1);
    while(incut.length()>0&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(0,incut.length()-1);
//    if(incut.substr(0,1)=="("&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(1,incut.length()-2); //this is to make (cut) to cut.
   
    if(Verbosity>0) {
        dout << "Now interpreting cut " << incut << endl;
    }
    /*
    if(incut=="ch1*ch2<0") return "OS";
    if(incut=="id1==id2") return "SF";
    if(incut=="id1!=id2") return "OF";
    */
    if(incut=="ch1*ch2<0") return "";
    if(incut=="(mll>55&&mll<70)||(mll>112&&mll<160)") return "SB";
    if(incut=="(mll>61&&mll<70)||(mll>112&&mll<190)") return "SB'";
    if(incut=="id1==id2") return "";
    if(incut=="id1!=id2") return "";
    if(incut=="mll>2") return "";

    if(incut=="mll>0") return ""; // my typical "fake cut"

    if(incut=="passed_triggers||!is_data") return "Triggers";
    if(incut=="pfjzb[0]>-998") return "";

   
    if(incut=="id1==0") return "ee";
    if(incut=="id1==1") return "#mu#mu";
    if(incut=="abs(mll-91.2)<20") return "|m_{l^{+}l^{-}}-m_{Z}|<20";
    if(incut=="pfJetGoodID[0]") return "";
    if(incut=="pfJetGoodID[1]") return "";
    if((int)incut.find("pfJetGoodNum")>-1) {
        //do_jet_cut(incut,permittednJets);
        stringstream result;
        result << "nJets" << incut.substr(12,incut.size()-12);
/*        dout << "Dealing with a jet cut: " << incut << endl;
        stringstream result;
        result << "nJets" << incut.substr(12,incut.size()-12);
        isJetCut=true;
        if(exactjetcut(incut,nJets))
//      nJets=atleastvalue((incut.substr(12,incut.size()-12)),morethanlessthan);
        return result.str();*/
        return result.str();
    }
    return incut;
}

string interpret_nJet_range(int *nJets) {
  for (int i=0;i<20;i++) dout << i << " : " << nJets[i] << endl;
  return "hello";
}

string interpret_cuts(vector<string> &cutparts)
{
    stringstream nicecut;
    int nJets;
    bool isJetCut;
    int finalJetCut=-1;
    int permittednJets[20];
    for(int ijet=0;ijet<20;ijet++) permittednJets[ijet]=1;
    int morethanlessthan=0;//-1: less than, 0: exactly, 1: more than
    for(int icut=0;icut<cutparts.size();icut++)
    {
        if(icut==0) nicecut<<interpret_cut(cutparts[icut],isJetCut,permittednJets);
        else {
          string nice_this_cut = interpret_cut(cutparts[icut],isJetCut,permittednJets);//blublu
          if(nice_this_cut.length()>0&&nicecut.str().length()>0) {
            if(!isJetCut) nicecut<<"&&"<<nice_this_cut;
            else {
              if(nJets>finalJetCut) finalJetCut=nJets;
          }
          }
          if(nice_this_cut.length()>0&&nicecut.str().length()==0) {
            if(!isJetCut) {
              nicecut<<nice_this_cut;
            }
            else {
              if(nJets>finalJetCut) finalJetCut=nJets;
            }
          }
        }
    }
    if(finalJetCut>-1) {
      if(nicecut.str().length()==0) {
        nicecut << "nJets#geq" << finalJetCut;
      }
      else
      {
        nicecut << "&&nJets#geq " << finalJetCut;
      }
    }
    
//    dout << "The nJet allowed range is given by: " << interpret_nJet_range(permittednJets) << endl;
    
    return nicecut.str();
}

string decipher_cut(TCut originalcut,TCut ignorethispart)
{
    string incut=(const char*)originalcut;
    string ignore=(const char*)ignorethispart;
    
    if(ignore.length()>0 && incut.find(ignore)!=string::npos) incut=incut.replace(incut.find(ignore),ignore.length(),"");
    
    vector<string>cutparts;
    splitupcut(incut,cutparts);
    string write_cut=interpret_cuts(cutparts);
    return write_cut;
}

//********************** ABOVE : CUT INTERPRETATION **************************//

Double_t GausRandom(Double_t mu, Double_t sigma) {
  return gRandom->Gaus(mu,sigma);// real deal
  //return mu;//debugging : no smearing.
}

int functionalhistocounter=0;
TH1F * makehistofromfunction(TF1 *f1,TH1F *model) {
  TH1F *histo = (TH1F*)model->Clone();
  functionalhistocounter++;
  stringstream histoname;
  histoname << "histo_based_on_function_" << f1->GetName() << "__"<<functionalhistocounter;
  histo->SetTitle(histoname.str().c_str());
  histo->SetName(histoname.str().c_str());
  int nbins=histo->GetNbinsX();
  float low=histo->GetBinLowEdge(1);
  float hi=histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX());
  
  for(int i=0;i<=nbins;i++) {
    histo->SetBinContent(i,(f1->Integral(histo->GetBinLowEdge(i),histo->GetBinLowEdge(i)+histo->GetBinWidth(i)))/histo->GetBinWidth(i));
    histo->SetBinError(i,TMath::Sqrt(histo->GetBinContent(i)));
  }

  return histo;
}

float hintegral(TH1 *histo, float low, float high) {
  float sum=0;
  for(int i=1;i<histo->GetNbinsX();i++) {
    if((histo->GetBinLowEdge(i)>=low)&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))<=high) sum+=histo->GetBinContent(i);
    //now on to the less clear cases!
    if(histo->GetBinLowEdge(i)<low&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))>low) {
      //need to consider this case still ... the bin is kind of in range but not sooooo much.
    }
    if(histo->GetBinLowEdge(i)<high&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))>high) {
      //need to consider this case still ... the bin is kind of in range but not sooooo much.
    }
    
  }
  return sum;
}

void DrawPrelim_OLD() {
  TText *writeline1 = write_text(0.77,0.87,"CMS Preliminary 2011");
  stringstream lumitext;
  lumitext<<"#sqrt{s}=7, L="<< lumi/1000<<" fb^{-1}";
  TText *writeline2 = write_text(0.77,0.83,lumitext.str().c_str());
  writeline1->SetTextSize(0.03);
  writeline2->SetTextSize(0.03);
  writeline1->Draw();
  writeline2->Draw();
}

void DrawPrelim() {
  string barn="pb";
  float writelumi=lumi;
  if(writelumi>=1000)
  {
    writelumi/=1000;
    barn="fb";
  }
  
  stringstream prelimtext;
  //prelimtext << "CMS Preliminary 2011 , #sqrt{s}= 7 TeV, L= O(1) fb^{-1}"; //temporary replacement
  prelimtext << "CMS Preliminary 2011 , #sqrt{s}= 7 TeV, L= "
             << std::setprecision(2) <<writelumi<<" "<<barn<<"^{-1}";
  TPaveText *eventSelectionPaveText = new TPaveText(0.27, 0.93,0.77, 1.0,"blNDC");
  eventSelectionPaveText->SetFillStyle(4000);
  eventSelectionPaveText->SetBorderSize(0.1);
  eventSelectionPaveText->SetFillColor(kWhite);
  eventSelectionPaveText->SetTextFont(42);
  eventSelectionPaveText->SetTextSize(0.0415);
  eventSelectionPaveText->AddText(prelimtext.str().c_str());
  eventSelectionPaveText->Draw();
}

string newjzbexpression(string oldexpression,float shift) {
  stringstream ss;
  if(shift>0) ss<<"("<<oldexpression<<"+"<<shift<<")";
  if(shift<0) ss<<"("<<oldexpression<<shift<<")";
  if(shift==0) ss<<oldexpression;
  return ss.str();
}
  
double Round(double num, unsigned int dig)
{
    num *= pow(10, dig);
    if (num >= 0)
        num = floor(num + 0.5);
    else
        num = ceil(num - 0.5);
    num/= pow(10, dig);
    return num;
}

// The two functions below are for distributed processing

int get_job_number(float ipoint, float Npoints,float Njobs) {
  float pointposition=(ipoint/Npoints);
  int njob=floor(pointposition*Njobs);
  if(njob>=Njobs) njob--;
//  cout << "Looking at point " << ipoint << " out of " << Npoints << " which is at position " << pointposition << " corresponding to  " << pointposition*Njobs << " --> JOB " << njob << endl;
  return njob;
}
  

bool do_this_point(int ipoint, int Npoints, int jobnumber, int Njobs) {
  if(get_job_number(ipoint,Npoints,Njobs)==jobnumber) return true;
  return false;
}

Double_t  DoIntegral(TH1F *histo, Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2, Double_t & error ,
                          Option_t *option, Bool_t doError)
{
   // internal function compute integral and optionally the error  between the limits
   // specified by the bin number values working for all histograms (1D, 2D and 3D)

   Int_t nbinsx = histo->GetNbinsX();
   if (binx1 < 0) binx1 = 0;
   if (binx2 > nbinsx+1 || binx2 < binx1) binx2 = nbinsx+1;
   if (histo->GetDimension() > 1) {
      Int_t nbinsy = histo->GetNbinsY();
      if (biny1 < 0) biny1 = 0;
      if (biny2 > nbinsy+1 || biny2 < biny1) biny2 = nbinsy+1;
   } else {
      biny1 = 0; biny2 = 0;
   }
   if (histo->GetDimension() > 2) {
      Int_t nbinsz = histo->GetNbinsZ();
      if (binz1 < 0) binz1 = 0;
      if (binz2 > nbinsz+1 || binz2 < binz1) binz2 = nbinsz+1;
   } else {
      binz1 = 0; binz2 = 0;
   }

   //   - Loop on bins in specified range
   TString opt = option;
   opt.ToLower();
   Bool_t width   = kFALSE;
   if (opt.Contains("width")) width = kTRUE;


   Double_t dx = 1.;
   Double_t dy = 1.;
   Double_t dz = 1.;
   Double_t integral = 0;
   Double_t igerr2 = 0;
   for (Int_t binx = binx1; binx <= binx2; ++binx) {
      if (width) dx = histo->GetXaxis()->GetBinWidth(binx);
      for (Int_t biny = biny1; biny <= biny2; ++biny) {
         if (width) dy = histo->GetYaxis()->GetBinWidth(biny);
         for (Int_t binz = binz1; binz <= binz2; ++binz) {
            if (width) dz = histo->GetZaxis()->GetBinWidth(binz);
            Int_t bin  = histo->GetBin(binx, biny, binz);
            if (width) integral += histo->GetBinContent(bin)*dx*dy*dz;
            else       integral += histo->GetBinContent(bin);
            if (doError) {
               if (width)  igerr2 += histo->GetBinError(bin)*histo->GetBinError(bin)*dx*dx*dy*dy*dz*dz;
               else        igerr2 += histo->GetBinError(bin)*histo->GetBinError(bin);
            }
         }
      }
   }

   if (doError) error = TMath::Sqrt(igerr2);
   return integral;
}

Double_t IntegralAndError(TH1F *histo, Int_t binx1, Int_t binx2, Double_t & error, Option_t *option)
{
   //Return integral of bin contents in range [binx1,binx2] and its error
   // By default the integral is computed as the sum of bin contents in the range.
   // if option "width" is specified, the integral is the sum of
   // the bin contents multiplied by the bin width in x.
   // the error is computed using error propagation from the bin errors assumming that
   // all the bins are uncorrelated
   return DoIntegral(histo,binx1,binx2,0,-1,0,-1,error,option,kTRUE);
}
Revision:	1.20
Committed:	Mon Aug 15 14:27:40 2011 UTC (13 years, 8 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.19:	+1 -1 lines
Log Message:	Minimal change to keep the directory from getting too crowded (plots now saved in Plots subdirectory)