OOTMultiplicityPlotMacros.cc

Go to the documentation of this file.

#include "OOTMultiplicityPlotMacros.h"
#include "DPGAnalysis/SiStripTools/interface/CommonAnalyzer.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TF1.h"
#include "TProfile.h"
#include <iostream>
#include <algorithm>
#include <cmath>

OOTSummary* ComputeOOTFractionvsFill(TFile* ff, const char* itmodule, const char* ootmodule, const char* etmodule, const char* hname, OOTSummary* ootsumm) {

  if(ootsumm==nullptr) {
    ootsumm = new OOTSummary; 
  }    
  
  if(ff) {
    CommonAnalyzer ca(ff,"",itmodule);
    std::vector<unsigned int> runs = ca.getFillList();
    std::sort(runs.begin(),runs.end());
    for(unsigned int i=0;i<runs.size();++i) {

      char runlabel[100];
      sprintf(runlabel,"%d",runs[i]);

      OOTResult* res = ComputeOOTFraction(ff,itmodule,ootmodule,etmodule,runs[i],hname,true);

      if(res->ngoodbx!=res->hratio->GetEntries()) std::cout << "Inconsistency in number of good bx" << std::endl;
      ootsumm->hngoodbx->Fill(runlabel,res->ngoodbx);
      int ibin = ootsumm->hootfracsum->Fill(runlabel,res->ootfracsum);
      ootsumm->hootfracsum->SetBinError(ibin,res->ootfracsumerr);
      int bin = ootsumm->hootfrac->Fill(runlabel,res->ootfrac);
      ootsumm->hootfrac->SetBinError(bin,res->ootfracerr);
      delete res;

    }
  }
  else {std::cout << "File is not ok" << std::endl; }

  return ootsumm;
}
OOTSummary* ComputeOOTFractionvsRun(TFile* ff, const char* itmodule, const char* ootmodule, const char* etmodule, const char* hname, OOTSummary* ootsumm) {

  if(ootsumm==nullptr) {
    ootsumm =new OOTSummary;
  }

  if(ff) {
    CommonAnalyzer ca(ff,"",itmodule);
    std::vector<unsigned int> runs = ca.getRunList();
    std::sort(runs.begin(),runs.end());
    for(unsigned int i=0;i<runs.size();++i) {

      char runlabel[100];
      sprintf(runlabel,"%d",runs[i]);

      OOTResult* res = ComputeOOTFraction(ff,itmodule,ootmodule,etmodule,runs[i],hname);

      if(res->ngoodbx!=res->hratio->GetEntries()) std::cout << "Inconsistency in number of good bx" << std::endl;
      ootsumm->hngoodbx->Fill(runlabel,res->ngoodbx);
      int ibin = ootsumm->hootfracsum->Fill(runlabel,res->ootfracsum);
      ootsumm->hootfracsum->SetBinError(ibin,res->ootfracsumerr);
      int bin = ootsumm->hootfrac->Fill(runlabel,res->ootfrac);
      ootsumm->hootfrac->SetBinError(bin,res->ootfracerr);
      delete res;

    }
  }
  else {std::cout << "File is not ok" << std::endl; }

  return ootsumm;
}

