Plotting/Modules/SUSYScan.C

#include <iostream>
#include <vector>
#include <sys/stat.h>
#include <fstream>
#include <pthread.h>
#include <assert.h>
#include <cctype>

#include <TCut.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TColor.h>
#include <TPaveText.h>
#include <TRandom.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TSQLResult.h>
#include <TProfile.h>
#include <TKey.h>

#include "ShapeDroplet.C"
#include "ShapeLimit.C"

//#include "TTbar_stuff.C"
using namespace std;

using namespace PlottingSetup;

void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, float njobs, float jobnumber, bool systonly, bool effonly,bool shapeanalysis);

bool isValidShapeResult(ShapeDroplet droplet) {
  if(droplet.isValid==false) return false;
  if(droplet.observed>50000) {
    // error code of ShapeLimit/CreateModel.sh, invoked if there were too many errors during computation
    write_warning(__FUNCTION__,"Shape limit was ridiculously high ("+any2string(droplet.observed)+"); this is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
    return false; 
  }
  if(droplet.observed==-12345) {
        write_warning(__FUNCTION__,"Shape limit was ("+any2string(droplet.observed)+"); this is the error signature of CreateModel.sh / ShapeDroplet. Shape result will be discarded, and c&c will be carried out instead.");
        return false; 
  }
  if(droplet.expected>-1&&droplet.expected==0&&droplet.observed>5000) {
        write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; this is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
        return false; 
  }
  if(droplet.expected>-1&&droplet.expected/droplet.observed>50) {
        write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; such a high ratio between expected & observed limit is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
        return false;
  }
  if(droplet.expected>-1&&droplet.observed/droplet.expected>50) {
        write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; such a high ratio between observed & expected limit is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
        return false;
  }
  return true;
}

void susy_scan_axis_labeling(TH2F *histo) {
  histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}");
  histo->GetXaxis()->CenterTitle();
  histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
  histo->GetYaxis()->CenterTitle();
}

void prepare_scan_axis(TH2 *variablemap,int scantype) {
  //assuming regular SMS
  variablemap->GetXaxis()->SetTitle("m_{glu}");
  variablemap->GetYaxis()->SetTitle("m_{LSP}");
 
  if(scantype==mSUGRA) {
    variablemap->GetXaxis()->SetTitle("m_{0}");
    variablemap->GetYaxis()->SetTitle("m_{1/2}");
  }

  if(scantype==GMSB) {
    variablemap->GetXaxis()->SetTitle("m_{glu}");
    variablemap->GetYaxis()->SetTitle("m_{NLSP}");
  }
 
  variablemap->GetXaxis()->CenterTitle();
  variablemap->GetYaxis()->CenterTitle();
}

void set_SUSY_style() {
  Int_t MyPalette[100];
  Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
  Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
  Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
  Double_t stop[] = {0., .25, .50, .75, 1.0};
  Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
  for (int i=0;i<100;i++) MyPalette[i] = FI+i;
  
  gStyle->SetPalette(100, MyPalette);
 
  gStyle->SetPadRightMargin(0.15);

}

bool delete_any_cached_scans() {
  char hostname[1023];
  gethostname(hostname,1023);
  
  
  /*
  //This should only be called via CRAB
  dout << "   Deleting all cached files." << endl;
  gSystem->Exec("ls -ltrh | grep root");
  if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) 
  {
        vector<string> all_files;
        char currentpath[1024];
        TString directory=TString(getcwd(currentpath,1024));
        directory=directory+"/";
        process_directory(directory,all_files);
        for(int ifile=0;ifile<all_files.size();ifile++)
        {
                if(Contains(all_files[ifile],"_clean_splitup_")) {
                  dout << "      About to delete " << all_files[ifile] << endl;
                  gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
                }
        }
        return true;
  }
  */
  if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) {
    dout << "Going to purge files in local_storage that have been loaded!" << endl;
    for(int i=0;i<(int)SUSYScanSpace::loaded_files.size();i++) {
      gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
      dout << "      Purging : Deleted file " << SUSYScanSpace::loaded_files[i] << endl;
    }
  }
  return true;
}  

bool initialized_t2=false;

int srmcpretries=0;

void load_local_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
    stringstream filetoload;
    string samplename;
    filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
    if(scantype==mSUGRA) {
      //filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
      filetoload.str("");
      filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
      samplename="mSUGRA";
//      filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
    }
    if(scantype==SMS) {
      filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
      samplename="SMS";
    }
    if(scantype==GMSB) {
      filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
      samplename="GMSB";
    }
    
    if(scantype==SMS && Contains(PlottingSetup::ScanSampleDirectory,"T1lh")) {
      filetoload.str("");
      filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/SMS_MassGlu_" << mglu << "__MassLSP_" << mlsp << ".root";
      samplename="SMS";
    }
    
    if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
        dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
        if((scansample.collection).size()>=1) {
           scansample.RemoveLastSample();
        }
        scanfileindex=(scansample.collection).size();
        //New: Loading file when necessary, not before (avoiding high memory usage and startup times)
        
        scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
        dout << "Just added the following file: " << filetoload.str() << endl;
    } else {
        dout << "Last sample is the same as the current one. Recycling it." << endl;
        scanfileindex=(scansample.collection).size()-1;
    }
    dout << " Going to use the following file: " << filetoload.str() << endl;
}

void load_remote_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
    write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
    string samplename;
    PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
    stringstream copyfile;
    gSystem->Exec("mkdir -p local_storage");
    copyfile << "lcg-cp -b -T srmv2 srm://storage01.lcg.cscs.ch:8443/srm/managerv2?SFN=/pnfs/lcg.cscs.ch/cms/trivcat/store/user/buchmann/bussola/" << PlottingSetup::ScanSampleDirectory;
    if(scantype==mSUGRA) {
      copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
      samplename="mSUGRA";
    }
    if(scantype==SMS) {
      copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
      samplename="SMS";
    }
    if(scantype==GMSB) {
      copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
      samplename="GMSB";
    }
    stringstream newfilename;
