// // deplVandLumiExtract.C // macro to extract the depletion voltage measurements after production // and the luminosity information from the raw data // // Created by Christian Elsasser on 31.05.13. // University of Zurich, elsasser@cern.ch // Copyright only for commercial use // // General include #include <iostream> #include <fstream> #include <string.h> #include <vector> #include <map> #include <locale> // ROOT include #include "../../include/incBasicROOT.h" #include "../../include/incDrawROOT.h" #include "../../include/incHistROOT.h" #include "../../include/incIOROOT.h" #include "STDatabase.h" #include "../../include/basicROOTTools.h" #include "../../include/basicRooFitTools.h" #include "../../include/lhcbstyle.C" int deplVandLumiExtract(int, char**); int main(int argc, char* argv[]){ TStyle* style = lhcbStyle(); int argc_copy = argc; TApplication* theApp = new TApplication("Analysis", &argc, argv,0,-1); argc = theApp->Argc(); argv = theApp->Argv(); bool runBatch = false; for (int i = 1; i<argc; i++) { if (strcmp(argv[i],"-b")==0) { runBatch = true; gROOT->SetBatch(); } } style->SetPadLeftMargin(0.12); style->SetPadRightMargin(0.22); style->SetPadBottomMargin(0.15); style->SetPadTopMargin(0.05); int exit = deplVandLumiExtract(argc,argv); Printf("End"); if (!runBatch) { theApp->Run(!runBatch); } return exit; } // Executable method int deplVandLumiExtract(int argc, char* argv[]){ bool useRadius = false; TString det("TT"); for (int i=1; i<argc; i++) { if (strcmp(argv[i],"-d")==0 && i+1<argc) { det = argv[++i]; } } TString lay; if (strcmp(det.Data(),"TT")==0) { lay = "TTaU"; }else{ lay = "T3X2"; } Info("deplVandLumiExtrac","Extract V_depl for %s / %s",lay.Data(),det.Data()); const char* homeVar = std::getenv ("CCEHOME"); const char* diskVar = std::getenv ("DISK"); TString dn_lumi(Form("%s/data/lumi",homeVar)); TString fn_lumi(Form("%s/lumi.root", dn_lumi.Data())); TString dn_out(Form("%s/data/db",homeVar)); TString fn_vdepl(Form("%s/vdepl.root", dn_out.Data())); TString fn_flux(Form("%s/lumi.root", dn_out.Data())); TString fn_db(Form("%s/%s_db.txt", dn_out.Data(), det.Data())); Printf("%s",fn_db.Data()); std::ifstream f_db(fn_db.Data()); std::vector<double> v_sector; std::vector<double> v_vdepl; std::vector<double> v_lumi7; std::vector<double> v_lumi8; std::vector<double> v_lumi13; std::vector<double> v_lumi14; std::vector<double> v_time; std::vector<double> v_fill; while (f_db) { std::string line; if (!getline(f_db,line)) { break; } std::istringstream s_line(line); std::vector<std::string> v_line; while (s_line) { std::string bit; if (!getline(s_line,bit,'\t')) { break; } v_line.push_back(bit); } if (v_line.size()>6 && strcmp(v_line[4].c_str(),lay.Data())==0) { Info("deplVandLumiExtract","Using sector %s",v_line[0].c_str()); v_sector.push_back(atof(v_line[0].c_str())); v_vdepl.push_back(atof(v_line[6].c_str())); }else{ Warning("deplVandLumiExtract","Database entry of %s has too few entries",v_line[0].c_str()); } } for (int i=0; i<v_sector.size(); i++) { Info("deplVandLumiExtract","Sector %d has %6.2f V", (int)v_sector[i],v_vdepl[i]); } TFile* f_vdepl = new TFile(fn_vdepl.Data(),"UPDATE"); TVectorD vr_sector = getROOTVector<double>(v_sector); TVectorD vr_vdepl = getROOTVector<double>(v_vdepl); f_vdepl->Delete(Form("v_sector_%s;*", lay.Data())); f_vdepl->Delete(Form("v_vdepl_%s;*",lay.Data())); vr_sector.Write(Form("v_sector_%s", lay.Data())); vr_vdepl.Write(Form("v_vdepl_%s",lay.Data())); f_vdepl->Close(); Info("deplVandLumiExtract","Results written to %s",fn_vdepl.Data()); TFile* f_lumi = TFile::Open(fn_lumi.Data()); TVectorF vr_fillTime = *(TVectorF*)f_lumi->Get("v_fillTime_tot"); TVectorF vr_stopTime = *(TVectorF*)f_lumi->Get("v_stopTime_tot"); TVectorF vr_beamEner = *(TVectorF*)f_lumi->Get("v_beamEner_tot"); TVectorF vr_peakLumi = *(TVectorF*)f_lumi->Get("v_peakLumi_tot"); if (vr_fillTime.GetNoElements()!=vr_beamEner.GetNoElements() || vr_fillTime.GetNoElements()!=vr_peakLumi.GetNoElements()) { Error("deplVandLumiExtract","Luminosity vectors do not have the same length"); return EXIT_FAILURE; } double lumi7 = 0.0; double lumi8 = 0.0; double lumi13 = 0.0; double lumi14 = 0.0; std::map<int,double> timeMap = STDatabase::runTimeMap; for (int i=0; i<vr_fillTime.GetNoElements(); i++) { time_t raw_time = (long)vr_fillTime(i); if (vr_beamEner(i)==7.0) { lumi7 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i); }else if(vr_beamEner(i)==8.0){ lumi8 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i); }else if(vr_beamEner(i)==13.0){ lumi13 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i); }else if(vr_beamEner(i)==14.0){ lumi14 += (double)(vr_stopTime(i)-vr_fillTime(i))*vr_peakLumi(i); }else{ Warning("deplVandLumiExtract","Unknown center of mass energy for %s",ctime(&raw_time)); } std::map<int,double>::iterator iter = timeMap.begin(); for (; iter!=timeMap.end(); ++iter) { if (vr_fillTime(i)>iter->second) { v_fill.push_back((double)iter->first); v_lumi7.push_back(lumi7); v_lumi8.push_back(lumi8); v_lumi13.push_back(lumi13); v_lumi14.push_back(lumi14); Info("deplVandLumiExtract","At fill %d %s",iter->first,ctime(&raw_time)); Printf(" Luminosity @ 7 TeV: %7.2f fb^{-1}",lumi7/1e9); Printf(" Luminosity @ 8 TeV: %7.2f fb^{-1}",lumi8/1e9); Printf(" Luminosity @ 13 TeV: %7.2f fb^{-1}",lumi13/1e9); Printf(" Luminosity @ 14 TeV: %7.2f fb^{-1}",lumi14/1e9); v_time.push_back((double)iter->second); timeMap.erase(iter); } } } std::map<int,double>::iterator iter = timeMap.begin(); for (; iter!=timeMap.end(); ++iter) { time_t fill_time = (long)iter->second; v_fill.push_back((double)iter->first); v_lumi7.push_back(lumi7); v_lumi8.push_back(lumi8); v_lumi13.push_back(lumi13); v_lumi14.push_back(lumi14); Info("deplVandLumiExtract","At fill %d: %s",iter->first,ctime(&fill_time)); Printf(" Luminosity @ 7 TeV: %7.2f fb^{-1}",lumi7/1e9); Printf(" Luminosity @ 8 TeV: %7.2f fb^{-1}",lumi8/1e9); Printf(" Luminosity @ 13 TeV: %7.2f fb^{-1}",lumi13/1e9); Printf(" Luminosity @ 14 TeV: %7.2f fb^{-1}",lumi14/1e9); v_time.push_back((double)iter->second); timeMap.erase(iter); } TFile* f_flux = new TFile(fn_flux.Data(),"UPDATE"); TVectorD vr_time = getROOTVector<double>(v_time); TVectorD vr_lumi7 = getROOTVector<double>(v_lumi7); TVectorD vr_lumi8 = getROOTVector<double>(v_lumi8); TVectorD vr_lumi13 = getROOTVector<double>(v_lumi13); TVectorD vr_lumi14 = getROOTVector<double>(v_lumi14); TVectorD vr_fill = getROOTVector<double>(v_fill); f_flux->Delete("v_time;*"); f_flux->Delete("v_lumi7;*"); f_flux->Delete("v_lumi8;*"); f_flux->Delete("v_lumi13;*"); f_flux->Delete("v_lumi14;*"); f_flux->Delete("v_fill;*"); vr_time.Write("v_time"); vr_lumi7.Write("v_lumi7"); vr_lumi8.Write("v_lumi8"); vr_lumi13.Write("v_lumi13"); vr_lumi14.Write("v_lumi14"); vr_fill.Write("v_fill"); f_flux->Close(); Info("deplVandLumiExtract","Results written to %s",fn_flux.Data()); return EXIT_SUCCESS; }