OOTResult* ComputeOOTFraction(TFile* ff, const char* itmodule, const char* ootmodule, const char* etmodule, const int run, const char* hname, const bool& perFill) {

  if(perFill)  {std::cout << "Processing fill " << run << std::endl;} else {std::cout << "Processing run " << run << std::endl;}

  char itpath[100];
  char ootpath[100];
  char etpath[100];
  if(perFill) {
    sprintf(itpath,"%s/fill_%d",itmodule,run);
    sprintf(ootpath,"%s/fill_%d",ootmodule,run);
    sprintf(etpath,"%s/fill_%d",etmodule,run);
  }
  else {
    sprintf(itpath,"%s/run_%d",itmodule,run);
    sprintf(ootpath,"%s/run_%d",ootmodule,run);
    sprintf(etpath,"%s/run_%d",etmodule,run);
  }

  OOTResult* res = new OOTResult;

  std::vector<int> filledbx = FillingSchemeFromProfile(ff,itpath,hname);
  res->nfilledbx = filledbx.size();

  if(!perFill) {
    std::vector<int> filledbxtest = FillingScheme(ff,etpath);
    if(filledbx.size() != filledbxtest.size()) std::cout << "Inconsistency in number of filled BX " 
                             << run << " " << filledbx.size() << " " << filledbxtest.size() << std::endl;
  }

  TProfile* itmult = nullptr;
  TProfile* ootmult = nullptr;
  res->hratio = new TH1F("ratio","ratio",200,0.,2.);
  
  float rclzb=0;
  float rclrandom=0;
  float errclzb=0;
  float errclrandom=0;
  float nzb=0;
  float nrandom=0;
  
  if(ff) {
    if(ff->cd(itpath)) {
      itmult = (TProfile*)gDirectory->Get(hname);
    }
    else { std::cout << "In time path is not ok" << std::endl;}
    if(ff->cd(ootpath)) {
      ootmult = (TProfile*)gDirectory->Get(hname);
    }
    else { std::cout << "out of time path is not ok" << std::endl;}
    if(itmult && ootmult) {
      ootmult->SetLineColor(kRed);      ootmult->SetMarkerColor(kRed);
      //      ootmult->Draw();
      //      itmult->Draw("same");
      for(std::vector<int>::const_iterator fbx=filledbx.begin();fbx!=filledbx.end();++fbx) {
    nzb += itmult->GetBinEntries(*fbx);
    nrandom += ootmult->GetBinEntries(*fbx+1);
      }
      for(std::vector<int>::const_iterator fbx=filledbx.begin();fbx!=filledbx.end();++fbx) {
    if(nzb > 0 && nrandom > 0) {
      rclzb += (itmult->GetBinContent(*fbx)*itmult->GetBinEntries(*fbx))/nzb;
      errclzb += (itmult->GetBinError(*fbx)*itmult->GetBinEntries(*fbx))*(itmult->GetBinError(*fbx)*itmult->GetBinEntries(*fbx))/(nzb*nzb);
      rclrandom += (ootmult->GetBinContent(*fbx+1)*ootmult->GetBinEntries(*fbx+1))/nrandom;
      errclrandom += (ootmult->GetBinError(*fbx+1)*ootmult->GetBinEntries(*fbx+1))*(ootmult->GetBinError(*fbx+1)*ootmult->GetBinEntries(*fbx+1))
        /(nrandom*nrandom);
    }
    if(itmult->GetBinContent(*fbx)==0) { std::cout << "No cluster in filled BX! " << *fbx << std::endl;}
    else if(ootmult->GetBinEntries(*fbx+1)==0) {/* std::cout << "No entry in OOT BX " << *fbx+1 << std::endl; */}
    else {
      float rat= ootmult->GetBinContent(*fbx+1)/itmult->GetBinContent(*fbx);
      res->hratio->Fill(rat);
    }     
      }
    }
    else {std::cout<<"histograms not found"<<std::endl;}
  }
  else { std::cout << "Input file pointer is not ok" <<std::endl;}

  res->ngoodbx = res->hratio->GetEntries();

  if(nzb>0 && nrandom>0 &&rclzb>0 && rclrandom > 0) {
    res->ootfracsum = rclrandom/rclzb;
    res->ootfracsumerr = rclrandom/rclzb*sqrt(errclzb/(rclzb*rclzb)+errclrandom/(rclrandom*rclrandom));
  }
  if(res->ngoodbx) {
    res->hratio->Fit("gaus","Q0","",.01,1.99);
    if(res->hratio->GetFunction("gaus")) {
      res->ootfrac = res->hratio->GetFunction("gaus")->GetParameter(1);
      res->ootfracerr = res->hratio->GetFunction("gaus")->GetParError(1);
    }
    else {std::cout << "Missing fitting function" << std::endl;}
  }
  else {std::cout << "No filled BX or strange filling scheme" <<std::endl;} 

  return res;
}

std::vector<int> FillingScheme(TFile* ff, const char* path, const float thr) {

  TH1F* bx=nullptr;
  std::vector<int> filledbx;
  if(ff) {
    if(ff->cd(path)) {
      bx = (TH1F*)gDirectory->Get("bx");
      if(bx) {
    //  bx->Draw();
    std::cout << "Number of entries " << bx->GetEntries() << " threshold " << thr/3564.*bx->GetEntries() << std::endl;
    for(int i=1;i<bx->GetNbinsX()+1;++i) {
      if(bx->GetBinContent(i)>thr/3564.*bx->GetEntries()) {
        if(!filledbx.empty() && i == filledbx[filledbx.size()-1]+1) {
          std::cout << "This is not a 50ns run ! " << std::endl;
          filledbx.clear();
          return filledbx;
        }
        filledbx.push_back(i);
      }
    }
      }
      else { std::cout << "Histogram not found" << std::endl;}
    }
    else { std::cout << "module path is not ok" << std::endl;}
  }
  else { std::cout << "Input file pointer is not ok" <<std::endl;}
  
  //  std::cout << filledbx.size() << " filled bunch crossings" << std::endl;
  //  for(std::vector<int>::const_iterator fbx=filledbx.begin();fbx!=filledbx.end();++fbx) { std::cout << *fbx << std::endl;}
  return filledbx;
}
std::vector<int> FillingSchemeFromProfile(TFile* ff, const char* path, const char* hname, const float thr) {

  TProfile* bx=nullptr;
  std::vector<int> filledbx;
  if(ff) {
    if(ff->cd(path)) {
      bx = (TProfile*)gDirectory->Get(hname);
      if(bx) {
    //  bx->Draw();
    std::cout << "Number of entries " << bx->GetEntries() << " threshold " << thr/3564.*bx->GetEntries() << std::endl;
    for(int i=1;i<bx->GetNbinsX()+1;++i) {
      if(bx->GetBinEntries(i)>thr/3564.*bx->GetEntries()) {
        if(!filledbx.empty() && i == filledbx[filledbx.size()-1]+1) {
          std::cout << "This is not a 50ns run ! " << std::endl;
          filledbx.clear();
          return filledbx;
        }
        filledbx.push_back(i);
      }
    }
      }
      else { std::cout << "Histogram not found" << std::endl;}
    }
    else { std::cout << "module path is not ok" << std::endl;}
  }
  else { std::cout << "Input file pointer is not ok" <<std::endl;}
  
  //  std::cout << filledbx.size() << " filled bunch crossings" << std::endl;
  //  for(std::vector<int>::const_iterator fbx=filledbx.begin();fbx!=filledbx.end();++fbx) { std::cout << *fbx << std::endl;}
  return filledbx;
}