//    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
    if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
    if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
    
    SUSYScanSpace::SUSYscantype=scantype;
    
    if((scansample.collection).size()==0||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
        dout << "The last sample is NOT the same one as the current one, possibly popping off last one and downloading as well as adding the new one." << endl;
        while((scansample.collection).size()>0) {
           scansample.RemoveLastSample();
        }
        dout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
        scanfileindex=(scansample.collection).size();
        //New: Loading file when necessary, not before (avoiding high memory usage and startup times)
        if(1) {
          if(!doesROOTFileExist(newfilename.str())) {
            dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
            int retcode = gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
            if(retcode==0) SUSYScanSpace::loaded_files.push_back(newfilename.str()); // adding it to the list of loaded files (will be deleted once we finish running
          }
          if(doesROOTFileExist(newfilename.str())) {
                initialized_t2=true;
                dout << "Going to add the following sample: " << newfilename.str() << endl;
                gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
                scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
          } else {
                srmcpretries++;
                if(srmcpretries<5) {
                        dout << "The file could not be loaded correctly - retrying!" << endl;
                        sleep(5);
                        load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
                } else {
                        dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
                        assert(srmcpretries<5);
                        dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
                }
          }
        }
    } else {
        dout << "Last sample is the same as the current one. Recycling it." << endl;
        scanfileindex=(scansample.collection).size()-1;
    }
}  

void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {

// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!

  SUSYScanSpace::SUSYscantype=scantype;
  dout << "Going to load file with Xzone=" << a << " and  Yzone=" << b << " for scantype (mSugra: " << (bool)(scantype==mSUGRA) << " , SMS: " << (bool)(scantype==SMS) << " , GMSB: " << (bool)(scantype==GMSB) << ")" << endl;
  char hostname[1023];
  gethostname(hostname,1023);
  if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) load_local_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
  else load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
}  

/*float get_xs(float &altxs, float mglu, float mlsp, string massgluname, string massLSPname, map <  pair<float, float>, map<string, float>  >  &xsec, string mcjzb, bool requireZ)  {
  altxs=0;
  for(int iproc=1;iproc<=10;iproc++) altxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
  return altxs;
}*/

void establish_SUSY_limits(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, TFile *fsyst, int ibin,float njobs=-1, float jobnumber=-1) {
  bool runninglocally=true;
  if(njobs>-1&&jobnumber>-1) {
    runninglocally=false;
    dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
  } else {
    dout << "Running locally " << endl;
    dout << "This will take a really, really long time - if  you want to see the results within hours instead of weeks try running the worker script on the grid (DistributedModelCalculation/Limits/)" << endl;
  }
 
  string massgluname="MassGlu";
  string massLSPname="MassLSP";
  jzbSel=jzb_cut[ibin];
  geqleq="geq";
  automatized=true;
  mcjzbexpression=mcjzb;

  string prefix="SMS_";
  // up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
  int scantype=SMS;
  TIter nextkey(fsyst->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)nextkey()))
    {
      TObject *obj = key->ReadObj();
      if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
      if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
    }
   
    map <  pair<float, float>, map<string, float>  >  xsec;

  if(scantype==mSUGRA) {
    massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
    massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
    mglustart=m0start;
    xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
    
    mgluend=m0end;
    mglustep=m0step;
    mLSPstart=m12start;
    mLSPend=m12end;
    mLSPstep=m12step;
    prefix="mSUGRA_";
    dout << "mSUGRA scan has been set up." << endl;
  } 
  if(scantype==SMS) {
    dout << "SMS scan has been set up." << endl;
      write_warning(__FUNCTION__,"Don't have the correct XS file yet");
  }
  if(scantype==GMSB) {
    string massgluname="MassGMSBGlu";
    string massLSPname="MassGMSBChi";
    prefix="GMSB_";
    dout << "GMSB scan has been set up." << endl;
      write_warning(__FUNCTION__,"Don't have the correct XS file yet");
  }
 
  set_SUSY_style();
std::cout << __LINE__ << std::endl;
 
  TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
 
  int Npoints=0;
  for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
    for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
  }
  TH2F *mceff     = (TH2F*)fsyst->Get((prefix+"efficiencymap"+any2string(jzb_cut[ibin])).c_str());
//write_warning(__FUNCTION__,"Currently the efficiencymap name was switched to Pablo's convention. This NEEDS to be switched BACK!");TH2F *mceff     = (TH2F*)fsyst->Get(("efficiency_jzbdiff"+any2string(jzb_cut[ibin])).c_str());
  TH2F *fullerr   = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
  TH2F *NEvents   = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
  TH2F *lflip     = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
  TH2F *XSmap     = (TH2F*)fsyst->Get((prefix+"XS"+any2string(jzb_cut[ibin])).c_str());
 
  TH2F *limitmap  = new TH2F((prefix+"limitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
  TH2F *explimitmap  = new TH2F((prefix+"explimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit
  TH2F *exp1plimitmap  = new TH2F((prefix+"exp1plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 1 sigma
  TH2F *exp2plimitmap  = new TH2F((prefix+"exp2plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 2 sigma
  TH2F *exp1mlimitmap  = new TH2F((prefix+"exp1mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 1 sigma
  TH2F *exp2mlimitmap  = new TH2F((prefix+"exp2mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 2 sigma
  
  TH2F *exclmap  = new TH2F((prefix+"exclusionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *flipmap  = new TH2F((prefix+"limitflipmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *asymptoticmap  = new TH2F((prefix+"asymptoticmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);

   
  if(!fullerr || !mceff || !NEvents) {
    write_error(__FUNCTION__,"The supplied systematics file did not contain the correct histograms - please check the file");
    dout << "mc eff address: " << mceff << " , error address: " << fullerr << " , NEvents address: " << NEvents << endl;
    delete limcanvas;
    return;
  }
  
  bool doexpected=true;
 
  int ipoint=-1;
  for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
    for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
    {
      ipoint++;
      if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
      int GlobalBin=flipmap->FindBin(mglu,mlsp); 
      float currmceff=mceff->GetBinContent(mceff->FindBin(mglu,mlsp));
      float currtoterr=(fullerr->GetBinContent(fullerr->FindBin(mglu,mlsp)))*currmceff;
      float nevents=NEvents->GetBinContent(NEvents->FindBin(mglu,mlsp));
      int flipped = (int)(lflip->GetBinContent(lflip->FindBin(mglu,mlsp)));
      flipmap->SetBinContent(GlobalBin,flipped);
      dout << "Looking at point " << ipoint << " / " << Npoints << " with masses " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << endl;
      dout << "Found : MCeff=" << currmceff << " and total error=" << currtoterr << " and Nevents=" << nevents << endl;
      string plotfilename=(string)(TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
      if(currmceff<=0||currtoterr<=0||nevents==0) {
        dout << "   Nothing to work with, skipping this point." << endl;
        continue;
      }
      vector<float> sigmas;
      if(scantype!=mSUGRA) {
        sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
        write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
      }
      else {
        //this is a bit trickier - what we want to do is compute the limit fast, and if it is in [0.5 xs, 2xs], compute it again but "correctly"!
        vector<float> asigmas;
        asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
        float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
/*      if(strength>0.5&&strength<2) {
          sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
          asymptoticmap->SetBinContent(GlobalBin,0);
        } else {*/
          sigmas=asigmas;
          asymptoticmap->SetBinContent(GlobalBin,1);
/*      }*/
        exclmap->SetBinContent(GlobalBin,strength);
      }
        
      if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
      if(sigmas.size()>1) {
        explimitmap->SetBinContent(GlobalBin,sigmas[1]);
        exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
        exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
        exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
        exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
      }

      dout << "An upper limit has been added for this point ( " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << " ) at " << sigmas[0] << endl;
    }
  }
 
  prepare_scan_axis(limitmap,scantype);
  gSystem->Exec("mkdir -p output");
  TFile *outputfile=new TFile(("output/DistributedLimitsFromSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for different JZB cuts and don't want to erase the previous data :-)
  outputfile->cd();
  limitmap->Write();
  XSmap->Write();
  flipmap->Write();
  asymptoticmap->Write();
  if(doexpected) {
    explimitmap->Write();
    exp1plimitmap->Write();
    exp1mlimitmap->Write();
    exp2plimitmap->Write();
    exp2mlimitmap->Write();
  }
  if(scantype==mSUGRA) exclmap->Write();
  outputfile->Close();
  delete limcanvas;
}

void establish_SUSY_limits(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, string fsystfilename, float njobs=-1, float jobnumber=-1) {
  dout << "Starting the SUSY scan from systematics file now with all " << jzb_cut.size() << " bin(s)" << " and using " << fsystfilename << endl;
  TFile *fsyst = new TFile(fsystfilename.c_str());
  if(!fsyst) {
    write_error(__FUNCTION__,"You provided an invalid systematics file. Please change it ... ");
    return;
  }
  for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
    establish_SUSY_limits(mcjzb,datajzb,jzb_cut,requireZ, peakerror, fsyst, ibin,njobs, jobnumber);
  }
}


string SimpleFormat(string input) {
  string output;
  for (unsigned int i=0; i<input.length();i++) {
      if (isalnum(input[i])) //If current character is alpha-numeric
       output+= input[i]; //Add the character to strConverted
 
  }
  return output;
}


vector<float> scan_SUSY_parameter_space_for_one_cut(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, int ibin,float njobs=-1, float jobnumber=-1, bool systematicsonly=false,bool efficiencyonly=false, bool shapeanalysis=false,string specialfilename="") {
  //note: njobs==-2 means that we're running a benchmark point
  dout << "Starting the SUSY scan now with " << jzb_cut[ibin] << " (bin #" << ibin << " )" << endl;
  dout << "   mcjzb: " << mcjzb << endl;
  dout << "   datajzb: " << datajzb << endl;
  dout << "   jzb_cuts: ";
  for(int ib=0;ib<jzb_cut.size()-1;ib++) dout << jzb_cut[ib] << ", ";
  dout << jzb_cut[jzb_cut.size()-1] << endl;
  dout << "   requireZ: " << requireZ << endl;
  dout << "   Peak error: " << peakerror << endl;
  dout << "   Data Peak e: " << peakerrordata << endl;
  dout << "   NJobs : " << njobs << endl;
  dout << "   Job num: " << jobnumber << endl;
  dout << "   Syst only? " << systematicsonly << endl;
  dout << "   Eff only? " << efficiencyonly << endl;
  dout << "   Shape analysis?"  << shapeanalysis << endl;
  dout << "   Special file name? \"" << specialfilename<< "\"" << endl;
  
  vector<float> empty;
  
  if(shapeanalysis&&ibin>0) return empty; // pointless since we're doing a shape analysis :-) 

  bool runninglocally=true;
  if(njobs>-1&&jobnumber>-1) {
    runninglocally=false;
    dout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
  } else {
    dout << "Running locally " << endl;
    if(njobs>-2) dout << "If you are running more than one point and you want to see the results THIS month try running scripts from DistributedModelCalculation" << endl;
  }
 
  jzbSel=jzb_cut[ibin];
  if(shapeanalysis) jzbSel=0;
  geqleq="geq";
  automatized=true;
  mcjzbexpression=mcjzb;
 
  string massgluname="MassGlu";
  string massLSPname="MassLSP";
 
  string prefix="SMS_";
  // up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
  int scantype=SMS;
write_info(__FUNCTION__,PlottingSetup::ScanSampleDirectory);
  if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
  if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;

  if(scantype==mSUGRA) {
    massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
    massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
    mglustart=m0start;
    mgluend=m0end;
    mglustep=m0step;
    mLSPstart=m12start;
    mLSPend=m12end;
    mLSPstep=m12step;
    prefix="mSUGRA_";
    dout << "mSUGRA scan has been set up." << endl;
  }
  if(scantype==SMS) {
    dout << "SMS scan has been set up." << endl;
  }
  if(scantype==GMSB) {
    dout << "GMSB scan has been set up." << endl;
    prefix="GMSB_";
    massgluname="MassGMSBGlu"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
    massLSPname="MassGMSBChi"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
  }
 
  Int_t MyPalette[100];
  Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
  Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
  Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
  Double_t stop[] = {0., .25, .50, .75, 1.0};
  Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
  for (int i=0;i<100;i++) MyPalette[i] = FI+i;
  
  gStyle->SetPalette(100, MyPalette);
 
  TH2F *rawlimitmap  = new TH2F((prefix+"rawlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
  TH2F *rawelimitmap  = new TH2F((prefix+"rawelimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
  
  TH2F *limitmap  = new TH2F((prefix+"limitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
  TH2F *explimitmap  = new TH2F((prefix+"explimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit
  TH2F *exp1plimitmap  = new TH2F((prefix+"exp1plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 1 sigma
  TH2F *exp2plimitmap  = new TH2F((prefix+"exp2plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 2 sigma
  TH2F *exp1mlimitmap  = new TH2F((prefix+"exp1mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 1 sigma
  TH2F *exp2mlimitmap  = new TH2F((prefix+"exp2mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 2 sigma
  
  TH2F *exclmap        = new TH2F((prefix+"exclusionmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *sysjesmap      = new TH2F((prefix+"sysjes"+any2string(jzbSel)).c_str(),   "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *sysjsumap      = new TH2F((prefix+"sysjsu"+any2string(jzbSel)).c_str(),   "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *sysresmap      = new TH2F((prefix+"sysresmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymap  = new TH2F((prefix+"efficiencymap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymapmet  =  new TH2F((prefix+"efficiencymapmet"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymapAcc  =  new TH2F((prefix+"efficiencymapAcc"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymapID   =  new TH2F((prefix+"efficiencymapID"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymapJets =  new TH2F((prefix+"efficiencymapJets"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymapSignal = new TH2F((prefix+"efficiencymapSignal"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *efficiencymapContam = new TH2F((prefix+"efficiencymapContam"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-
mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *noscefficiencymap = new TH2F((prefix+"noscefficiencymap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *Neventsmap     = new TH2F((prefix+"Neventsmap"+any2string(jzbSel)).c_str(),"",   (int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *ipointmap      = new TH2F((prefix+"ipointmap"+any2string(jzbSel)).c_str(),"",   (int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *syspdfmap      = new TH2F((prefix+"syspdfmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *systotmap      = new TH2F((prefix+"systotmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *sysstatmap     = new TH2F((prefix+"sysstatmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *XSmap     = new TH2F((prefix+"XS"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *absXSmap     = new TH2F((prefix+"absXS"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *FilterEff     = new TH2F((prefix+"FilterEfficiency"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  
  TH2F *ScanIdentifier = new TH2F(("ScanIdentifier"+any2string(jzbSel)+SimpleFormat(PlottingSetup::ScanSampleDirectory)).c_str(),("ScanIdentifier"+any2string(jzbSel)+SimpleFormat(PlottingSetup::ScanSampleDirectory)).c_str(),1,0,1,1,0,1);
  
  
  TH2F *asymptoticmap  = new TH2F((prefix+"asymptoticmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  // only used for asymptoticmap

  TH2F *imposedxmap;
  TH2F *realxmap;

  TH2F *centermap;
  TH2F *widthmap;

  TH2F *timemap        = new TH2F((prefix+"timemap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  
  TH2F *flipmap        = new TH2F((prefix+"flipmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  
  if(ibin==0) {
    centermap    = new TH2F((prefix+"centermap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
    widthmap       = new TH2F((prefix+"widthmap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
    
    imposedxmap    = new TH2F((prefix+"imposedxmap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
    realxmap       = new TH2F((prefix+"realxmap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  }

  
  if(shapeanalysis) jzbSel=jzb_cut[ibin];
  vector<int> susy_Zdecay_origin;
  vector<TH2F*> h_susy_Zdecay_origin;
  // we only generate these histograms IF we're dealing with the lowest cut! (we don't need them n times, they're the same for all of them)
  for(int i=1000022;i<1000039;i++) {
    if(ibin>0) continue;
    susy_Zdecay_origin.push_back(i);
    TH2F *susy_dec = new TH2F((prefix+"SUSY_Decay_From_"+any2string(i)+"_to_Z").c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
    h_susy_Zdecay_origin.push_back(susy_dec);
  }
  for(int i=2000006;i<2000016;i++) {
    if(ibin>0) continue;
    susy_Zdecay_origin.push_back(i);
    TH2F *susy_dec = new TH2F((prefix+"SUSY_Decay_From_"+any2string(i)+"_to_Z").c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
    h_susy_Zdecay_origin.push_back(susy_dec);
  }
  
  TH2F *PrimaryZSource    = new TH2F((prefix+"PrimaryZSource").c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  TH2F *SecondaryZSource    = new TH2F((prefix+"SecondaryZSource").c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);

write_warning(__FUNCTION__,"CURRENTLY SWITCHING AUTOMATIZED MODE OFF!");automatized=false;
 
  gStyle->SetPadRightMargin(0.15);
 
  TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
 
  
  bool doexpected=true; 
  
  int Npoints=0;
  for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
    for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
  }
  
  map <  pair<float, float>, map<string, float>  >  xsec;
  if(scantype==mSUGRA) {
    xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
  }

  int ipoint=-1;
  for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
    for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
    {
      ipoint++;
      int GlobalBin=efficiencymap->FindBin(mglu,mlsp);
      if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
      
      int flipped=0;
      float result=-987,resulterr=-987;
      float m0trees = PlottingSetup::mgluend;
      float m12trees = PlottingSetup::mLSPend;
      int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
      int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
      int scanfileindex=PlottingSetup::ScanYzones*a+b;
      if(njobs>-2) load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
      else {
        if((scansample.collection).size()>=1) scansample.RemoveLastSample();
        scanfileindex=(scansample.collection).size();
        scansample.AddSample(specialfilename,"LMpoint",1,1,false,true,scanfileindex,kRed);
      }
      
      clock_t start,finish;
      start = clock(); // starting the clock to measure how long the computation takes!
      stringstream addcut;
      addcut << "(TMath::Abs("<<massgluname<<"-"<<mglu<<")<5)&&(TMath::Abs("<<massLSPname<<"-"<<mlsp<<")<5)";
      if(njobs==-2) {
        dout << "Processing a signal sample, not a scan : setting addcut to \"eventNum>0||eventNum<1\" (empty cut)" << endl;
        addcut.str("(eventNum>0||eventNum<1)");
      }

      (scansample.collection)[scanfileindex].events->Draw("eventNum",addcut.str().c_str(),"goff");
      float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
      vector<vector<float> > systematics;
      if(nevents<10) {
        dout << "This point ("<<ipoint<<" / " << Npoints << ") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") does not contain enough events and will be \033[1;31m skipped \033[0m."<< endl;
        continue;
      } else {
        dout << "OK! This point ("<<ipoint<<") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") contains " << nevents << " and will therefore \033[1;32m not be skipped \033[0m. " << endl;
        if(PlottingSetup::RestrictToMassPeak==true) dout << "Analysis type: on-peak" << endl;
        if(PlottingSetup::RestrictToMassPeak==false) dout << "Analysis type: offpeak" << endl;
        if(shapeanalysis) dout << "Shape analysis? Yes." << endl;
        else dout << "Shape analysis? No." << endl;
      }
      float hitrecord=0;
      int hitrecordholder=0;
      float hitsecondplace=0;
      int hitsecondplaceholder=0;
      for(int imo=0;imo<(int)susy_Zdecay_origin.size()&&ibin==0&&njobs>-2;imo++) {
        if(ibin>0) continue;
        int MotherParticlePDGid = susy_Zdecay_origin[imo];
        (scansample.collection)[scanfileindex].events->Draw("eventNum",("("+addcut.str()+")&&(abs(SourceOfZ-"+any2string(MotherParticlePDGid)+")<0.5)").c_str(),"goff");
        float hits = (scansample.collection)[scanfileindex].events->GetSelectedRows();
        hits/=nevents;
        h_susy_Zdecay_origin[imo]->SetBinContent(GlobalBin,hits);
        if(hits>hitrecord) {
          hitsecondplace=hitrecord;
          hitsecondplaceholder=hitrecordholder;
          hitrecord=hits;
          hitrecordholder=imo;
        }
        if(hits<hitrecord&&hits>hitsecondplace) {
          hitsecondplace=hits;
          hitsecondplaceholder=imo;
        }
      }
      if(ibin==0) {
        PrimaryZSource->SetBinContent(GlobalBin,hitrecordholder);
        SecondaryZSource->SetBinContent(GlobalBin,hitsecondplaceholder);
      }
      
      if(ibin==0&&njobs>-2) {
        TH1F *realxhisto = new TH1F("realxhisto","realxhisto",1000,0,1);
        (scansample.collection)[scanfileindex].events->Draw("realx>>realxhisto",addcut.str().c_str(),"goff");
        realxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
        realxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
        (scansample.collection)[scanfileindex].events->Draw("imposedx>>realxhisto",addcut.str().c_str(),"goff");
        imposedxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
        imposedxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
        stringstream specialexpression;
        specialexpression << mcjzb << ">>jzbhisto";
        TH1F *jzbhisto = new TH1F("jzbhisto","jzbhisto",100,-500,500);
        (scansample.collection)[scanfileindex].events->Draw(specialexpression.str().c_str(),addcut.str().c_str(),"goff");
        centermap->SetBinContent(GlobalBin,jzbhisto->GetMean());
        widthmap->SetBinContent(GlobalBin,jzbhisto->GetRMS());
        delete realxhisto;
        delete jzbhisto;
      }


      if(scantype==mSUGRA) {
        float absxs=0;
        for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
        TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
        TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
        float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
        FilterEff->SetBinContent(GlobalBin,filtereff);
        FilterEffFile->Close();
        absXSmap->SetBinContent(GlobalBin,absxs);
        absxs*=filtereff;
        XSmap->SetBinContent(GlobalBin,absxs);
      }
        
      if(nevents!=0) {
          Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
          if(result<0&&allowflipping) {
                write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
                flipped=1;
                MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
          }
          efficiencymap->SetBinContent(GlobalBin,result);
          efficiencymap->SetBinError(GlobalBin,resulterr);
          flipmap->SetBinContent(GlobalBin,flipped);
          noscefficiencymap->SetBinContent(GlobalBin,effwosigcont.getValue());
          noscefficiencymap->SetBinError(GlobalBin,effwosigcont.getError());
          efficiencymapContam->SetBinContent(GlobalBin,effwosigcont.getValue()-result);
          efficiencymapContam->SetBinError(GlobalBin,TMath::Sqrt(effwosigcont.getError()+resulterr));
          //Partial efficiencies
          float resultAcc, resultID, resultJets, resultSignal, resultmet;
          float resulterrAcc, resulterrID, resulterrJets, resulterrSignal, resulterrmet;
          MCPartialefficiency((scansample.collection)[scanfileindex].events,resultAcc,resulterrAcc,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 1);
          MCPartialefficiency((scansample.collection)[scanfileindex].events,resultID,resulterrID,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 2);
          MCPartialefficiency((scansample.collection)[scanfileindex].events,resultJets,resulterrJets,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 3);
          MCPartialefficiency((scansample.collection)[scanfileindex].events,resultSignal,resulterrSignal,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 4);
          MCPartialefficiency((scansample.collection)[scanfileindex].events,resultmet,resulterrmet,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 5);
          efficiencymapAcc->SetBinContent(GlobalBin,resultAcc);
          efficiencymapAcc->SetBinError(GlobalBin,resulterrAcc);
          efficiencymapmet->SetBinContent(GlobalBin,resultmet);
          efficiencymapmet->SetBinError(GlobalBin,resulterrmet);
          efficiencymapID->SetBinContent(GlobalBin,resultID);
          efficiencymapID->SetBinError(GlobalBin,resulterrID);
          efficiencymapJets->SetBinContent(GlobalBin,resultJets);
          efficiencymapJets->SetBinError(GlobalBin,resulterrJets);
          efficiencymapSignal->SetBinContent(GlobalBin,resultSignal); 
          efficiencymapSignal->SetBinError(GlobalBin,resulterrSignal);
          finish = clock();
      }
      Neventsmap->SetBinContent(GlobalBin,nevents);
      ipointmap->SetBinContent(GlobalBin,ipoint);
      if(efficiencyonly) continue;
      
      if(!shapeanalysis) {
        do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
//      float JZBcutat = systematics[0][0];
        float mceff    = systematics[0][1];
        float mcefferr = systematics[0][2];//MC stat error
        float toterr   = systematics[0][4];
        float sys_jes  = systematics[0][5]; // Jet Energy Scale
        float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
        float sys_res  = systematics[0][7]; // resolution
//      float mcwoscef = systematics[0][8]; // efficiency without signal contamination
//      float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
        float sys_pdf   = 0;
        if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
        dout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
        if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
          dout << "Limits can't be calculated for this configuration (" << massgluname <<"="<<mglu<<" , " << massLSPname << "="<<mlsp << ") as either the efficiency or its error are not positive numbers! (mceff="<<mceff<<"and toterr="<<toterr<<")"<< endl;
          continue;
        } else {
          // Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
          if(!systematicsonly&&!efficiencyonly) {
            dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
            vector<float> sigmas;
            string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
            sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
//          do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
            if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
              limitmap->SetBinContent(GlobalBin,sigmas[0]);
              if(sigmas.size()>1) {
                explimitmap->SetBinContent(GlobalBin,sigmas[1]);
                exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
                exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
                exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
                exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
              }
              sysjesmap->SetBinContent(GlobalBin,sys_jes);
              sysjsumap->SetBinContent(GlobalBin,sys_jsu);
              sysresmap->SetBinContent(GlobalBin,sys_res);
              syspdfmap->SetBinContent(GlobalBin,sys_pdf);
              systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
              sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
              dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
            }//end of if sigma is positive
          }//end of not systematics only condition
          if(systematicsonly) {
            sysjesmap->SetBinContent(GlobalBin,sys_jes);
            sysjsumap->SetBinContent(GlobalBin,sys_jsu);
            sysresmap->SetBinContent(GlobalBin,sys_res);
            syspdfmap->SetBinContent(GlobalBin,sys_pdf);
            systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
            sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
          }
        }//efficiency is valid
      } 
      if(shapeanalysis) {
        //shape analysis
        dout << "This is the shape analysis. Have fun with it!" << endl;
        float quickmceff,quickmcefferr;
        MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
        if(quickmceff==0) {
          write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
          continue;
        }
        stringstream PointName;
        PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
        SUSYScanSpace::SavedMLSP=mlsp;
        SUSYScanSpace::SavedMGlu=mglu;
        SUSYScanSpace::SavedMLSPname=massLSPname;
        SUSYScanSpace::SavedMGluname=massgluname;
        float referenceXS=0;
        if(scantype==mSUGRA) referenceXS=XSmap->GetBinContent(GlobalBin);
        else referenceXS=getSMSxs(mlsp,mglu);
        float low_fullCLs=referenceXS*0.5;
        float high_fullCLs=referenceXS*2.5;
        if(scantype!=mSUGRA) {
          low_fullCLs=-1;
          high_fullCLs=10e10;
        }
        ShapeDroplet droplet = LimitsFromShapes(low_fullCLs,high_fullCLs, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
        if(!isValidShapeResult(droplet)) {
          write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - ignoring this point!!!");
//        scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, (float)Npoints, (float)ipoint, systematicsonly, efficiencyonly,false); // don't uncomment this line; it would lead to multiple limitmap histos ...
          if(njobs==-2) {
            delete limcanvas;
            return empty;
          }
          continue;

        }

        float excludedXSobs = droplet.observed/droplet.SignalIntegral;
        absXSmap->SetBinContent(GlobalBin,referenceXS);
        write_info(__FUNCTION__,"Established limit at "+any2string(excludedXSobs)+" , ref XS is "+any2string(referenceXS));
        if(njobs==-2) {
          //this is the non-scan case
          vector<float> sigmas;
          sigmas.push_back(droplet.observed/droplet.SignalIntegral);
          sigmas.push_back(droplet.expected/droplet.SignalIntegral);
          sigmas.push_back(droplet.expectedPlus1Sigma/droplet.SignalIntegral);
          sigmas.push_back(droplet.expectedMinus1Sigma/droplet.SignalIntegral);
          sigmas.push_back(droplet.expectedPlus2Sigma/droplet.SignalIntegral);
          sigmas.push_back(droplet.expectedMinus2Sigma/droplet.SignalIntegral);
          delete limcanvas;
          return sigmas;
        }
          
        rawlimitmap->SetBinContent(GlobalBin,droplet.observed);
        rawelimitmap->SetBinContent(GlobalBin,droplet.expected);
        limitmap->SetBinContent(GlobalBin,droplet.observed/droplet.SignalIntegral);
        explimitmap->SetBinContent(GlobalBin,droplet.expected/droplet.SignalIntegral);
        exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma/droplet.SignalIntegral);
        exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma/droplet.SignalIntegral);
        exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma/droplet.SignalIntegral);
        exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma/droplet.SignalIntegral);
        dout <<"Saved the following shape limit: " << droplet.observed/droplet.SignalIntegral << " (observed) , " << droplet.expected/droplet.SignalIntegral << " (expected) " << endl;
        sysjesmap->SetBinContent(GlobalBin,droplet.JES);
        sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
        syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
        systotmap->SetBinContent(GlobalBin,droplet.toterr);
        sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
      }
  finish = clock();
  timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
    }
  }
  prepare_scan_axis(limitmap,scantype);
  prepare_scan_axis(sysjesmap,scantype);
  prepare_scan_axis(sysjsumap,scantype);
  prepare_scan_axis(sysresmap,scantype);
  prepare_scan_axis(syspdfmap,scantype);
  prepare_scan_axis(systotmap,scantype);
  prepare_scan_axis(sysstatmap,scantype);
  prepare_scan_axis(timemap,scantype);
 
  if(!systematicsonly&&!efficiencyonly) {
    TObject* old=gDirectory->GetList()->FindObject("limcanvas");
    if(old) gDirectory->GetList()->Remove(old);
    TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
    limitmap->Draw("COLZ");
    string titletobewritten="Limits in LSP-Glu plane";
    if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
    if(scantype==GMSB) titletobewritten="Limits in Chi-Glu plane";
    TText *title = write_title(titletobewritten);
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/Limits_JZB_geq"+any2string(jzb_cut[ibin]));
    } else {
      TFile *outputfile;
      if(systematicsonly) outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for different JZB cuts and don't want to erase the previous data :-)
    else outputfile=new TFile(("output/DistributedLimits_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for
      limitmap->Write();
      rawlimitmap->Write();
      rawelimitmap->Write();
      if(doexpected) {
        explimitmap->Write();
        exp1plimitmap->Write();
        exp1mlimitmap->Write();
        exp2plimitmap->Write();
        exp2mlimitmap->Write();
      }
      sysjesmap->Write();
      sysjsumap->Write();
      sysresmap->Write();
      efficiencymap->Write();
      if(ibin==0) centermap->Write();
      if(ibin==0) widthmap->Write();
      flipmap->Write();
      efficiencymapmet->Write();
      efficiencymapAcc->Write();
      efficiencymapID->Write();
      efficiencymapJets->Write();
      efficiencymapSignal->Write();
      noscefficiencymap->Write();
      efficiencymapContam->Write();
      syspdfmap->Write();
      systotmap->Write();
      sysstatmap->Write();
      XSmap->Write();
      absXSmap->Write();
      FilterEff->Write();
      if(shapeanalysis) asymptoticmap->Write();
      if(ibin==0) imposedxmap->Write();
      if(ibin==0) realxmap->Write();
      Neventsmap->Write();
      ScanIdentifier->Write();
      ipointmap->Write();
      timemap->Write();
      for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
        h_susy_Zdecay_origin[i]->Write();
      }
      outputfile->Close();
    }
  }
  if(systematicsonly) { // systematics only :
    limcanvas->cd();
    sysjesmap->Draw("COLZ");
    string titletobewritten="Jet Energy scale in LSP-Glu plane";
    if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
    if(scantype==GMSB) titletobewritten="Limits in Chi-Glu plane";
    TText *title = write_title(titletobewritten);
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/JES_geq"+any2string(jzb_cut[ibin]));
    }
    sysjsumap->Draw("COLZ");
    titletobewritten="JZB Scale Uncertainty in LSP-Glu plane";
    if(scantype==mSUGRA) titletobewritten="JZB Scale Uncertainty in m_{1/2}-m_{0} plane";
    if(scantype==GMSB) titletobewritten="JZB Scale Uncertainty in Chi-Glu plane";
    TText *title2 = write_title(titletobewritten);
    title2->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/JSU_geq"+any2string(jzb_cut[ibin]));
    }
    sysresmap->Draw("COLZ");
    titletobewritten="Resolution in LSP-Glu plane";
    if(scantype==mSUGRA) titletobewritten="Resolution in m_{1/2}-m_{0} plane";
    if(scantype==GMSB) titletobewritten="Resolution in Chi-Glu plane";
    TText *title3 = write_title(titletobewritten);
    title3->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/Resolution_geq"+any2string(jzb_cut[ibin]));
    }
   
    if(!runninglocally) {
      TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
      sysjesmap->Write();
      sysjsumap->Write();
      sysresmap->Write();
      flipmap->Write();
      if(ibin==0) centermap->Write();
      if(ibin==0) widthmap->Write();
      efficiencymap->Write();
      efficiencymapmet->Write();
      efficiencymapAcc->Write();
      efficiencymapID->Write();
      efficiencymapJets->Write();
      efficiencymapSignal->Write();
      noscefficiencymap->Write();
      efficiencymapContam->Write();
      Neventsmap->Write();
      ScanIdentifier->Write();
      ipointmap->Write();
      syspdfmap->Write();
      systotmap->Write();
      sysstatmap->Write();
      absXSmap->Write();
      XSmap->Write();
      FilterEff->Write();
      if(shapeanalysis) asymptoticmap->Write();
      if(ibin==0) imposedxmap->Write();
      if(ibin==0) realxmap->Write();
      timemap->Write();
      for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
        h_susy_Zdecay_origin[i]->Write();
      }
      outputfile->Close();
    }
  }//end of systematics only
  if(efficiencyonly) {
    limcanvas->cd();
    string titletobewritten="Efficiencies in LSP-Glu plane";
    if(scantype==mSUGRA) titletobewritten="Efficiencies in m_{1/2}-m_{0} plane";
    if(scantype==GMSB) titletobewritten="Efficiencies in Chi-Glu plane";
    TText *title = write_title(titletobewritten);
    efficiencymap->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/Efficiency_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    efficiencymapmet->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/EfficiencyMet_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    efficiencymapAcc->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/EfficiencyAcc_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    efficiencymapID->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/EfficiencyID_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    efficiencymapJets->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/EfficiencyJets_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    efficiencymapSignal->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/EfficiencySignal_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    noscefficiencymap->Draw("COLZ");
    title->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/NoEfficiency_geq"+any2string(jzb_cut[ibin]));
    }
    limcanvas->cd();
    sysjesmap->Draw("COLZ");
    titletobewritten="N(events) in LSP-Glu plane";
    if(scantype==mSUGRA) titletobewritten="N(events) in m_{1/2}-m_{0} plane";
    if(scantype==GMSB) titletobewritten="Resolution in Chi-Glu plane";
    TText *title2 = write_title(titletobewritten);
    title2->Draw();
    if(runninglocally) {
      CompleteSave(limcanvas,"SUSYScan/Nevents_geq"+any2string(jzb_cut[ibin]));
    }
    if(!runninglocally) {
      TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
      ipointmap->Write();
      Neventsmap->Write();
      ScanIdentifier->Write();
      noscefficiencymap->Write();
      efficiencymapContam->Write();
      efficiencymap->Write();
      if(ibin==0) centermap->Write();
      if(ibin==0) widthmap->Write();
      flipmap->Write();
      efficiencymapmet->Write();
      efficiencymapAcc->Write();
      efficiencymapID->Write();
      efficiencymapJets->Write();
      efficiencymapSignal->Write();
      timemap->Write();
      outputfile->Close();
    }
  }//end of efficiencies only

 for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
    delete h_susy_Zdecay_origin[i];
  }
  
  delete PrimaryZSource;
  delete SecondaryZSource;

 
  delete limitmap;
  delete exclmap;
  delete sysjesmap;
  delete sysjsumap;
  delete sysresmap;
  delete noscefficiencymap;
  delete efficiencymapContam;
  delete efficiencymap;
  delete efficiencymapmet;
  delete efficiencymapAcc;
  delete efficiencymapID;
  delete efficiencymapJets;
  delete efficiencymapSignal;
  delete Neventsmap;
  delete ScanIdentifier;
  delete ipointmap;
  delete syspdfmap;
  delete systotmap;
  delete sysstatmap;
  delete XSmap;
  delete absXSmap;
  TObject* old=gDirectory->GetList()->FindObject("limcanvas");
  if(old) gDirectory->GetList()->Remove(old);
  return empty;
}

void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, float njobs=-1, float jobnumber=-1, bool systonly=false, bool effonly=false,bool shapeanalysis=false) {
  dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
  dout << "   mcjzb: " << mcjzb << endl;
  dout << "   datajzb: " << datajzb << endl;
  dout << "   jzb_cut: ";
  for(int ib=0;ib<jzb_cut.size()-1;ib++) dout << jzb_cut[ib] << ", ";
  dout << jzb_cut[jzb_cut.size()-1] << endl;
  dout << "   requireZ: " << requireZ << endl;
  dout << "   Peak error: " << peakerror << endl;
  dout << "   Data Peak e: " << peakerrordata << endl;
  dout << "   NJobs : " << njobs << endl;
  dout << "   Job num: " << jobnumber << endl;
  dout << "   Syst only?" << systonly << endl;
  dout << "   Eff only? " << effonly << endl;
  dout << "   Shape analysis?" << shapeanalysis << endl;
  
  for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
    scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
  }
  delete_any_cached_scans();// tidy up ...
}


void compute_upper_limits_from_shapes(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata) {
  dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
  vector<vector<float> > limits;
  vector<string> samplename;
  for(int isample=0;isample<signalsamples.collection.size();isample++) {
    limits.push_back(scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, 0, -2, -2,false,false,true,signalsamples.collection[isample].filename));// the "-2" stands for "load from customized path"
  }
  for(int isample=0;isample<signalsamples.collection.size();isample++) {
    dout << signalsamples.collection[isample].samplename << " :" << endl;
    dout << "       observed :   " << limits[isample][0] << endl;
    dout << "       expected :   " << limits[isample][1] << endl;
    dout << "       exp. band:   [ " << limits[isample][5] << " , " << limits[isample][3] << " , \033[1;34m " << limits[isample][1] << "\033[0m , " << limits[isample][2] << " , " << limits[isample][4] << "] " << endl;
  }
}
Revision:	1.19
Committed:	Wed May 30 11:45:32 2012 UTC (12 years, 11 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.18:	+38 -1 lines
Log Message:	Being more verbose about the SUSY scan parameters with which the algorithm was called