diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0698a9e --- /dev/null +++ b/.gitignore @@ -0,0 +1,41 @@ +# Ignore temporary files +*~ +**/*~ + +# Ignore compiled files +*.pyc +**/*.pyc +*.so +**/*.so +*.o +**/*.o + +# Ignore executables +# Find them with: +# find . \( ! -regex '.*/\..*' \) -type f | xargs -n 1 file | egrep "ELF" +STAging/macros/flukaExtract +STAging/macros/lumiExtract +STAging/macros/tempExtract +macros/runCond/flukaExtract +macros/runCond/fluenceExtract +macros/runCond/lumiExtract +macros/runCond/tempExtract +macros/runCond/deplVandLumiExtrac +macros/CCEScan/pulseFit +macros/CCEScan/landauFit +macros/CCEScan/ratioFit_spline +macros/CCEScan/Hamburg +macros/CCEScan/plotSector +macros/CCEScan/plotSector_gcc3.4.14 +macros/CCEScan/deplVphiN +macros/CCEScan/deplVphiN_normalized +macros/CCEScan/voltageFit_noPulse +macros/CCEScan/voltageFit_spline +macros/CCEScan/ratioFit +macros/CCEScan/voltageFit + +# Ignore temporary scripts +macros/CCEScan/batch/temp_scripts + +# Ignore backups +data/leakageCurrent/summer_2015 diff --git a/STAging/cmt/requirements~ b/STAging/cmt/requirements~ deleted file mode 100755 index 6c0543f..0000000 --- a/STAging/cmt/requirements~ +++ /dev/null @@ -1,44 +0,0 @@ -#============================================================================ -# Created : 2013-05-31 -# Maintainer : Christian elsasser -#============================================================================ -package STAging -version v1r0 - -#============================================================================ -# Structure, i.e. directories to process. -#============================================================================ -branches cmt doc src options python macros data bin -#============================================================================ -# Used packages. Specify the version, * at the end specifies 'any revision' -# Put as many lines as needed, with all packages, without the '#' -#============================================================================ -# use Package v1r* Hat -use GSL v* LCG_Interfaces -use GaudiAlg v* -use GaudiUtils v* -use RecEvent v* Event -use MCEvent v* Event -use TrackEvent v* Event -use PatKernel v* Tf -use TrackFitEvent v* Tr -use TrackInterfaces v* Tr -use TrackKernel v* Tr -use STKernel v* ST -use STTELL1Event v* ST -use PartProp v* Kernel -use CaloDet v* Det -use PhysEvent v* Event -use MCEvent v* Event - -#============================================================================ -# Component library building rule -#============================================================================ -library STAging ../src/*.cpp -import=AIDA -#============================================================================ -# define component library link options -#============================================================================ -apply_pattern install_python_modules -#apply_pattern component_library library=STAging - -alias write_st_xml_cond $STVETRAANALYSISROOT/python/write_st_xml_cond.py \ No newline at end of file diff --git a/STAging/macros/flukaExtract b/STAging/macros/flukaExtract deleted file mode 100755 index 91e015a..0000000 --- a/STAging/macros/flukaExtract +++ /dev/null Binary files differ diff --git a/STAging/macros/flukaExtract.C~ b/STAging/macros/flukaExtract.C~ deleted file mode 100755 index 2a69339..0000000 --- a/STAging/macros/flukaExtract.C~ +++ /dev/null @@ -1,475 +0,0 @@ -// -// flukaExtract.C -// -// -// Created by Christian Elsasser on 31.05.13. -// University of Zurich, elsasser@cern.ch -// Copyright only for commercial use -// - -// General include -#include -#include -#include -#include -#include -#include - - -// ROOT include -#include "incBasicROOT.h" -#include "incDrawROOT.h" -#include "incHistROOT.h" -#include "incIOROOT.h" - - - - -#include "basicROOTTools.h" -#include "basicRooFitTools.h" - -#include "lhcbstyle.C" - - - - - - - - -int flukaExtract(int, char**); - -int writeHisto(TString, TH2F*); -int setUncertainty(TH2F*); -int scale(TH2F*, TH2F*); -int merge(TH2F*, TH2F*); -int storeUncertainty(TH2F*); -int configProfile(TH1D*); - -int main(int argc, char* argv[]){ - - TStyle* style = lhcbStyle(); - - int argc_copy = argc; - - TApplication theApp("Analysis", &argc, argv); - argc = theApp.Argc(); - argv = theApp.Argv(); - - style->SetPadLeftMargin(0.12); - style->SetPadRightMargin(0.22); - style->SetPadBottomMargin(0.15); - style->SetPadTopMargin(0.05); - - int exit = flukaExtract(argc,argv); - - Printf("End"); - theApp.Run(kTRUE); - - return exit; -} - - -// Executable method -int flukaExtract(int argc, char* argv[]){ - - const char* homeVar = std::getenv ("HOME"); - TString dn_lumi(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/fluka",homeVar)); - - TString fn_output(Form("%s/fluka.root", - dn_lumi.Data())); - - - TString fn_TT_old(Form("%s/TTSistrip_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "14TeV")); - TString fn_T3_old(Form("%s/IT3_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "14TeV")); - - - TString fn_TT_8(Form("%s/TTSistrip_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "8TeV")); - TString fn_T3A_8(Form("%s/IT3As_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "8TeV")); - TString fn_T3B_8(Form("%s/IT3bottom_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "8TeV")); - TString fn_T3C_8(Form("%s/IT3Cs_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "8TeV")); - TString fn_T3T_8(Form("%s/IT3top_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "8TeV")); - - TString fn_TT_7(Form("%s/TTSistrip_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "7TeV")); - TString fn_T3A_7(Form("%s/IT3As_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "7TeV")); - TString fn_T3B_7(Form("%s/IT3bottom_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "7TeV")); - TString fn_T3C_7(Form("%s/IT3Cs_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "7TeV")); - TString fn_T3T_7(Form("%s/IT3top_%s_Si1MeVne_cm-2_per_coll.txt", - dn_lumi.Data(), - "7TeV")); - - - - double minRange = -60.0; - double maxRange = 60.0; - double binRange = 1.0; - int nBins = 120; - - TH2F* h_tt_old = new TH2F("h_tt_old","Fluka map", - nBins,minRange,maxRange, - nBins,minRange,maxRange); - h_tt_old->SetXTitle("x position [cm]"); - h_tt_old->SetYTitle("y position [cm]"); - h_tt_old->SetZTitle("#phi_{N,1MeV} [cm^{-2} per coll]"); - h_tt_old->GetZaxis()->SetTitleOffset(1.30); - h_tt_old->SetMinimum(1e-4); - h_tt_old->SetMaximum(1e+0); - - TH2F* h_tt_7 = (TH2F*)h_tt_old->Clone("h_tt_7"); - TH2F* h_tt_8 = (TH2F*)h_tt_old->Clone("h_tt_8"); - TH2F* h_it_old = (TH2F*)h_tt_old->Clone("h_it_old"); - TH2F* h_it_7 = (TH2F*)h_tt_7->Clone("h_it_7"); - TH2F* h_it_8 = (TH2F*)h_tt_8->Clone("h_it_8"); - - writeHisto(fn_TT_old,h_tt_old); - writeHisto(fn_T3_old,h_it_old); - - - - writeHisto(fn_TT_8, h_tt_8); - writeHisto(fn_T3A_8,h_it_8); - writeHisto(fn_T3B_8,h_it_8); - writeHisto(fn_T3C_8,h_it_8); - writeHisto(fn_T3T_8,h_it_8); - - writeHisto(fn_TT_7, h_tt_7); - writeHisto(fn_T3A_7,h_it_7); - writeHisto(fn_T3B_7,h_it_7); - writeHisto(fn_T3C_7,h_it_7); - writeHisto(fn_T3T_7,h_it_7); - - // Mirror old one - int nBinsX = h_tt_old->GetNbinsX(); - int nBinsY = h_tt_old->GetNbinsY(); - for (int i=nBinsX/2+1; i<=nBinsX; i++) { - for (int j=nBinsY/2+1; j<=nBinsY; j++) { - double tt_val = h_tt_old->GetBinContent(i,j); - double it_val = h_it_old->GetBinContent(i,j); - h_tt_old->SetBinContent(nBinsX-i+1,j,tt_val); - h_tt_old->SetBinContent(nBinsX-i+1,nBinsY-j+1,tt_val); - h_tt_old->SetBinContent(i,nBinsY-j+1,tt_val); - h_it_old->SetBinContent(nBinsX-i+1,j,it_val); - h_it_old->SetBinContent(nBinsX-i+1,nBinsY-j+1,it_val); - h_it_old->SetBinContent(i,nBinsY-j+1,it_val); - } - } - - setUncertainty(h_tt_8); - setUncertainty(h_tt_7); - setUncertainty(h_tt_old); - setUncertainty(h_it_8); - setUncertainty(h_it_7); - setUncertainty(h_it_old); - - TH2F* h_tt_old_8 = (TH2F*)h_tt_old->Clone("h_tt_old_8"); - TH2F* h_tt_old_7 = (TH2F*)h_tt_old->Clone("h_tt_old_7"); - TH2F* h_it_old_8 = (TH2F*)h_it_old->Clone("h_it_old_8"); - TH2F* h_it_old_7 = (TH2F*)h_it_old->Clone("h_it_old_7"); - - scale(h_tt_8,h_tt_old_8); - scale(h_it_8,h_it_old_8); - scale(h_tt_7,h_tt_old_7); - scale(h_it_7,h_it_old_7); - - merge(h_tt_8,h_tt_old_8); - merge(h_tt_7,h_tt_old_7); - - storeUncertainty(h_tt_8); - storeUncertainty(h_tt_7); - storeUncertainty(h_it_7); - storeUncertainty(h_it_7); - - - TH1D* h_xProf_7 = h_tt_7->ProjectionX("h_xProf_7", - h_tt_7->GetNbinsY()/2, - h_tt_7->GetNbinsY()/2); - TH1D* h_yProf_7 = h_tt_7->ProjectionY("h_yProf_7", - h_tt_7->GetNbinsX()/2, - h_tt_7->GetNbinsX()/2); - - TH1D* h_xProf_8 = h_tt_8->ProjectionX("h_xProf_8", - h_tt_8->GetNbinsY()/2, - h_tt_8->GetNbinsY()/2); - TH1D* h_yProf_8 = h_tt_8->ProjectionY("h_yProf_8", - h_tt_8->GetNbinsX()/2, - h_tt_8->GetNbinsX()/2); - - configProfile(h_xProf_7); - configProfile(h_yProf_7); - configProfile(h_xProf_8); - configProfile(h_yProf_8); - - - TCanvas* c_tt_8 = new TCanvas("c_tt_8","c_tt_8",1000,800); - c_tt_8->SetLogz(); - h_tt_8->Draw("colz"); - - TCanvas* c_tt_7 = new TCanvas("c_tt_7","c_tt_7",1000,800); - c_tt_7->SetLogz(); - h_tt_7->Draw("colz"); - - TCanvas* c_it_8 = new TCanvas("c_it_8","c_it_8",1000,800); - c_it_8->SetLogz(); - h_it_8->Draw("colz"); - - TCanvas* c_it_7 = new TCanvas("c_it_7","c_it_7",1000,800); - c_it_7->SetLogz(); - h_it_7->Draw("colz"); - - - gStyle->SetPadLeftMargin(0.20); - gStyle->SetPadRightMargin(0.10); - gStyle->SetPadBottomMargin(0.15); - gStyle->SetPadTopMargin(0.05); - - - TCanvas* c_tt_pX7 = new TCanvas("c_tt_pX7","c_tt_pX7",1000,800); - c_tt_pX7->SetLogy(); - h_xProf_7->Draw("H"); - - TCanvas* c_tt_pX8 = new TCanvas("c_tt_pX8","c_tt_pX8",1000,800); - c_tt_pX8->SetLogy(); - h_xProf_8->Draw("H"); - - TCanvas* c_tt_pY7 = new TCanvas("c_tt_pY7","c_tt_pY7",1000,800); - c_tt_pY7->SetLogy(); - h_yProf_7->Draw("H"); - - TCanvas* c_tt_pY8 = new TCanvas("c_tt_pY8","c_tt_pY8",1000,800); - c_tt_pY8->SetLogy(); - h_yProf_8->Draw("H"); - - - TFile* f_output = new TFile(fn_output.Data(),"RECREATE"); - h_tt_7->Write("h_tt_7"); - h_tt_8->Write("h_tt_8"); - - h_it_7->Write("h_t3_7"); - h_it_8->Write("h_t3_8"); - - f_output->Close(); - - - TH1F* h_func = new TH1F("h_func","h_func",60,0.0,60.0); - h_func->SetXTitle("r [cm]"); - h_func->SetYTitle("#simga_{phi_{N,1MeV}}/phi_{N,1MeV}"); - h_func->SetMaximum(1.0); - h_func->SetMinimum(0.0); - - double alpha = 0.918; - double r1 = 0.7e-2; - - TF1* f_relunc = new TF1("f_relunc","[0]*TMath::Power(x,[1])",0.0,100.0); - f_relunc->SetParameters(r1,alpha); - f_relunc->SetLineColor(kRed); - - TCanvas* c_relunc = new TCanvas("c_relunc","c_relunc",1000,800); - h_func->Draw(""); - f_relunc->Draw("Same"); - - - return 0; - -} - - - - - - - -int writeHisto(TString f_name,TH2F* h_write){ - std::ifstream f_input(f_name.Data()); - - std::vector > a_input; - double binWidthX = 0.0; - double binWidthY = 0.0; - while (f_input) { - - - std::string line; - if (!getline(f_input,line)) { - Error("flukaExtract","Problem in line! Leave read-in!"); - break; - } - - - std::string binSizeTag = "* Bin size: "; - - if (line[0]=='*') { - std::size_t p = line.find(binSizeTag); - if (p!=std::string::npos) { - line.replace(line.find(binSizeTag),binSizeTag.length(),""); - std::istringstream s_line(line); - std::string bit; - if (!getline(s_line,bit,',')) { - Error("flukaExtract","Bin width X not available!"); - return 1; - } - binWidthX = atof(bit.c_str()); - if (!getline(s_line,bit,',')) { - Error("flukaExtract","Bin width Y not available!"); - return 1; - } - binWidthY = atof(bit.c_str()); - - Info("flukaExtract","Bin size (%4.2f,%4.2f)",binWidthX,binWidthY); - continue; - } - - Info("flukaExtract","Header line ignoring!"); - continue; - } - - std::istringstream s_line(line); - std::vector v_line; - double x,y,flux; - for (int i=0; i<3 && s_line; i++) { - std::string bit; - if (!getline(s_line,bit,'\t')) { - break; - } - if (bit.compare("")==0) { - break; - } - if (i==0) { - x = atof(bit.c_str()); - }else if(i==1){ - y = atof(bit.c_str()); - }else if(i==2){ - flux = atof(bit.c_str()); - for (int xcoord=0; 1.0*xcoordFindBin(x+0.5+1.0*xcoord,y+0.5+1.0*ycoord); - h_write->SetBinContent(bin,flux); - } - } - } - - - } - - } - return 0; -} - -// Determines uncertainty on fluence -int setUncertainty(TH2F* h_write){ - // Determination of uncertainty - double alpha = 0.918; - double r1 = 0.7e-2; - for (int i=1; i<=h_write->GetNbinsX(); i++) { - for (int j=1; j<=h_write->GetNbinsY(); j++) { - double x = h_write->GetXaxis()->GetBinCenter(i); - double y = h_write->GetYaxis()->GetBinCenter(j); - double r = TMath::Sqrt(x*x+y*y); - double u = r1*TMath::Power(r,alpha); - h_write->SetBinError(i,j,u); - } - } - return 0; -} - -// Scales old FLUKA map by chi^2 fit -int scale(TH2F* h_ref, TH2F* h_tar){ - if (h_tar->GetNbinsX()!=h_ref->GetNbinsX() || - h_tar->GetNbinsY()!=h_ref->GetNbinsY()) { - Error("scale","Histograms do not have the same size!"); - return 1; - } - double xy = 0.0, xx = 0.0; - for (int i=1; i<=h_tar->GetNbinsX(); i++) { - for (int j=1; j<=h_tar->GetNbinsY(); j++) { - if (TMath::Abs(h_tar->GetBinContent(i,j))<1e-9 || - TMath::Abs(h_ref->GetBinContent(i,j))<1e-9) { - continue; - } - xy+=h_tar->GetBinContent(i,j)/(TMath::Power(h_ref->GetBinError(i,j),2)*h_ref->GetBinContent(i,j)); - xx+=TMath::Power(h_tar->GetBinContent(i,j),2)/(TMath::Power(h_ref->GetBinError(i,j)*h_ref->GetBinContent(i,j),2)); - } - } - double alpha = xy/xx; - double chi2 = 0.0; - int ndf = 0; - for (int i=1; i<=h_tar->GetNbinsX(); i++) { - for (int j=1; j<=h_tar->GetNbinsY(); j++) { - h_tar->SetBinContent(i,j,alpha*h_tar->GetBinContent(i,j)); - if (TMath::Abs(h_tar->GetBinContent(i,j))<1e-9 || - TMath::Abs(h_ref->GetBinContent(i,j))<1e-9) { - continue; - } - chi2 += TMath::Power((h_tar->GetBinContent(i,j)-h_ref->GetBinContent(i,j))/ - (h_ref->GetBinContent(i,j)*h_ref->GetBinError(i,j)),2); - ndf++; - } - } - Info("scale","alpha = %4.2f",alpha); - Info("scale","chi^2 = %4.2f/%d",chi2,(ndf-1)); - return 0; -} - -// Adds values for empty bins from old FLUKA map -int merge(TH2F* h_tar, TH2F* h_sou){ - if (h_tar->GetNbinsX()!=h_sou->GetNbinsX() || - h_tar->GetNbinsY()!=h_sou->GetNbinsY()) { - Error("merge","Histograms do not have the same size!"); - return 1; - } - for (int i=1; i<=h_tar->GetNbinsX(); i++) { - for (int j=1; j<=h_tar->GetNbinsY(); j++) { - if (TMath::Abs(h_tar->GetBinContent(i,j))>1e-9) { - continue; - } - h_tar->SetBinContent(i,j,h_sou->GetBinContent(i,j)); - h_tar->SetBinError(i,j,h_sou->GetBinError(i,j)); - } - } - return 0; -} - -// Computes absolute uncertainty from relative uncertainty -int storeUncertainty(TH2F* h){ - for (int i=1; i<=h->GetNbinsX(); i++) { - for (int j=1; j<=h->GetNbinsY(); j++) { - h->SetBinError(i,j,h->GetBinError(i,j)*h->GetBinContent(i,j)); - } - } - return 0; -} - -// Configure Profiles -int configProfile(TH1D* h){ - h->SetYTitle("#phi_{N,1MeV} [cm^{-2} per coll]"); - h->SetLineColor(kRed); - h->SetMarkerColor(kRed); - h->GetYaxis()->SetTitleOffset(1.3); - h->SetMaximum(1.e+0); - h->SetMinimum(1.e-4); - return 0; -} - - diff --git a/STAging/macros/lumiExtract b/STAging/macros/lumiExtract deleted file mode 100755 index 7ba6a54..0000000 --- a/STAging/macros/lumiExtract +++ /dev/null Binary files differ diff --git a/STAging/macros/tempExtract b/STAging/macros/tempExtract deleted file mode 100755 index ba24f51..0000000 --- a/STAging/macros/tempExtract +++ /dev/null Binary files differ diff --git a/STAging/options/base/HitIntegrator.py~ b/STAging/options/base/HitIntegrator.py~ deleted file mode 100755 index b243cbd..0000000 --- a/STAging/options/base/HitIntegrator.py~ +++ /dev/null @@ -1,87 +0,0 @@ -""" -############################################################################## -# # -# Hit integration over ROI per sector # -# # -# @package Tell1/Vetra # -# @author Christan Elsasser # -# @date 09/06/2012 # -# # -############################################################################## -""" - -from GaudiKernel.ProcessJobOptions import importOptions - -from Configurables import GaudiSequencer -from TrackSys.Configuration import * -from Configurables import ( TrackSys, RecSysConf ) -from Configurables import TrackContainerCopy -from Configurables import Velo__VeloTrackMonitor -from Configurables import EventClockSvc -from Configurables import Vetra -from Configurables import VetraRecoConf -from Configurables import ST__STDumpADCs - - - -generalOutputLevel = 2 - -Vetra().MainSequence = [ GaudiSequencer('HitIntegrator')] -#GaudiSequencer('VetraSequencer').IgnoreFilterPassed = True - -## TrackSys().ExcludedLayers = ["TTbX"] -## TrackSys().TrackPatRecAlgorithms = TrackSys().DefaultPatRecAlgorithms -## VetraRecoConf().Sequence = RecSysConf().DefaultTrackingSubdets # = ["Decoding", "VELO","TT","IT","OT","Tr","Vertex"] - -## from Configurables import CondDB, CondDBAccessSvc, CondDBTimeSwitchSvc -## connection="sqlite_file:$STSQLDDDBROOT/db/STCOND.db/COND" -## CondDB().addLayer( CondDBAccessSvc("COND", ConnectionString = connection) ) - - - -## importOptions('$HOME/cmtuser/Vetra_v12r2/ST/STVetraAnalysis/options/TTEmulator.py') -#importOptions('$HOME/cmtuser/Vetra_v12r2/ST/STVetraAnalysis/options/ITEmulator.py') - - -from Configurables import STHitIntegrator - -HitIntegrator = STHitIntegrator() -HitIntegrator.DetType = "TT" -HitIntegrator.Layer = "TTaU" -HitIntegrator.MaxX = 800.0 -HitIntegrator.MaxY = 800.0 -HitIntegrator.OutputLevel = 3 - - -#GaudiSequencer('HitIntegrator').Members = [ GaudiSequencer('EmuTTSeq'), GaudiSequencer('EmuTTSeqPrev'), GaudiSequencer('EmuTTSeqNext'),HitIntegrator ] -#GaudiSequencer('HitIntegrator').IgnoreFilterPassed = True#ReturnOK = True - - -GaudiSequencer('HitIntegrator').Members = [HitIntegrator ] -GaudiSequencer('HitIntegrator').IgnoreFilterPassed = True#ReturnOK = True - - - -#path = "/afs/cern.ch/work/e/elsasser/Aging/TT/TTbX" -path = "." -histoname = "Histo.root" -tuplename = "Tuple.root" - -Vetra().EvtMax = 0 -#Vetra().PrintFreq = 100 -Vetra().OutputLevel = INFO -Vetra().HistogramFile = histoname -Vetra().TupleFile = tuplename -vetra = Vetra() - -def doMyChanges(): - from Configurables import VetraInit, FilterByBankType - bankFilter = FilterByBankType() - bankFilter.BankNames = [ "TTFull", "ITFull" ] - from Configurables import LoKi__ODINFilter - odinFilter = LoKi__ODINFilter('BBODINFilter', Code = " LHCb.ODIN.BeamCrossing == ODIN_BXTYP ") - GaudiSequencer( 'Init' ).Members = [ VetraInit(), odinFilter, bankFilter ] -appendPostConfigAction(doMyChanges) - - - diff --git a/STAging/options/base/ITAging.py~ b/STAging/options/base/ITAging.py~ deleted file mode 100755 index 1fe17db..0000000 --- a/STAging/options/base/ITAging.py~ +++ /dev/null @@ -1,85 +0,0 @@ -""" -############################################################################## -# # -# Configuration for ST Aging CCE scan Analys for IT # -# # -# @package ST/STAging # -# @author Ch. Elsasser # -# @date 15/06/2013 # -# # -############################################################################## -""" - -from GaudiKernel.ProcessJobOptions import importOptions - -from Configurables import GaudiSequencer -from TrackSys.Configuration import * -from Configurables import ( TrackSys, RecSysConf ) -from Configurables import TrackContainerCopy -from Configurables import Velo__VeloTrackMonitor -from Configurables import EventClockSvc -from Configurables import Vetra -from Configurables import VetraRecoConf -from Configurables import ST__STDumpADCs - - - -generalOutputLevel = 3 - -Vetra().MainSequence = [ 'ProcessPhase/Reco', GaudiSequencer('MoniITSeq'), GaudiSequencer('MoniSTSeq')] -#GaudiSequencer('VetraSequencer').IgnoreFilterPassed = True - -TrackSys().ExcludedLayers = ["T3X2"] -TrackSys().TrackPatRecAlgorithms = TrackSys().DefaultPatRecAlgorithms -VetraRecoConf().Sequence = RecSysConf().DefaultTrackingSubdets # = ["Decoding", "VELO","TT","IT","OT","Tr","Vertex"] - -from Configurables import CondDB, CondDBAccessSvc, CondDBTimeSwitchSvc -connection="sqlite_file:$STSQLDDDBROOT/db/STCOND.db/COND" -CondDB().addLayer( CondDBAccessSvc("COND", ConnectionString = connection) ) - -importOptions('$HOME/cmtuser/Vetra_v13r2/ST/STVetraAnalysis/options/ITEmulator.py') - - -from Configurables import STNZSResolution, STADCTrackMonitor - -#ODINChecker = STODINCheck() -#ODINChecker.OutputLevel = generalOutputLevel -#ODINChecker.ODINData = "DAQ/ODIN"; - -#GaudiSequencer('ODINPreChecker').Members = [ODINChecker] - -itRes = STNZSResolution("ITResolution") -itRes.OutputLevel = 3 -itRes.DetType = "IT" -itRes.UseNZSdata = True -itRes.InputData = "/Event/Raw/IT/LCMSADCs"; - -TrackMonitor = STADCTrackMonitor() -TrackMonitor.OutputLevel = generalOutputLevel -TrackMonitor.DetType = "IT" -TrackMonitor.RawDataCent = "/Event/Raw/IT/LCMSADCs"; -TrackMonitor.ODINData = "DAQ/ODIN"; -TrackMonitor.Layer = "T3X2" -TrackMonitor.MaxNeighbours = 3 - -GaudiSequencer('MoniITSeq').Members = [ GaudiSequencer('EmuITSeq'), GaudiSequencer('EmuITSeqPrev'), GaudiSequencer('EmuITSeqNext'),TrackMonitor] -GaudiSequencer('MoniITSeq').IgnoreFilterPassed = True#ReturnOK = True - -histoname = "Histo.root" -tuplename = "Tuple.root" - -Vetra().EvtMax = -1 -Vetra().OutputLevel = INFO -Vetra().HistogramFile = histoname -Vetra().TupleFile = tuplename -vetra = Vetra() - - -def doMyChanges(): - from Configurables import VetraInit, FilterByBankType - bankFilter = FilterByBankType() - bankFilter.BankNames = [ "TTFull", "ITFull" ] - from Configurables import LoKi__ODINFilter - odinFilter = LoKi__ODINFilter('BBODINFilter', Code = " LHCb.ODIN.BeamCrossing == ODIN_BXTYP ") - GaudiSequencer( 'Init' ).Members = [ VetraInit(), odinFilter, bankFilter ] -appendPostConfigAction(doMyChanges) diff --git a/STAging/options/base/TTAging.py~ b/STAging/options/base/TTAging.py~ deleted file mode 100755 index d0e66b4..0000000 --- a/STAging/options/base/TTAging.py~ +++ /dev/null @@ -1,86 +0,0 @@ -""" -############################################################################## -# # -# Configuration for ST Aging CCE scan analysis for TT # -# # -# @package ST/STAging # -# @author Ch. Elsasser # -# @date 15/06/2013 # -# # -############################################################################## -""" - -from GaudiKernel.ProcessJobOptions import importOptions - -from Configurables import GaudiSequencer -from TrackSys.Configuration import * -from Configurables import ( TrackSys, RecSysConf ) -from Configurables import TrackContainerCopy -from Configurables import Velo__VeloTrackMonitor -from Configurables import EventClockSvc -from Configurables import Vetra -from Configurables import VetraRecoConf -from Configurables import ST__STDumpADCs - - - -generalOutputLevel = 3 - -Vetra().MainSequence = [ 'ProcessPhase/Reco', GaudiSequencer('MoniTTSeq'), GaudiSequencer('MoniSTSeq')] - - -TrackSys().ExcludedLayers = ["TTaU"] -TrackSys().TrackPatRecAlgorithms = TrackSys().DefaultPatRecAlgorithms -VetraRecoConf().Sequence = RecSysConf().DefaultTrackingSubdets # = ["Decoding", "VELO","TT","IT","OT","Tr","Vertex"] - -from Configurables import CondDB, CondDBAccessSvc, CondDBTimeSwitchSvc -connection="sqlite_file:$STSQLDDDBROOT/db/STCOND.db/COND" -CondDB().addLayer( CondDBAccessSvc("COND", ConnectionString = connection) ) - - -importOptions('$HOME/cmtuser/Vetra_v13r2/ST/STVetraAnalysis/options/TTEmulator.py') - -from Configurables import STNZSResolution, STADCTrackMonitor - -#ODINChecker = STODINCheck() -#ODINChecker.OutputLevel = generalOutputLevel -#ODINChecker.ODINData = "DAQ/ODIN"; - -#GaudiSequencer('ODINPreChecker').Members = [ODINChecker] - -ttRes = STNZSResolution("TTResolution") -ttRes.OutputLevel = 4 -ttRes.DetType = "TT" -ttRes.UseNZSdata = True -ttRes.InputData = "/Event/Raw/TT/LCMSADCs"; - -TrackMonitor = STADCTrackMonitor() -TrackMonitor.OutputLevel = generalOutputLevel -TrackMonitor.DetType = "TT" -TrackMonitor.RawDataCent = "/Event/Raw/TT/LCMSADCs"; -TrackMonitor.ODINData = "DAQ/ODIN"; -TrackMonitor.Layer = "TTaU" -TrackMonitor.MaxNeighbours = 3 - - - -GaudiSequencer('MoniTTSeq').Members = [ GaudiSequencer('EmuTTSeq'), GaudiSequencer('EmuTTSeqPrev'), GaudiSequencer('EmuTTSeqNext'),TrackMonitor] -GaudiSequencer('MoniTTSeq').IgnoreFilterPassed = True#ReturnOK = True - -histoname = "Histo.root" -tuplename = "Tuple.root" - -Vetra().EvtMax = -1 -Vetra().OutputLevel = INFO -Vetra().HistogramFile = histoname -Vetra().TupleFile = tuplename -vetra = Vetra() - -def doMyChanges(): - from Configurables import VetraInit, FilterByBankType - bankFilter = FilterByBankType() - bankFilter.BankNames = [ "TTFull", "ITFull" ] - from Configurables import LoKi__ODINFilter - odinFilter = LoKi__ODINFilter('BBODINFilter', Code = " LHCb.ODIN.BeamCrossing == ODIN_BXTYP ") - GaudiSequencer( 'Init' ).Members = [ VetraInit(), odinFilter, bankFilter ] -appendPostConfigAction(doMyChanges) diff --git a/STAging/options/base/TTClusters.py~ b/STAging/options/base/TTClusters.py~ deleted file mode 100755 index 2cef0ea..0000000 --- a/STAging/options/base/TTClusters.py~ +++ /dev/null @@ -1,78 +0,0 @@ - """ -############################################################################## -# # -# Configuration for ST Cluster analysis for TT # -# # -# @package ST/STAging # -# @author Ch. Elsasser # -# @date 15/06/2013 # -# # -############################################################################## -""" - -from GaudiKernel.ProcessJobOptions import importOptions - -from Configurables import GaudiSequencer -from TrackSys.Configuration import * -from Configurables import ( TrackSys, RecSysConf ) -from Configurables import TrackContainerCopy -from Configurables import Velo__VeloTrackMonitor -from Configurables import EventClockSvc -from Configurables import Vetra -from Configurables import VetraRecoConf -from Configurables import ST__STDumpADCs - - - -generalOutputLevel = 3 - -Vetra().MainSequence = [ 'ProcessPhase/Reco', GaudiSequencer('MoniTTSeq'), GaudiSequencer('MoniSTSeq')] - - - -TrackSys().TrackPatRecAlgorithms = TrackSys().DefaultPatRecAlgorithms -VetraRecoConf().Sequence = RecSysConf().DefaultTrackingSubdets # = ["Decoding", "VELO","TT","IT","OT","Tr","Vertex"] - -from Configurables import CondDB, CondDBAccessSvc, CondDBTimeSwitchSvc -connection="sqlite_file:$STSQLDDDBROOT/db/STCOND.db/COND" -CondDB().addLayer( CondDBAccessSvc("COND", ConnectionString = connection) ) - - -importOptions('$HOME/cmtuser/Vetra_v13r2/ST/STVetraAnalysis/options/TTEmulator.py') - -from Configurables import STNZSResolution, STClusterTrackMonitor - -#ODINChecker = STODINCheck() -#ODINChecker.OutputLevel = generalOutputLevel -#ODINChecker.ODINData = "DAQ/ODIN"; - -#GaudiSequencer('ODINPreChecker').Members = [ODINChecker] - -TrackMonitor = STClusterTrackMonitor() -TrackMonitor.OutputLevel = 4 -TrackMonitor.DetType = "TT" -TrackMonitor.ODINData = "DAQ/ODIN"; - - - - -GaudiSequencer('MoniTTSeq').Members = [ GaudiSequencer('EmuTTSeq'), GaudiSequencer('EmuTTSeqPrev'), GaudiSequencer('EmuTTSeqNext'),TrackMonitor] -GaudiSequencer('MoniTTSeq').IgnoreFilterPassed = True#ReturnOK = True - -histoname = "Histo.root" -tuplename = "Tuple.root" - -Vetra().EvtMax = 8000 -Vetra().OutputLevel = INFO -Vetra().HistogramFile = histoname -Vetra().TupleFile = tuplename -vetra = Vetra() - -def doMyChanges(): - from Configurables import VetraInit, FilterByBankType - bankFilter = FilterByBankType() - bankFilter.BankNames = [ "TTFull", "ITFull" ] - from Configurables import LoKi__ODINFilter - odinFilter = LoKi__ODINFilter('BBODINFilter', Code = " LHCb.ODIN.BeamCrossing == ODIN_BXTYP ") - GaudiSequencer( 'Init' ).Members = [ VetraInit(), odinFilter, bankFilter ] -appendPostConfigAction(doMyChanges) diff --git a/STAging/options/base/filetest.py~ b/STAging/options/base/filetest.py~ deleted file mode 100755 index 0cbaaef..0000000 --- a/STAging/options/base/filetest.py~ +++ /dev/null @@ -1,10 +0,0 @@ -from Gaudi.Configuration import * -from Configurables import LHCbApp -EventSelector().PrintFreq = 100 -#LHCbApp().SkipEvents = 3000 -EventSelector().Input += [ - "DATAFILE='/tmp/elsasser/ST.raw' SVC='LHCb::MDFSelector'" - ] -ApplicationMgr().HistogramPersistency = 'ROOT' -HistogramPersistencySvc().OutputFile = 'run-69333.root' -Vetra().EvtMax = 1000 diff --git a/STAging/options/ganga/ITsubmit.py~ b/STAging/options/ganga/ITsubmit.py~ deleted file mode 100755 index da174c4..0000000 --- a/STAging/options/ganga/ITsubmit.py~ +++ /dev/null @@ -1,69 +0,0 @@ -""" -############################################################################## -# # -# Ganga submission script to extract CCE scans for IT # -# # -# @package ST/STAging # -# @author Ch. Elsasser # -# @date 27/06/2013 # -# # -############################################################################## -""" - - - -# List runs belonging to fills - -runDict = { - 3478 : [135603, 135604, 135605], # 2013-01-21 - 3108 : [129494, 129518, 129519], # 2012-09-28 - 2797 : [120006, 120009, 120010], # 2012-07-02 - 2472 : [111267, 111270, 111271], # 2012-04-05 - 2252 : [104090, 104091, 104092, 104093, 104106, 104107], # 2011-10-25 - 2083 : [101357, 101359], # 2011-09-07 - 1944 : [ 95930, 95932, 95933, 95934, 95936], # 2011-07-14 - 1616 : [ 87144, 87146, 87147, 87151, 87152, 87153], # 2011-03-14 - 1020 : [ 69505, 69506, 69507, 69508, 69509, 69510, - 69511, 69513, 69517, 69521, 69525, 69526, - 69528, 69531, 69555, 69556, 69558, 69570, - 69572, 69573, 69574, 69575, 69576, 69577, - 69578, 69579, 69580, 69583, 69586, - 69589, 69592, 69593, 69594, 69595, 69596, - 69597, 69599, 69602, 69603, 69604, 69605, - 69606, 69607, 69608, 69610, 69611], # 2010-04-04 - 0 : [129494] -} - -fills=[1616,1944,2083,2252,2472,2797,3108,3478] - -#fills=[3108,3478] -#fills=[1020] -#fills=[0] - -for k in fills: - runs=runDict[k] - print "=====================================" - print "Fill: "+str(k) - for i in runs: - print "=====================================" - j=Job(application='Vetra',name="IT-"+str(i)) - print "Run: "+str(i) - j.application.optsfile="/afs/cern.ch/user/e/elsasser/cmtuser/Vetra_v13r2/ST/STAging/options/base/ITAging.py" - j.application.version='v15r0' - bkq = BKQuery(type='Run', dqflag='All') - bkq.path = "/"+str(i)+"/Real Data/90000001/RAW" - print "BK Query path: "+bkq.path - j.inputdata = bkq.getDataset() - print "Input data contain "+str(len(j.inputdata.getFileNames()))+" file(s)" - j.backend=Dirac() - j.backend.settings['CPUTime']=2*172800 - j.splitter=SplitByFiles(filesPerJob=5, - ignoremissing=True, - bulksubmit=True) - j.outputfiles=[DiracFile('*.root')] - j.postprocessors=[Notifier(address='elsasser@cern.ch')] - j.submit() - print "Job submitted!" - print "=====================================" - print "" - print "" diff --git a/STAging/options/ganga/TTClusterSubmit.py~ b/STAging/options/ganga/TTClusterSubmit.py~ deleted file mode 100755 index c42b2d3..0000000 --- a/STAging/options/ganga/TTClusterSubmit.py~ +++ /dev/null @@ -1,51 +0,0 @@ -""" -############################################################################## -# # -# Ganga submission script to analyze clusters for TT # -# # -# @package ST/STAging # -# @author Ch. Elsasser # -# @date 27/06/2013 # -# # -############################################################################## -""" - - - -# List runs belonging to fills - -runDict = { - 0 : [117200] -} - -#fills=[1020,1616,1944,2083,2252,2472,2797,3108,3478] -fills=[0] -#fills=[1020] - -for k in fills: - runs=runDict[k] - print "=====================================" - print "Fill: "+str(k) - for i in runs: - print "=====================================" - j=Job(application='Vetra',name="TT-"+str(i)) - print "Run: "+str(i) - j.application.optsfile="/afs/cern.ch/user/e/elsasser/cmtuser/Vetra_v13r2/ST/STAging/options/base/TTClusters.py" - j.application.version='v13r2' - bkq = BKQuery(type='Run', dqflag='All') - bkq.path = "/"+str(i)+"/Real Data/96000000/RAW" - print "BK Query path: "+bkq.path - j.inputdata = bkq.getDataset() - print "Input data contain "+str(len(j.inputdata.getFileNames()))+" file(s)" - j.backend=Dirac() - j.backend.settings['CPUTime']=2*172800 - j.splitter=SplitByFiles(filesPerJob=5, - ignoremissing=True, - bulksubmit=True) - j.outputfiles=[DiracFile('*.root')] - j.postprocessors=[Notifier(address='elsasser@cern.ch')] - j.submit() - print "Job submitted!" - print "=====================================" - print "" - print "" diff --git a/STAging/options/ganga/TTsubmit.py~ b/STAging/options/ganga/TTsubmit.py~ deleted file mode 100755 index bb973d0..0000000 --- a/STAging/options/ganga/TTsubmit.py~ +++ /dev/null @@ -1,67 +0,0 @@ -""" -############################################################################## -# # -# Ganga submission script to extract CCE scans for TT # -# # -# @package ST/STAging # -# @author Ch. Elsasser # -# @date 27/06/2013 # -# # -############################################################################## -""" - - - -# List runs belonging to fills - -runDict = { - 3478 : [135603, 135604, 135605], # 2013-01-21 - 3108 : [129494, 129518, 129519], # 2012-09-28 - 2797 : [120006, 120009, 120010], # 2012-07-02 - 2472 : [111267, 111270, 111271], # 2012-04-05 - 2252 : [104090, 104091, 104092, 104093, 104106, 104107], # 2011-10-25 - 2083 : [101357, 101359], # 2011-09-07 - 1944 : [ 95930, 95932, 95933, 95934, 95936], # 2011-07-14 - 1616 : [ 87144, 87146, 87147, 87151, 87152, 87153], # 2011-03-14 - 1020 : [ 69505, 69506, 69507, 69508, 69509, 69510, - 69511, 69513, 69517, 69521, 69525, 69526, - 69528, 69531, 69555, 69556, 69558, 69570, - 69572, 69573, 69574, 69575, 69576, 69577, - 69578, 69579, 69580, 69583, 69586, - 69589, 69592, 69593, 69594, 69595, 69596, - 69597, 69599, 69602, 69603, 69604, 69605, - 69606, 69607, 69608, 69610, 69611], # 2010-04-04 - 0 : [129494] -} - -fills=[1020,1616,1944,2083,2252,2472,2797,3108,3478] -#fills=[3108,3478] -#fills=[1020] - -for k in fills: - runs=runDict[k] - print "=====================================" - print "Fill: "+str(k) - for i in runs: - print "=====================================" - j=Job(application='Vetra',name="TT-"+str(i)) - print "Run: "+str(i) - j.application.optsfile="/afs/cern.ch/user/e/elsasser/cmtuser/Vetra_v13r2/ST/STAging/options/base/TTAging.py" - j.application.version='v15r0' - bkq = BKQuery(type='Run', dqflag='All') - bkq.path = "/"+str(i)+"/Real Data/90000001/RAW" - print "BK Query path: "+bkq.path - j.inputdata = bkq.getDataset() - print "Input data contain "+str(len(j.inputdata.getFileNames()))+" file(s)" - j.backend=Dirac() - j.backend.settings['CPUTime']=2*172800 - j.splitter=SplitByFiles(filesPerJob=5, - ignoremissing=True, - bulksubmit=True) - j.outputfiles=[DiracFile('*.root')] - j.postprocessors=[Notifier(address='elsasser@cern.ch')] - j.submit() - print "Job submitted!" - print "=====================================" - print "" - print "" diff --git a/STAging/src/STADCTrackMonitor.cpp~ b/STAging/src/STADCTrackMonitor.cpp~ deleted file mode 100755 index 1046d3b..0000000 --- a/STAging/src/STADCTrackMonitor.cpp~ +++ /dev/null @@ -1,610 +0,0 @@ -// Include files -#include "STADCTrackMonitor.h" - -//event -#include "Event/Track.h" -#include "Event/TrackFitResult.h" -#include "Event/State.h" -#include "Kernel/LHCbID.h" -#include "Event/ODIN.h" -#include "Event/VeloCluster.h" -#include "Event/VeloPhiMeasurement.h" -#include "Event/VeloRMeasurement.h" -#include "Event/KalmanFitResult.h" -#include "Kernel/HitPattern.h" -#include "TrackInterfaces/IHitExpectation.h" -#include "TrackInterfaces/IVeloExpectation.h" - -// Det -#include "STDet/DeSTDetector.h" -#include "STDet/DeITDetector.h" -#include "VeloDet/DeVelo.h" -#include "Event/FitNode.h" -#include "Event/RecVertex.h" - - -// ST DAQ -#include "Kernel/STDAQDefinitions.h" -#include "Kernel/STBoardMapping.h" -#include "Kernel/LHCbConstants.h" -#include "Kernel/ISTReadoutTool.h" -#include "Kernel/STDAQDefinitions.h" - -// STTELL1Event -#include "Event/STTELL1Data.h" -#include "Kernel/STTell1Board.h" - -// gsl -#include "gsl/gsl_math.h" - -#include "GaudiKernel/AlgFactory.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/ToStream.h" - -#include -#include -#include -#include "Map.h" - -// Boost -#include -#include -#include -#include - -using namespace boost::lambda; -using namespace LHCb; -using namespace Gaudi; - - -DECLARE_ALGORITHM_FACTORY( STADCTrackMonitor ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -STADCTrackMonitor::STADCTrackMonitor(const std::string& name, - ISvcLocator* pSvcLocator ) : - TrackMonitorTupleBase( name , pSvcLocator ) -{ - declareProperty("DetType",m_detType,"TT"); - declareProperty("Layer",m_layer,"TTaX"); - declareProperty("RawDataCent",m_rawDataLocCent,"/Event/Raw/TT/LCMSADCs"); - - declareProperty("ODINData",m_odinLoc,"DAQ/ODIN"); - - declareProperty("MaxNeighbours",m_maxNeighbours=3); - - - // TT layer mapping - m_layerMap.insert(std::pair("TTaX",9)); - m_layerMap.insert(std::pair("TTaU",10)); - m_layerMap.insert(std::pair("TTbV",17)); - m_layerMap.insert(std::pair("TTbX",18)); - - // IT layer mapping - m_layerMap.insert(std::pair("T1X1",9)); - m_layerMap.insert(std::pair("T1U", 10)); - m_layerMap.insert(std::pair("T1V", 11)); - m_layerMap.insert(std::pair("T1X2",12)); - - m_layerMap.insert(std::pair("T2X1",17)); - m_layerMap.insert(std::pair("T2U", 18)); - m_layerMap.insert(std::pair("T2V", 19)); - m_layerMap.insert(std::pair("T2X2",20)); - - m_layerMap.insert(std::pair("T3X1",25)); - m_layerMap.insert(std::pair("T3U", 26)); - m_layerMap.insert(std::pair("T3V", 27)); - m_layerMap.insert(std::pair("T3X2",28)); -} - -//============================================================================= -// Destructor -//============================================================================= -STADCTrackMonitor::~STADCTrackMonitor() {}; - -StatusCode STADCTrackMonitor::initialize() -{ - // Mandatory initialization of GaudiAlgorith - StatusCode sc = TrackMonitorTupleBase::initialize(); - - - - - if ( sc.isFailure() ) { return sc; } - m_veloDet = getDet( DeVeloLocation::Default ) ; - - - // Definition of readout tool - if ( m_detType=="IT" ){ - m_readoutTool = tool("ITReadoutTool"); - }else{ - m_readoutTool = tool("TTReadoutTool"); - } - - - m_veloExpectation = tool("VeloExpectation"); - m_ttExpectation = tool("TTHitExpectation"); - m_itExpectation = tool("ITHitExpectation"); - m_otExpectation = tool("OTHitExpectation"); - - m_extrapolator = tool("TrackParabolicExtrapolator"); - m_STDet = getDet(DeSTDetLocation::location(m_detType)); - - - return StatusCode::SUCCESS; -} - -//============================================================================= -// Execute -//============================================================================= -StatusCode STADCTrackMonitor::execute() -{ - - - // Create Tuple - Tuples::Tuple m_tuple_hit = GaudiTupleAlg::nTuple("HitInfo/"+m_layer, "HitInfo", CLID_ColumnWiseTuple ); - Tuples::Tuple m_tuple_noi = GaudiTupleAlg::nTuple("NoiseInfo/"+m_layer, "NoiseInfo", CLID_ColumnWiseTuple ); - - // Check if ODIN is available, otherwise skip event - bool isEverythingHere = true; - - if(!exist(m_odinLoc)){ - warning() << "ODIN data are not available!" << endmsg; - isEverythingHere = false; - } - - if(!isEverythingHere){ - warning() << "Skip event!" << endmsg; - return StatusCode::SUCCESS; - } - - // get ODIN - const LHCb::ODIN* odin = get(m_odinLoc); - int calibStepValue = odin->calibrationStep(); - debug() << "Value of ODIN calibration step is " << calibStepValue << "." << endmsg; - - // Check if track container is available,... - if(!exist(inputContainer())){ - warning() << "Tracks are not available!" << endmsg; - isEverythingHere = false; - } - - // Check if data of this bunch crossing is available,... - if(!exist(m_rawDataLocCent)){ - warning() << "Central raw data are not available!" << endmsg; - isEverythingHere = false; - } - - // ... otherwise skip event, unless next/previous bunch crossing is excluded - if(!isEverythingHere){ - warning() << "Skip event!" << endmsg; - return StatusCode::SUCCESS; - } - - - - // get the input data - LHCb::Track::Range tracks = get(inputContainer()); - - // define data prefix - std::vector dataPrefix; - - // get raw data - std::vector rawData; - - // get raw data locations - rawData.push_back(get(m_rawDataLocCent)); - dataPrefix.push_back("C"); - - - // Cast calibration step to string - std::string calibStep = boost::lexical_cast( calibStepValue ); - - //for(MultiplicityMap::iterator it = m_multiplicityMap.begin() ; - // it != m_multiplicityMap.end(); ++it) it->second = 0 ; - - - // Extract further event info - int runN = odin->runNumber(); - int evtN = odin->eventNumber(); - int BCID = odin->bunchId(); - int BCType = odin->bunchCrossingType(); - - - // Vectors to store boards and strip assigned to tracks - std::vector boardVec; - std::vector stripVec; - - - // Loop over tracks in event - BOOST_FOREACH( const LHCb::Track* track, tracks) { - std::string type = histoDirName(*track) ; - m_multiplicityMap[type] += 1; - - // Extract track momentum - double trackMom = track->p(); - - - debug() << "Track type is " << track->type() << " with momentum " << trackMom << " MeV/c." << endmsg; - - // Check if track is long (type()=3) TO DO: change to Track::Long - if(track->type()==Track::Long){ - debug() << "It's a long track!" << endmsg; - - // Extract track info - double px = track->momentum().x(); - double py = track->momentum().y(); - double pz = track->momentum().z(); - double pt = track->pt(); - double trChi2 = track->chi2(); - int ndof = track->nDoF(); - double ghostP = track->ghostProbability(); - int nMeas = track->nMeasurements(); - int nIDs = track->nLHCbIDs(); - - - - - // Get extrapolated hit position - StateVector* hit = 0; - if(m_detType=="TT"){ - hit = GetTTHit(track); - }else{ - hit = GetITHit(track); - } - - if(!hit){ - warning() << "Extrapolated hit not found! Skip event!" << endmsg; - return StatusCode::SUCCESS; - } - - // Extract hit - double xHit = hit->x(); - double yHit = hit->y(); - double zHit = hit->z(); - double dxdz = hit->tx(); - double dydz = hit->ty(); - double pHit = hit->p(); - - // Get hits extrapolated from track - std::vector ids(0); - if(m_detType=="TT"){ - m_ttExpectation->collect(*track,ids); - }else{ - m_itExpectation->collect(*track,ids); - } - - delete hit; - - std::vector< LHCb::STChannelID > flt_ids; - - // Check if hit is in right detector and layer - BOOST_FOREACH( LHCb::LHCbID chan, ids) { - if (chan.isST() == true){ - if (m_detType == "IT" ){ - if (chan.isIT() == true && m_layerMap[m_layer]==chan.stID().uniqueLayer()) flt_ids.push_back(chan.stID()); - debug() << "Filtered id is: " << chan.stID() << endmsg; - debug() << "Unique layer is: " << chan.stID().uniqueLayer() << endmsg; - debug() << "Unique sector is: " << chan.stID().uniqueSector() << endmsg; - - } - if (m_detType == "TT"){ - if (chan.isTT() == true && m_layerMap[m_layer]==chan.stID().uniqueLayer()) flt_ids.push_back(chan.stID()); - debug() << "Filtered id is: " << chan.stID() << endmsg; - debug() << "Unique layer is: " << chan.stID().uniqueLayer() << endmsg; - debug() << "Unique sector is: " << chan.stID().uniqueSector() << endmsg; - } - } - } - - // Loop over hits - BOOST_FOREACH( LHCb::STChannelID chan, flt_ids){ - - int sector = chan.uniqueSector(); - - double isf; - - STDAQ::chanPair aPair = m_readoutTool->offlineChanToDAQ(chan,isf); - debug() << "Board " << aPair.first <<", Channel " << aPair.second << endmsg; - unsigned int stripID = aPair.second; - unsigned int boardID = aPair.first.id(); - - if(stripID > 3071){ - warning() << "Strip ID is too high (>3071) or too low (<0)" << endmsg; - continue; - } - - // Loop over bunch crossings - std::vector::iterator prefixIter = dataPrefix.begin(); - BOOST_FOREACH( LHCb::STTELL1Datas* dataLoop,rawData){ - // Loop over TELL1 boards - LHCb::STTELL1Datas::const_iterator iterBoard = dataLoop->begin(); - for(;iterBoard != dataLoop->end(); ++iterBoard){ - int nStripsBeetle = LHCbConstants::nStripsInBeetle; - - // Check if hit is associated to TELL1 board - if ((*iterBoard)->TELL1ID() == boardID){ - // Register TELL1 board as used - boardVec.push_back(boardID); - - - std::string secName = boost::lexical_cast( sector ); - - const int stripInHybrid = stripID%(4*nStripsBeetle); - - // Filling track info - debug() << "Filling track info" << endmsg; - m_tuple_hit->column("p", trackMom); - m_tuple_hit->column("pt", pt); - m_tuple_hit->column("px", px); - m_tuple_hit->column("py", py); - m_tuple_hit->column("pz", pz); - m_tuple_hit->column("TrChi2", trChi2); - m_tuple_hit->column("TrNDoF", ndof); - m_tuple_hit->column("GhostP", ghostP); - m_tuple_hit->column("nMeas", nMeas); - m_tuple_hit->column("nLHCbIDs", nIDs); - - // Filling hit info - debug() << "Filling hit info" << endmsg; - m_tuple_hit->column("xHit", xHit); - m_tuple_hit->column("yHit", yHit); - m_tuple_hit->column("zHit", zHit); - m_tuple_hit->column("tx", dxdz); - m_tuple_hit->column("ty", dydz); - m_tuple_hit->column("pHit", pHit); - - // Filling ODIN info - debug() << "Filling ODIN info" << endmsg; - m_tuple_hit->column("odinStep", calibStepValue); - m_tuple_hit->column("runNumber",runN); - m_tuple_hit->column("evtNumber",evtN); - m_tuple_hit->column("bunchID", BCID); - m_tuple_hit->column("BCType", BCType); - - - // Filling sector Number - debug() << "Filling sector number" << endmsg; - m_tuple_hit->column("sector", sector); - - - int summedADCValue = 0; - // Loop over cluster sizes - for(int nNeigh = 0; nNeigh<=m_maxNeighbours; nNeigh++){ - std::string nStrips = boost::lexical_cast( nNeigh*2+1 ); - if (stripInHybrid > nNeigh-1 && stripInHybrid < (4*nStripsBeetle)-nNeigh){ - debug() << "Fill ADC value for " << nStrips << " strips" << endmsg; - summedADCValue= GetADCValue((*iterBoard),stripID,nStripsBeetle, - summedADCValue,nNeigh); - m_tuple_hit->column("val"+nStrips,summedADCValue); - }else{ - debug() << "Too close to boarder for next higher neighbour number!" << endmsg; - m_tuple_hit->column("val"+nStrips,-1000); - } - } - - // Write - debug() << "Write signal tuple" << endmsg; - StatusCode status = m_tuple_hit->write(); - if(status.isFailure()){ - return Error("Cannot fill signal ntuple"); - } - } - } - ++prefixIter; - } - } - } - } - - - - - // Retrieve noise - - - // Take random channel in TELL1 board - std::vector::iterator prefixIter = dataPrefix.begin(); - srand(time(NULL)); - unsigned int stripID = rand()%3072; - // Loop over bunch crossings - BOOST_FOREACH( LHCb::STTELL1Datas* dataLoop,rawData){ - // Loop over TELL1 boards - LHCb::STTELL1Datas::const_iterator iterBoard = dataLoop->begin(); - for(;iterBoard != dataLoop->end(); ++iterBoard){ - // Check if board is not registered to have hit - if (!(*iterBoard)) { - error() << "Board has not been found!" << endmsg; - continue; - } - - if(boardVec.end()==(find(boardVec.begin(),boardVec.end(),(*iterBoard)->TELL1ID()))){ - debug() << "TELL1 Board " << (*iterBoard)->TELL1ID() << " not used" << endmsg; - int nStripsBeetle = LHCbConstants::nStripsInBeetle; - - debug() << "Channel " << stripID << endmsg; - if(stripID > 3071){ - warning() << "Strip ID is too high (>3071) or too low (<0)" << endmsg; - continue; - } - - - - // Check if TELL1 board belongs to considered layer - debug() << "Check if TELL1 board belongs to considered layer" << endmsg; - debug() << stripID << endmsg; - std::vector sectorVec = m_readoutTool->sectorIDs(STTell1ID((*iterBoard)->TELL1ID())); - if (stripID/(4*nStripsBeetle)>= sectorVec.size() || m_layerMap[m_layer]!=sectorVec[stripID/(4*nStripsBeetle)].uniqueLayer() || - ((m_detType=="TT" && sectorVec[stripID/(4*nStripsBeetle)].isIT()) || - (m_detType=="IT" && sectorVec[stripID/(4*nStripsBeetle)].isTT()))){ - debug() << "Channel is not available!" << endmsg; - debug() << stripID/(4*nStripsBeetle) << ":" << sectorVec.size() << endmsg; - continue; - } - - - // Filling ODIN info - debug() << "Filling ODIN info" << endmsg; - m_tuple_noi->column("odinStep", calibStepValue); - m_tuple_noi->column("runNumber",runN); - m_tuple_noi->column("evtNumber",evtN); - m_tuple_noi->column("bunchID", BCID); - m_tuple_noi->column("BCType", BCType); - - // Filling sector Number - debug() << "Filling sector number" << endmsg; - double sector = sectorVec[stripID/(4*nStripsBeetle)].uniqueSector(); - m_tuple_noi->column("sector", sector); - - std::string secName = boost::lexical_cast( sectorVec[stripID/(4*nStripsBeetle)].uniqueSector() ); - debug() << "Sector "<< secName << " used for noise!" << endmsg; - const int stripInHybrid = stripID%(4*nStripsBeetle); - - int summedADCValue = 0; - - for(int nNeigh = 0; nNeigh<=m_maxNeighbours; nNeigh++){ - std::string nStrips = boost::lexical_cast( nNeigh*2+1 ); - if (stripInHybrid > nNeigh-1 && stripInHybrid < (4*nStripsBeetle)-nNeigh){ - debug() << "Fill ADC value for " << nStrips << " strips" << endmsg; - summedADCValue= GetADCValue((*iterBoard),stripID,nStripsBeetle,summedADCValue,nNeigh); - m_tuple_noi->column("val"+nStrips,summedADCValue); - }else{ - debug() << "Too close to boarder for next higher neighbour number!" << endmsg; - m_tuple_noi->column("val"+nStrips,-1000); - } - } - // Write - debug() << "Write noise tuple" << endmsg; - StatusCode status = m_tuple_noi->write(); - if(status.isFailure()){ - return Error("Cannot fill noise ntuple"); - } - - - } - } - ++prefixIter; - } - - - - - - - - - debug() << "Increment counters" << endmsg; - - // fill counters.... - counter("#Tracks") += tracks.size(); - for(MultiplicityMap::const_iterator it = m_multiplicityMap.begin() ; - it != m_multiplicityMap.end(); ++it) { - counter("#"+it->first) += it->second; - } - - - - - debug() << "End of execute!" << endmsg; - return StatusCode::SUCCESS; -}; - -void STADCTrackMonitor::findRefStates(const LHCb::Track& track, - const LHCb::State*& firstMeasurementState, - const LHCb::State*& lastMeasurementState) const -{ - firstMeasurementState = lastMeasurementState = 0 ; - LHCb::Track::ConstNodeRange nodes = track.nodes() ; - for( LHCb::Track::ConstNodeRange::const_iterator inode = nodes.begin() ; - inode != nodes.end() ; ++inode) - if( (*inode)->type()==LHCb::Node::HitOnTrack ) { - if( firstMeasurementState==0 || - (*inode)->z() < firstMeasurementState->z() ) - firstMeasurementState = &((*inode)->state()) ; - if( lastMeasurementState==0 || - (*inode)->z() > lastMeasurementState->z() ) - lastMeasurementState = &((*inode)->state()) ; } -} - - int STADCTrackMonitor::GetADCValue(LHCb::STTELL1Data* board, - unsigned int stripID, - int nStripsBeetle, - int summedADCValue, - int nNeigh - ) - { - if(nNeigh==0){ - summedADCValue += board->data()[(stripID)/nStripsBeetle][(stripID)%nStripsBeetle]; - }else{ - summedADCValue += board->data()[(stripID-nNeigh)/nStripsBeetle][(stripID-nNeigh)%nStripsBeetle]; - summedADCValue += board->data()[(stripID+nNeigh)/nStripsBeetle][(stripID+nNeigh)%nStripsBeetle]; - } - return summedADCValue; - } - - - -StateVector* STADCTrackMonitor::GetTTHit(const LHCb::Track* track) -{ - - const DeSTDetector::Layers& layers = m_STDet->layers(); - - int layerIndex = m_layerMap[m_layer]-9; - - if(layerIndex>3 || layerIndex<0) - { - error() << "Layer index for IT out of bound: " << layerIndex << endmsg; - return 0; - } - - double zTT = layers[layerIndex]->globalCentre().z(); - - const State& TTState = track->closestState(zTT); - - - StateVector* stateVectorTT = new StateVector(TTState.position(), TTState.slopes()); - - // loop over the sectors - const DeSTDetector::Sectors& sectorVector = m_STDet->sectors(); - DeSTDetector::Sectors::const_iterator iterS = sectorVector.begin(); - for (; iterS != sectorVector.end(); ++iterS){ - if((*iterS)->elementID().uniqueLayer()==m_layerMap[m_layer]){ - m_extrapolator->propagate((*stateVectorTT),(*iterS)->globalCentre().z()); - return stateVectorTT; - } - } - return 0; -} - -StateVector* STADCTrackMonitor::GetITHit(const LHCb::Track* track) -{ - const DeSTDetector::Layers& layers = m_STDet->layers(); - - int layerIndex = m_layerMap[m_layer]-9+(m_layerMap[m_layer]-9)/8*4; - - if(layerIndex>11 || layerIndex<0) - { - error() << "Layer index for IT out of bound: " << layerIndex << endmsg; - return 0; - } - - double zIT = layers[layerIndex]->globalCentre().z(); - - const State& ITState = track->closestState(zIT); - StateVector* stateVectorIT = new StateVector(ITState.position(), ITState.slopes()); - - // loop over the sectors - const DeSTDetector::Sectors& sectorVector = m_STDet->sectors(); - DeSTDetector::Sectors::const_iterator iterS = sectorVector.begin(); - for (; iterS != sectorVector.end(); ++iterS){ - if((*iterS)->elementID().uniqueLayer()==m_layerMap[m_layer]){ - m_extrapolator->propagate((*stateVectorIT),(*iterS)->globalCentre().z()); - return stateVectorIT; - } - } - return 0; -} - - - diff --git a/STAging/src/STClusterTrackMonitor.cpp~ b/STAging/src/STClusterTrackMonitor.cpp~ deleted file mode 100755 index 3abd7f8..0000000 --- a/STAging/src/STClusterTrackMonitor.cpp~ +++ /dev/null @@ -1,610 +0,0 @@ -// Include files -#include "STADCTrackMonitor.h" - -//event -#include "Event/Track.h" -#include "Event/TrackFitResult.h" -#include "Event/State.h" -#include "Kernel/LHCbID.h" -#include "Event/ODIN.h" -#include "Event/VeloCluster.h" -#include "Event/VeloPhiMeasurement.h" -#include "Event/VeloRMeasurement.h" -#include "Event/KalmanFitResult.h" -#include "Kernel/HitPattern.h" -#include "TrackInterfaces/IHitExpectation.h" -#include "TrackInterfaces/IVeloExpectation.h" - -// Det -#include "STDet/DeSTDetector.h" -#include "STDet/DeITDetector.h" -#include "VeloDet/DeVelo.h" -#include "Event/FitNode.h" -#include "Event/RecVertex.h" - - -// ST DAQ -#include "Kernel/STDAQDefinitions.h" -#include "Kernel/STBoardMapping.h" -#include "Kernel/LHCbConstants.h" -#include "Kernel/ISTReadoutTool.h" -#include "Kernel/STDAQDefinitions.h" - -// STTELL1Event -#include "Event/STTELL1Data.h" -#include "Kernel/STTell1Board.h" - -// gsl -#include "gsl/gsl_math.h" - -#include "GaudiKernel/AlgFactory.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/ToStream.h" - -#include -#include -#include -#include "Map.h" - -// Boost -#include -#include -#include -#include - -using namespace boost::lambda; -using namespace LHCb; -using namespace Gaudi; - - -DECLARE_ALGORITHM_FACTORY( STADCTrackMonitor ) - -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -STADCTrackMonitor::STADCTrackMonitor(const std::string& name, - ISvcLocator* pSvcLocator ) : - TrackMonitorTupleBase( name , pSvcLocator ) -{ - declareProperty("DetType",m_detType,"TT"); - declareProperty("Layer",m_layer,"TTaX"); - declareProperty("RawDataCent",m_rawDataLocCent,"/Event/Raw/TT/LCMSADCs"); - - declareProperty("ODINData",m_odinLoc,"DAQ/ODIN"); - - declareProperty("MaxNeighbours",m_maxNeighbours=3); - - - // TT layer mapping - m_layerMap.insert(std::pair("TTaX",9)); - m_layerMap.insert(std::pair("TTaU",10)); - m_layerMap.insert(std::pair("TTbV",17)); - m_layerMap.insert(std::pair("TTbX",18)); - - // IT layer mapping - m_layerMap.insert(std::pair("T1X1",9)); - m_layerMap.insert(std::pair("T1U", 10)); - m_layerMap.insert(std::pair("T1V", 11)); - m_layerMap.insert(std::pair("T1X2",12)); - - m_layerMap.insert(std::pair("T2X1",17)); - m_layerMap.insert(std::pair("T2U", 18)); - m_layerMap.insert(std::pair("T2V", 19)); - m_layerMap.insert(std::pair("T2X2",20)); - - m_layerMap.insert(std::pair("T3X1",25)); - m_layerMap.insert(std::pair("T3U", 26)); - m_layerMap.insert(std::pair("T3V", 27)); - m_layerMap.insert(std::pair("T3X2",28)); -} - -//============================================================================= -// Destructor -//============================================================================= -STADCTrackMonitor::~STADCTrackMonitor() {}; - -StatusCode STADCTrackMonitor::initialize() -{ - // Mandatory initialization of GaudiAlgorith - StatusCode sc = TrackMonitorTupleBase::initialize(); - - - - - if ( sc.isFailure() ) { return sc; } - m_veloDet = getDet( DeVeloLocation::Default ) ; - - - // Definition of readout tool - if ( m_detType=="IT" ){ - m_readoutTool = tool("ITReadoutTool"); - }else{ - m_readoutTool = tool("TTReadoutTool"); - } - - - m_veloExpectation = tool("VeloExpectation"); - m_ttExpectation = tool("TTHitExpectation"); - m_itExpectation = tool("ITHitExpectation"); - m_otExpectation = tool("OTHitExpectation"); - - m_extrapolator = tool("TrackParabolicExtrapolator"); - m_STDet = getDet(DeSTDetLocation::location(m_detType)); - - - return StatusCode::SUCCESS; -} - -//============================================================================= -// Execute -//============================================================================= -StatusCode STADCTrackMonitor::execute() -{ - - - // Create Tuple - Tuples::Tuple m_tuple_hit = GaudiTupleAlg::nTuple("HitInfo/"+m_layer, "HitInfo", CLID_ColumnWiseTuple ); - Tuples::Tuple m_tuple_noi = GaudiTupleAlg::nTuple("NoiseInfo/"+m_layer, "NoiseInfo", CLID_ColumnWiseTuple ); - - // Check if ODIN is available, otherwise skip event - bool isEverythingHere = true; - - if(!exist(m_odinLoc)){ - warning() << "ODIN data are not available!" << endmsg; - isEverythingHere = false; - } - - if(!isEverythingHere){ - warning() << "Skip event!" << endmsg; - return StatusCode::SUCCESS; - } - - // get ODIN - const LHCb::ODIN* odin = get(m_odinLoc); - int calibStepValue = odin->calibrationStep(); - debug() << "Value of ODIN calibration step is " << calibStepValue << "." << endmsg; - - // Check if track container is available,... - if(!exist(inputContainer())){ - warning() << "Tracks are not available!" << endmsg; - isEverythingHere = false; - } - - // Check if data of this bunch crossing is available,... - if(!exist(m_rawDataLocCent)){ - warning() << "Central raw data are not available!" << endmsg; - isEverythingHere = false; - } - - // ... otherwise skip event, unless next/previous bunch crossing is excluded - if(!isEverythingHere){ - warning() << "Skip event!" << endmsg; - return StatusCode::SUCCESS; - } - - - - // get the input data - LHCb::Track::Range tracks = get(inputContainer()); - - // define data prefix - std::vector dataPrefix; - - // get raw data - std::vector rawData; - - // get raw data locations - rawData.push_back(get(m_rawDataLocCent)); - dataPrefix.push_back("C"); - - - // Cast calibration step to string - std::string calibStep = boost::lexical_cast( calibStepValue ); - - //for(MultiplicityMap::iterator it = m_multiplicityMap.begin() ; - // it != m_multiplicityMap.end(); ++it) it->second = 0 ; - - - // Extract further event info - int runN = odin->runNumber(); - int evtN = odin->eventNumber(); - int BCID = odin->bunchId(); - int BCType = odin->bunchCrossingType(); - - - // Vectors to store boards and strip assigned to tracks - std::vector boardVec; - std::vector stripVec; - - - // Loop over tracks in event - BOOST_FOREACH( const LHCb::Track* track, tracks) { - std::string type = histoDirName(*track) ; - m_multiplicityMap[type] += 1; - - // Extract track momentum - double trackMom = track->p(); - - - debug() << "Track type is " << track->type() << " with momentum " << trackMom << " MeV/c." << endmsg; - - // Check if track is long (type()=3) TO DO: change to Track::Long - if(track->type()==Track::Long){ - debug() << "It's a long track!" << endmsg; - - // Extract track info - double px = track->momentum().x(); - double py = track->momentum().y(); - double pz = track->momentum().z(); - double pt = track->pt(); - double trChi2 = track->chi2(); - int ndof = track->nDoF(); - double ghostP = track->ghostProbability(); - int nMeas = track->nMeasurements(); - int nIDs = track->nLHCbIDs(); - - - - - // Get extrapolated hit position - StateVector* hit = 0; - if(m_detType=="TT"){ - hit = GetTTHit(track); - }else{ - hit = GetITHit(track); - } - - if(!hit){ - warning() << "Extrapolated hit not found! Skip event!" << endmsg; - return StatusCode::SUCCESS; - } - - // Extract hit - double xHit = hit->x(); - double yHit = hit->y(); - double zHit = hit->z(); - - - // Get hits extrapolated from track - std::vector ids(0); - if(m_detType=="TT"){ - m_ttExpectation->collect(*track,ids); - }else{ - m_itExpectation->collect(*track,ids); - } - - - std::vector< LHCb::STChannelID > flt_ids; - - // Check if hit is in right detector and layer - BOOST_FOREACH( LHCb::LHCbID chan, ids) { - if (chan.isST() == true){ - if (m_detType == "IT" ){ - if (chan.isIT() == true && m_layerMap[m_layer]==chan.stID().uniqueLayer()) flt_ids.push_back(chan.stID()); - debug() << "Filtered id is: " << chan.stID() << endmsg; - debug() << "Unique layer is: " << chan.stID().uniqueLayer() << endmsg; - debug() << "Unique sector is: " << chan.stID().uniqueSector() << endmsg; - - } - if (m_detType == "TT"){ - if (chan.isTT() == true && m_layerMap[m_layer]==chan.stID().uniqueLayer()) flt_ids.push_back(chan.stID()); - debug() << "Filtered id is: " << chan.stID() << endmsg; - debug() << "Unique layer is: " << chan.stID().uniqueLayer() << endmsg; - debug() << "Unique sector is: " << chan.stID().uniqueSector() << endmsg; - } - } - } - - // Loop over hits - BOOST_FOREACH( LHCb::STChannelID chan, flt_ids){ - - int sector = chan.uniqueSector(); - - double isf; - - STDAQ::chanPair aPair = m_readoutTool->offlineChanToDAQ(chan,isf); - debug() << "Board " << aPair.first <<", Channel " << aPair.second << endmsg; - unsigned int stripID = aPair.second; - unsigned int boardID = aPair.first.id(); - - if(stripID > 3071){ - warning() << "Strip ID is too high (>3071) or too low (<0)" << endmsg; - continue; - } - - // Loop over bunch crossings - std::vector::iterator prefixIter = dataPrefix.begin(); - BOOST_FOREACH( LHCb::STTELL1Datas* dataLoop,rawData){ - // Loop over TELL1 boards - LHCb::STTELL1Datas::const_iterator iterBoard = dataLoop->begin(); - for(;iterBoard != dataLoop->end(); ++iterBoard){ - int nStripsBeetle = LHCbConstants::nStripsInBeetle; - - // Check if hit is associated to TELL1 board - if ((*iterBoard)->TELL1ID() == boardID){ - // Register TELL1 board as used - boardVec.push_back(boardID); - - - std::string secName = boost::lexical_cast( sector ); - - const int stripInHybrid = stripID%(4*nStripsBeetle); - - // Filling track info - debug() << "Filling track info" << endmsg; - m_tuple_hit->column("p", trackMom); - m_tuple_hit->column("pt", pt); - m_tuple_hit->column("px", px); - m_tuple_hit->column("py", py); - m_tuple_hit->column("pz", pz); - m_tuple_hit->column("TrChi2", trChi2); - m_tuple_hit->column("TrNDoF", ndof); - m_tuple_hit->column("GhostP", ghostP); - m_tuple_hit->column("nMeas", nMeas); - m_tuple_hit->column("nLHCbIDs", nIDs); - - // Filling hit info - debug() << "Filling hit info" << endmsg; - m_tuple_hit->column("xHit", xHit); - m_tuple_hit->column("yHit", yHit); - m_tuple_hit->column("zHit", zHit); - - - // Filling ODIN info - debug() << "Filling ODIN info" << endmsg; - m_tuple_hit->column("odinStep", calibStepValue); - m_tuple_hit->column("runNumber",runN); - m_tuple_hit->column("evtNumber",evtN); - m_tuple_hit->column("bunchID", BCID); - m_tuple_hit->column("BCType", BCType); - - - // Filling sector Number - debug() << "Filling sector number" << endmsg; - m_tuple_hit->column("sector", sector); - - - int summedADCValue = 0; - // Loop over cluster sizes - for(int nNeigh = 0; nNeigh<=m_maxNeighbours; nNeigh++){ - std::string nStrips = boost::lexical_cast( nNeigh*2+1 ); - if (stripInHybrid > nNeigh-1 && stripInHybrid < (4*nStripsBeetle)-nNeigh){ - debug() << "Fill ADC value for " << nStrips << " strips" << endmsg; - summedADCValue= GetADCValue((*iterBoard),stripID,nStripsBeetle, - summedADCValue,nNeigh); - m_tuple_hit->column("val"+nStrips,summedADCValue); - }else{ - debug() << "Too close to boarder for next higher neighbour number!" << endmsg; - m_tuple_hit->column("val"+nStrips,-1000); - } - } - - // Write - debug() << "Write signal tuple" << endmsg; - StatusCode status = m_tuple_hit->write(); - if(status.isFailure()){ - return Error("Cannot fill signal ntuple"); - } - } - } - ++prefixIter; - } - } - } - } - - - - - // Retrieve noise - - - // Take random channel in TELL1 board - std::vector::iterator prefixIter = dataPrefix.begin(); - srand(time(NULL)); - unsigned int stripID = rand()%3072; - // Loop over bunch crossings - BOOST_FOREACH( LHCb::STTELL1Datas* dataLoop,rawData){ - // Loop over TELL1 boards - LHCb::STTELL1Datas::const_iterator iterBoard = dataLoop->begin(); - for(;iterBoard != dataLoop->end(); ++iterBoard){ - // Check if board is not registered to have hit - if (!(*iterBoard)) { - error() << "Board has not been found!" << endmsg; - continue; - } - - if(boardVec.end()==(find(boardVec.begin(),boardVec.end(),(*iterBoard)->TELL1ID()))){ - debug() << "TELL1 Board " << (*iterBoard)->TELL1ID() << " not used" << endmsg; - int nStripsBeetle = LHCbConstants::nStripsInBeetle; - - debug() << "Channel " << stripID << endmsg; - if(stripID > 3071){ - warning() << "Strip ID is too high (>3071) or too low (<0)" << endmsg; - continue; - } - - - - // Check if TELL1 board belongs to considered layer - debug() << "Check if TELL1 board belongs to considered layer" << endmsg; - debug() << stripID << endmsg; - std::vector sectorVec = m_readoutTool->sectorIDs(STTell1ID((*iterBoard)->TELL1ID())); - if (stripID/(4*nStripsBeetle)>= sectorVec.size() || m_layerMap[m_layer]!=sectorVec[stripID/(4*nStripsBeetle)].uniqueLayer() || - ((m_detType=="TT" && sectorVec[stripID/(4*nStripsBeetle)].isIT()) || - (m_detType=="IT" && sectorVec[stripID/(4*nStripsBeetle)].isTT()))){ - debug() << "Channel is not available!" << endmsg; - debug() << stripID/(4*nStripsBeetle) << ":" << sectorVec.size() << endmsg; - continue; - } - - - // Filling ODIN info - debug() << "Filling ODIN info" << endmsg; - m_tuple_noi->column("odinStep", calibStepValue); - m_tuple_noi->column("runNumber",runN); - m_tuple_noi->column("evtNumber",evtN); - m_tuple_noi->column("bunchID", BCID); - m_tuple_noi->column("BCType", BCType); - - // Filling sector Number - debug() << "Filling sector number" << endmsg; - double sector = sectorVec[stripID/(4*nStripsBeetle)].uniqueSector(); - m_tuple_noi->column("sector", sector); - - std::string secName = boost::lexical_cast( sectorVec[stripID/(4*nStripsBeetle)].uniqueSector() ); - debug() << "Sector "<< secName << " used for noise!" << endmsg; - const int stripInHybrid = stripID%(4*nStripsBeetle); - - int summedADCValue = 0; - - for(int nNeigh = 0; nNeigh<=m_maxNeighbours; nNeigh++){ - std::string nStrips = boost::lexical_cast( nNeigh*2+1 ); - if (stripInHybrid > nNeigh-1 && stripInHybrid < (4*nStripsBeetle)-nNeigh){ - debug() << "Fill ADC value for " << nStrips << " strips" << endmsg; - summedADCValue= GetADCValue((*iterBoard),stripID,nStripsBeetle,summedADCValue,nNeigh); - m_tuple_noi->column("val"+nStrips,summedADCValue); - }else{ - debug() << "Too close to boarder for next higher neighbour number!" << endmsg; - m_tuple_noi->column("val"+nStrips,-1000); - } - } - // Write - debug() << "Write noise tuple" << endmsg; - StatusCode status = m_tuple_noi->write(); - if(status.isFailure()){ - return Error("Cannot fill noise ntuple"); - } - - - } - } - ++prefixIter; - } - - - - - - - - - debug() << "Increment counters" << endmsg; - - // fill counters.... - counter("#Tracks") += tracks.size(); - for(MultiplicityMap::const_iterator it = m_multiplicityMap.begin() ; - it != m_multiplicityMap.end(); ++it) { - counter("#"+it->first) += it->second; - } - - - - - debug() << "End of execute!" << endmsg; - return StatusCode::SUCCESS; -}; - -void STADCTrackMonitor::findRefStates(const LHCb::Track& track, - const LHCb::State*& firstMeasurementState, - const LHCb::State*& lastMeasurementState) const -{ - firstMeasurementState = lastMeasurementState = 0 ; - LHCb::Track::ConstNodeRange nodes = track.nodes() ; - for( LHCb::Track::ConstNodeRange::const_iterator inode = nodes.begin() ; - inode != nodes.end() ; ++inode) - if( (*inode)->type()==LHCb::Node::HitOnTrack ) { - if( firstMeasurementState==0 || - (*inode)->z() < firstMeasurementState->z() ) - firstMeasurementState = &((*inode)->state()) ; - if( lastMeasurementState==0 || - (*inode)->z() > lastMeasurementState->z() ) - lastMeasurementState = &((*inode)->state()) ; } -} - - int STADCTrackMonitor::GetADCValue(LHCb::STTELL1Data* board, - unsigned int stripID, - int nStripsBeetle, - int summedADCValue, - int nNeigh - ) - { - if(nNeigh==0){ - summedADCValue += board->data()[(stripID)/nStripsBeetle][(stripID)%nStripsBeetle]; - }else{ - summedADCValue += board->data()[(stripID-nNeigh)/nStripsBeetle][(stripID-nNeigh)%nStripsBeetle]; - summedADCValue += board->data()[(stripID+nNeigh)/nStripsBeetle][(stripID+nNeigh)%nStripsBeetle]; - } - return summedADCValue; - } - - - -StateVector* STADCTrackMonitor::GetTTHit(const LHCb::Track* track) -{ - - const DeSTDetector::Stations& stations = m_STDet->stations(); - double zTTa = stations[0]->globalCentre().z(); - double zTTb = stations[1]->globalCentre().z(); - - const State& TTaState = track->closestState(zTTa); - const State& TTbState = track->closestState(zTTb); - StateVector* stateVectorTT; - - if(std::string::npos != m_layer.find("TTa")){ - stateVectorTT = new StateVector(TTaState.position(), TTaState.slopes()); - } - else{ - stateVectorTT = new StateVector(TTbState.position(), TTbState.slopes()); - } - - - // loop over the sectors - const DeSTDetector::Sectors& sectorVector = m_STDet->sectors(); - DeSTDetector::Sectors::const_iterator iterS = sectorVector.begin(); - for (; iterS != sectorVector.end(); ++iterS){ - if((*iterS)->elementID().uniqueLayer()==m_layerMap[m_layer]){ - m_extrapolator->propagate((*stateVectorTT),(*iterS)->globalCentre().z()); - return stateVectorTT; - } - } - return 0; -} - -StateVector* STADCTrackMonitor::GetITHit(const LHCb::Track* track) -{ - const DeSTDetector::Stations& stations = m_STDet->stations(); - double zIT1 = stations[0]->globalCentre().z(); - double zIT2 = stations[1]->globalCentre().z(); - double zIT3 = stations[2]->globalCentre().z(); - - const State& IT1State = track->closestState(zIT1); - const State& IT2State = track->closestState(zIT2); - const State& IT3State = track->closestState(zIT3); - StateVector* stateVectorIT; - - if(std::string::npos != m_layer.find("IT1")){ - stateVectorIT = new StateVector(IT1State.position(), IT1State.slopes()); - } - else if(std::string::npos != m_layer.find("IT2")){ - stateVectorIT = new StateVector(IT2State.position(), IT2State.slopes()); - } - else{ - stateVectorIT = new StateVector(IT3State.position(), IT3State.slopes()); - } - - - // loop over the sectors - const DeSTDetector::Sectors& sectorVector = m_STDet->sectors(); - DeSTDetector::Sectors::const_iterator iterS = sectorVector.begin(); - for (; iterS != sectorVector.end(); ++iterS){ - if((*iterS)->elementID().uniqueLayer()==m_layerMap[m_layer]){ - m_extrapolator->propagate((*stateVectorIT),(*iterS)->globalCentre().z()); - return stateVectorIT; - } - } - return 0; -} - - - diff --git a/STAging/src/STClusterTrackMonitor.h~ b/STAging/src/STClusterTrackMonitor.h~ deleted file mode 100755 index f319318..0000000 --- a/STAging/src/STClusterTrackMonitor.h~ +++ /dev/null @@ -1,76 +0,0 @@ -#ifndef STCLUSTERTRACKMONITOR_H -#define STCLUSTERTRACKMONITOR_H 1 - -// Include files - -#include -#include - - -// from Gaudi -#include "TrackMonitorTupleBase.h" - -// from Kernel -#include "Kernel/ISTReadoutTool.h" -#include "Kernel/STToolBase.h" - -// STTELL1Event -#include "Event/STTELL1Data.h" - -namespace LHCb { - class Track ; - class State ; -} -class DeVelo ; -class DeOTDetector ; -class IHitExpectation; -class IVeloExpectation; -class ITrackExtraplator; -class DeSTDetector; - -/** @class STClusterTrackMonitor STClusterTrackMonitor.h "TrackCheckers/TrackMonitor" - * - * Class for monitoring of cluster association to tracks - * @author Ch. Elsasser - * @date 17-6-2013 - */ - -class STClusterTrackMonitor : public TrackMonitorTupleBase { - - public: - - /** Standard construtor */ - STClusterTrackMonitor( const std::string& name, ISvcLocator* pSvcLocator ); - - /** Destructor */ - virtual ~STClusterTrackMonitor(); - - /** Algorithm initialize */ - virtual StatusCode initialize(); - - /** Algorithm execute */ - virtual StatusCode execute(); - - private: - - - - std::string m_detType; - std::string m_odinLoc; - std::string m_trackLoc; - - - - - - - typedef std::map MultiplicityMap ; - MultiplicityMap m_multiplicityMap; - - std::map m_layerMap; - - -}; - - -#endif // STCLUSTERTRACKMONITOR_H diff --git a/STAging/src/STHitIntegrator.h~ b/STAging/src/STHitIntegrator.h~ deleted file mode 100755 index 2fd77b6..0000000 --- a/STAging/src/STHitIntegrator.h~ +++ /dev/null @@ -1,87 +0,0 @@ -#ifndef STHITINTEGRATOR_H -#define STHITINTEGRATOR_H 1 - -// Include files - -#include -#include - -// from Gaudi -#include "TrackMonitorBase.h" - -// from Kernel -#include "Kernel/ISTReadoutTool.h" - -// STTELL1Event -#include "Event/STTELL1Data.h" - -#include "GaudiAlg/GaudiTupleAlg.h" - - -namespace LHCb { - class Track ; - class State ; -} -class DeVelo ; -class DeOTDetector ; -class IHitExpectation; -class IVeloExpectation; -class ITrackExtraplator; -class DeSTDetector; - -/** @class STHitIntegrator STHitIntegrator.h - * - * Class for get active area of IT/TT - * @author Ch.Elsasser - * @date 19-5-2012 - */ - -class STHitIntegrator : public GaudiTupleAlg { - - public: - - /** Standard construtor */ - STHitIntegrator( const std::string& name, ISvcLocator* pSvcLocator ); - - /** Destructor */ - virtual ~STHitIntegrator(); - - /** Algorithm initialize */ - virtual StatusCode initialize(); - - - - private: - - const DeSTDetector* m_STDet; - - std::string m_detType; - std::string m_layer; - - double m_Ntestint; - - double m_maxX; - double m_maxY; - - - - - - - - - - - - - - typedef std::map MultiplicityMap ; - MultiplicityMap m_multiplicityMap; - - std::map m_layerMap; - - -}; - - -#endif // STHITINTEGRATOR_H diff --git a/include/lhcbstyle_C.so b/include/lhcbstyle_C.so deleted file mode 100755 index a0cc1d4..0000000 --- a/include/lhcbstyle_C.so +++ /dev/null Binary files differ diff --git a/macros/CCEScan/Hamburg b/macros/CCEScan/Hamburg deleted file mode 100755 index f795a0e..0000000 --- a/macros/CCEScan/Hamburg +++ /dev/null Binary files differ diff --git a/macros/CCEScan/deplVphiN b/macros/CCEScan/deplVphiN deleted file mode 100755 index e0cbe9e..0000000 --- a/macros/CCEScan/deplVphiN +++ /dev/null Binary files differ diff --git a/macros/CCEScan/deplVphiN.C~ b/macros/CCEScan/deplVphiN.C~ deleted file mode 100755 index d5cb169..0000000 --- a/macros/CCEScan/deplVphiN.C~ +++ /dev/null @@ -1,735 +0,0 @@ -// -// deplVphiN.C -// Script to plot the result of the depletion voltage measurement as a function of the fluence -// -// Created by Christian Elsasser on 31.05.13. -// University of Zurich, elsasser@cern.ch -// Copyright only for commercial use -// - -// General include -#include -#include -#include -#include -#include -#include -#include -#include - - -// ROOT include -#include "incBasicROOT.h" -#include "incDrawROOT.h" -#include "incGraphROOT.h" -#include "incHistROOT.h" -#include "incIOROOT.h" -#include "incRooFit.h" - - -#include "basicROOTTools.h" -#include "basicRooFitTools.h" - - -#include "lhcbstyle.C" - - -#include "deplTool.h" -#include "hamburgTool.h" - -double vTop = 280.0; -double vBottom = 125.0; -double v1Bound = 160.0; -double v2Bound = 220.0; - -if (det.EqualTo("IT")) { - vTop = 250.0; - vBottom = 75.0; - v1Bound = 150.0; - v2Bound = 220.0; -} - - -// -d detector (TT/IT) -// -s set 0: all; -// 1: TT [125-160] V / IT [ 75-150] V; -// 1: TT [160-220] V / IT [150-220] V; -// 1: TT [220-280] V / IT [220-250] V; -// Voltage values are the initial depletion voltage values measured after production - - -int deplVphiN(int argc, char* argv[]); - - -int main(int argc, char* argv[]){ - - TStyle* style = lhcbStyle(); - - int argc_copy = argc; - - TApplication* theApp = new TApplication("Analysis", &argc, argv); - argc = theApp->Argc(); - argv = theApp->Argv(); - - bool runBatch = false; - - for (int i = 1; iSetBatch(); - } - } - - - - int exit = deplVphiN(argc,argv); - - - - Printf("End"); - if (!runBatch) { - theApp->Run(!runBatch); - } - - - return exit; -} - - -// Executable method -int deplVphiN(int argc, char* argv[]){ - - TVectorD v_res(5); - - gStyle->SetPadLeftMargin(0.20); - gStyle->SetTitleOffset(1.5,"Y"); - gStyle->SetTitleOffset(1.0,"X"); - - const char* diskVar = std::getenv ("DISK"); - const char* homeVar = std::getenv ("HOME"); - - - TString det("TT"); - TString lay("TTaU"); - int set = 0; - - - for (int i=1; iGet("v_fill"); - TVectorD* vp_lumi7 = (TVectorD*)f_fill->Get("v_lumi7"); - TVectorD* vp_lumi8 = (TVectorD*)f_fill->Get("v_lumi8"); - - if (!vp_fill || !vp_lumi7 || !vp_lumi8) { - Error("deplVphiN","Luminosity vectors are not available"); - return EXIT_FAILURE; - } - - TVectorD v_fill = *vp_fill; - TVectorD v_lumi7 = *vp_lumi7; - TVectorD v_lumi8 = *vp_lumi8; - - - // Get numbers of available read-out sectors and corresponding initial depletion voltage values - TVectorD* vp_sector = (TVectorD*)f_deplV->Get(Form("v_sector_%s",lay.Data())); - TVectorD* vp_vdepl = (TVectorD*)f_deplV->Get(Form("v_vdepl_%s",lay.Data())); - - if (!vp_sector || !vp_vdepl) { - Error("deplVphiN","Depletion Voltage vectors are not available"); - return EXIT_FAILURE; - } - - TVectorD v_sector = *vp_sector; - TVectorD v_vdepl = *vp_vdepl; - - Int_t effnR = nRTT; - - // No radius is taken into account in case of IT - if (det.EqualTo("IT")) { - effnR = 1; - } - - // Define the starting depletion voltage for the predictions - - double startV = 200.0; - if (set==1) { - startV = 160.0; - }else if(set==3){ - startV = 240.0; - } - if (det.EqualTo("IT")) { - startV = 100.0; - if (set==1) { - startV = 80.0; - }else if(set==3){ - startV = 120.0; - } - } - - - // Prediction calculation based on the stable damage part in the Hamburg model - int nTheoPoints = 1000; - double minTheo = 1e+6; - double maxTheo = 1e+20; - double binTheo = TMath::Power(maxTheo/minTheo,1.0/nTheoPoints); - - std::vector v_vdepl_v_theo; - std::vector v_vdepl_e_theo; - std::vector v_phiN_v_theo; - std::vector v_phiN_e_theo; - - Info("deplVphiN","Calculating theoretical prediction"); - for (int i=0; i v_func = STHamburg::getVeff(phiN_point,det.EqualTo("IT")); - Printf(" phi = %20.2f",phiN_point); - Printf(" V = (%6.2f+/-%6.2f) V",TMath::Abs(startV-v_func[0]),v_func[1]); - Printf("****"); - v_vdepl_v_theo.push_back(TMath::Abs(startV-v_func[0])); - v_vdepl_e_theo.push_back(v_func[1]); - v_phiN_v_theo.push_back(phiN_point); - v_phiN_e_theo.push_back(0.0); - } - - // Convert vectors related to predictions into root vectors - - TVectorD vr_vdepl_v_theo = getROOTVector(v_vdepl_v_theo); - TVectorD vr_vdepl_e_theo = getROOTVector(v_vdepl_e_theo); - TVectorD vr_phiN_v_theo = getROOTVector(v_phiN_v_theo); - TVectorD vr_phiN_e_theo = getROOTVector(v_phiN_e_theo); - - - - TGraphErrors* ge_theo = new TGraphErrors(vr_phiN_v_theo, - vr_vdepl_v_theo, - vr_phiN_e_theo, - vr_vdepl_e_theo); - TGraph* ge_theo_noerr = new TGraph(vr_phiN_v_theo, - vr_vdepl_v_theo); - ge_theo->SetFillColor(kGray+2); - - TMultiGraph* mg_frame = new TMultiGraph("mg_frame","Frame multigraph"); - - - - TCanvas* c_frame = new TCanvas("c_frame","Frame Canvas",800,600); - c_frame->SetLogx(); - mg_frame->Draw(""); - - // Vectors to store measurement points - - std::vector v_vdepl_v_norm_low; - std::vector v_vdepl_v_norm_mid; - std::vector v_vdepl_v_norm_high; - std::vector v_vdepl_v_r_low; - std::vector v_vdepl_v_r_mid; - std::vector v_vdepl_v_r_high; - std::vector v_vdepl_e_norm_low; - std::vector v_vdepl_e_norm_mid; - std::vector v_vdepl_e_norm_high; - std::vector v_vdepl_e_r_low; - std::vector v_vdepl_e_r_mid; - std::vector v_vdepl_e_r_high; - std::vector v_phiN_v_norm_low; - std::vector v_phiN_v_norm_mid; - std::vector v_phiN_v_norm_high; - std::vector v_phiN_v_r_low; - std::vector v_phiN_v_r_mid; - std::vector v_phiN_v_r_high; - std::vector v_phiN_e_norm_low; - std::vector v_phiN_e_norm_mid; - std::vector v_phiN_e_norm_high; - std::vector v_phiN_e_r_low; - std::vector v_phiN_e_r_mid; - std::vector v_phiN_e_r_high; - - - double vTop = 280.0; - double vBottom = 125.0; - double v1Bound = 160.0; - double v2Bound = 220.0; - - if (det.EqualTo("IT")) { - vTop = 250.0; - vBottom = 75.0; - v1Bound = 150.0; - v2Bound = 220.0; - } - - double vMax = vTop; - double vMin = vBottom; - - if (set==1) { - vMax = v1Bound; - }else if(set==2){ - vMax = v2Bound; - vMin = v1Bound; - }else if(set==3){ - vMin = v2Bound; - } - - double mean_sys = 0.0; - int nMeas = 0; - - // Loop over all fills, sectors and radii - - for (int i=0; i0.0) { - s_rad = Form("%d",(int)rTT[k]); - vn_sector_flux += "_r"+s_rad; - vn_flux7_val += "_r"+s_rad; - vn_flux8_val += "_r"+s_rad; - vn_flux7_err += "_r"+s_rad; - vn_flux8_err += "_r"+s_rad; - vn_data += "_r"+s_rad; - } - - - // Extract them from the file - TVectorD* vp_sector_flux = (TVectorD*)f_sector->Get(vn_sector_flux.Data()); - TVectorD* vp_flux7_val = (TVectorD*)f_sector->Get(vn_flux7_val.Data()); - TVectorD* vp_flux8_val = (TVectorD*)f_sector->Get(vn_flux8_val.Data()); - TVectorD* vp_flux7_err = (TVectorD*)f_sector->Get(vn_flux7_err.Data()); - TVectorD* vp_flux8_err = (TVectorD*)f_sector->Get(vn_flux8_err.Data()); - - if (!vp_sector_flux || - !vp_flux7_val || !vp_flux7_err || - !vp_flux8_val || !vp_flux8_err) { - Error("deplVphiN","Flux information for radius = %6.2f mm",rTT[k]>0.0?rTT[k]:2000.0); - continue; - } - - TVectorD* vp_data = (TVectorD*)f_data->Get(vn_data.Data()); - - if (!vp_data){ - Warning("deplVphiN","No vector found for %s",vn_data.Data()); - continue; - } - - TVectorD v_sector_flux = *vp_sector_flux; - TVectorD v_flux7_val = *vp_flux7_val; - TVectorD v_flux7_err = *vp_flux7_err; - TVectorD v_flux8_val = *vp_flux8_val; - TVectorD v_flux8_err = *vp_flux8_err; - - - TVectorD v_data = *vp_data; - - int i_flux = getClosestIndex(v_sector_flux,(double)sector); - - double flux7_val = v_flux7_val(i_flux); - double flux7_err = v_flux8_err(i_flux); - double flux8_val = v_flux7_val(i_flux); - double flux8_err = v_flux8_err(i_flux); - - double lumi7 = v_lumi7(i); - double lumi8 = v_lumi8(i); - - double vdepl_prod = v_vdepl(j); - - double vdepl_val = v_data(0); - double vdepl_err = v_data(1); - double vdepl_sys = v_data(2); - - double vdepl_tot = vdepl_err; - - double phiN_val = STTool::FLUKAConvFac; - double phiN_err = STTool::FLUKAConvFac; - - phiN_val *= (flux7_val*lumi7+flux8_val*lumi8); - phiN_err *= (flux7_err*lumi7+flux8_err*lumi8); - - - // Define maximal and minimal Vdepl values to be drawn - double plotMin = 40.0; - double plotMax = 290.0; - - if (det.EqualTo("TT")) { - plotMin = 80.0; - } - - - // Remove points with too high uncertainty - if (vdepl_err>20.0 ){ - Warning("deplVphiN","Too high or low uncertainty on V_depl"); - continue; - } - - - - if (vdepl_val>plotMax || vdepl_valvMax || vdepl_prod<=vMin) { - continue; - } - - // Add systematic uncertainty to get mean value - mean_sys += vdepl_sys; - nMeas++; - - // Remove excluded sectors - if (STTool::IsExcluded(sector)) { - Warning("deplVphiN","Sector excluded"); - continue; - } - - Info("deplVphiN","Result for %s/%s sector %d, fill %d, radius %7.2f mm:", - det.Data(),lay.Data(),sector,fill,rTT[k]>0.0?rTT[k]:2000.0); - Printf(" phi_N = (%10.2f+/-%7.2f)*10^12 cm^{-2}", - phiN_val/1e12,phiN_err/1e12); - Printf(" V_depl = (%6.2f+/-%6.2f(stat)+/-%6.2f(syst)) V", - vdepl_val,vdepl_err,vdepl_sys); - Printf(" V_depl,prod = %6.2f V", - vdepl_prod); - - // Define data point in plot - TMarker* m_res = new TMarker(phiN_val,vdepl_val,20); - TLine* l_hor = new TLine(phiN_val-phiN_err,vdepl_val, - phiN_val+phiN_err,vdepl_val); - TLine* l_ver = new TLine(phiN_val,vdepl_val-vdepl_err, - phiN_val,vdepl_val+vdepl_err); - - // Add them to the proper vector - - if (det.EqualTo("TT")) { - if (vdepl_prod>v2Bound) { - if (rTT[k]>0.0) { - v_vdepl_v_r_high.push_back(vdepl_val); - v_vdepl_e_r_high.push_back(vdepl_tot); - v_phiN_v_r_high.push_back(phiN_val); - v_phiN_e_r_high.push_back(phiN_err); - }else{ - v_vdepl_v_norm_high.push_back(vdepl_val); - v_vdepl_e_norm_high.push_back(vdepl_tot); - v_phiN_v_norm_high.push_back(phiN_val); - v_phiN_e_norm_high.push_back(phiN_err); - } - }else if(vdepl_prod>v1Bound) { - if (rTT[k]>0.0) { - v_vdepl_v_r_mid.push_back(vdepl_val); - v_vdepl_e_r_mid.push_back(vdepl_tot); - v_phiN_v_r_mid.push_back(phiN_val); - v_phiN_e_r_mid.push_back(phiN_err); - }else{ - v_vdepl_v_norm_mid.push_back(vdepl_val); - v_vdepl_e_norm_mid.push_back(vdepl_tot); - v_phiN_v_norm_mid.push_back(phiN_val); - v_phiN_e_norm_mid.push_back(phiN_err); - } - }else{ - if (rTT[k]>0.0) { - v_vdepl_v_r_low.push_back(vdepl_val); - v_vdepl_e_r_low.push_back(vdepl_tot); - v_phiN_v_r_low.push_back(phiN_val); - v_phiN_e_r_low.push_back(phiN_err); - }else{ - v_vdepl_v_norm_low.push_back(vdepl_val); - v_vdepl_e_norm_low.push_back(vdepl_tot); - v_phiN_v_norm_low.push_back(phiN_val); - v_phiN_e_norm_low.push_back(phiN_err); - } - } - }else{ - if (vdepl_prod>v2Bound) { - v_vdepl_v_norm_high.push_back(vdepl_val); - v_vdepl_e_norm_high.push_back(vdepl_tot); - v_phiN_v_norm_high.push_back(phiN_val); - v_phiN_e_norm_high.push_back(phiN_err); - }else if(vdepl_prod>v1Bound) { - v_vdepl_v_norm_mid.push_back(vdepl_val); - v_vdepl_e_norm_mid.push_back(vdepl_tot); - v_phiN_v_norm_mid.push_back(phiN_val); - v_phiN_e_norm_mid.push_back(phiN_err); - }else{ - v_vdepl_v_norm_low.push_back(vdepl_val); - v_vdepl_e_norm_low.push_back(vdepl_tot); - v_phiN_v_norm_low.push_back(phiN_val); - v_phiN_e_norm_low.push_back(phiN_err); - } - } - } - } - } - - // Calculate mean systematic uncertainty - if (nMeas>0) { - mean_sys /= nMeas; - } - - // Add systmatic uncertainty from the extraction procedure of the depletion voltage - if (det.EqualTo("IT")) { - mean_sys = TMath::Sqrt(mean_sys*mean_sys+STTool::voltSyst_IT*STTool::voltSyst_IT); - } - if (det.EqualTo("TT")) { - mean_sys = TMath::Sqrt(mean_sys*mean_sys+STTool::voltSyst_TT*STTool::voltSyst_TT); - } - - - vr_vdepl_v_theo+=mean_sys; - - TGraph* ge_theo_high = new TGraph(vr_phiN_v_theo, - vr_vdepl_v_theo); - ge_theo_high->SetLineStyle(9); - - vr_vdepl_v_theo-=2*mean_sys; - - TGraph* ge_theo_low = new TGraph(vr_phiN_v_theo, - vr_vdepl_v_theo); - ge_theo_low->SetLineStyle(9); - - - // Transfrom STL vectors to ROOT vectors - TVectorD vr_vdepl_v_norm_low = getROOTVector(v_vdepl_v_norm_low); - TVectorD vr_vdepl_v_norm_mid = getROOTVector(v_vdepl_v_norm_mid); - TVectorD vr_vdepl_v_norm_high = getROOTVector(v_vdepl_v_norm_high); - TVectorD vr_vdepl_v_r_low = getROOTVector(v_vdepl_v_r_low); - TVectorD vr_vdepl_v_r_mid = getROOTVector(v_vdepl_v_r_mid); - TVectorD vr_vdepl_v_r_high = getROOTVector(v_vdepl_v_r_high); - TVectorD vr_vdepl_e_norm_low = getROOTVector(v_vdepl_e_norm_low); - TVectorD vr_vdepl_e_norm_mid = getROOTVector(v_vdepl_e_norm_mid); - TVectorD vr_vdepl_e_norm_high = getROOTVector(v_vdepl_e_norm_high); - TVectorD vr_vdepl_e_r_low = getROOTVector(v_vdepl_e_r_low); - TVectorD vr_vdepl_e_r_mid = getROOTVector(v_vdepl_e_r_mid); - TVectorD vr_vdepl_e_r_high = getROOTVector(v_vdepl_e_r_high); - TVectorD vr_phiN_v_norm_low = getROOTVector(v_phiN_v_norm_low); - TVectorD vr_phiN_v_norm_mid = getROOTVector(v_phiN_v_norm_mid); - TVectorD vr_phiN_v_norm_high = getROOTVector(v_phiN_v_norm_high); - TVectorD vr_phiN_v_r_low = getROOTVector(v_phiN_v_r_low); - TVectorD vr_phiN_v_r_mid = getROOTVector(v_phiN_v_r_mid); - TVectorD vr_phiN_v_r_high = getROOTVector(v_phiN_v_r_high); - TVectorD vr_phiN_e_norm_low = getROOTVector(v_phiN_e_norm_low); - TVectorD vr_phiN_e_norm_mid = getROOTVector(v_phiN_e_norm_mid); - TVectorD vr_phiN_e_norm_high = getROOTVector(v_phiN_e_norm_high); - TVectorD vr_phiN_e_r_low = getROOTVector(v_phiN_e_r_low); - TVectorD vr_phiN_e_r_mid = getROOTVector(v_phiN_e_r_mid); - TVectorD vr_phiN_e_r_high = getROOTVector(v_phiN_e_r_high); - - - - // Generate graphs - double markerSize = 1.0; - - TGraphErrors* ge_norm_low = new TGraphErrors(vr_phiN_v_norm_low, - vr_vdepl_v_norm_low, - vr_phiN_e_norm_low, - vr_vdepl_e_norm_low); - ge_norm_low->SetMarkerSize(markerSize); - ge_norm_low->SetMarkerStyle(23); - ge_norm_low->SetMarkerColor(kGreen+3); - ge_norm_low->SetLineColor(kGreen+3); - - TGraphErrors* ge_norm_mid = new TGraphErrors(vr_phiN_v_norm_mid, - vr_vdepl_v_norm_mid, - vr_phiN_e_norm_mid, - vr_vdepl_e_norm_mid); - ge_norm_mid->SetMarkerSize(markerSize); - ge_norm_mid->SetMarkerStyle(22); - ge_norm_mid->SetMarkerColor(kBlue); - ge_norm_mid->SetLineColor(kBlue); - - TGraphErrors* ge_norm_high = new TGraphErrors(vr_phiN_v_norm_high, - vr_vdepl_v_norm_high, - vr_phiN_e_norm_high, - vr_vdepl_e_norm_high); - ge_norm_high->SetMarkerSize(markerSize); - ge_norm_high->SetMarkerStyle(20); - ge_norm_high->SetMarkerColor(kRed); - ge_norm_high->SetLineColor(kRed); - - - TGraphErrors* ge_r_low = new TGraphErrors(vr_phiN_v_r_low, - vr_vdepl_v_r_low, - vr_phiN_e_r_low, - vr_vdepl_e_r_low); - ge_r_low->SetMarkerSize(markerSize); - ge_r_low->SetMarkerStyle(32); - ge_r_low->SetMarkerColor(kGreen+3); - ge_r_low->SetLineColor(kGreen+3); - - TGraphErrors* ge_r_mid = new TGraphErrors(vr_phiN_v_r_mid, - vr_vdepl_v_r_mid, - vr_phiN_e_r_mid, - vr_vdepl_e_r_mid); - ge_r_mid->SetMarkerSize(markerSize); - ge_r_mid->SetMarkerStyle(26); - ge_r_mid->SetMarkerColor(kBlue); - ge_r_mid->SetLineColor(kBlue); - - TGraphErrors* ge_r_high = new TGraphErrors(vr_phiN_v_r_high, - vr_vdepl_v_r_high, - vr_phiN_e_r_high, - vr_vdepl_e_r_high); - ge_r_high->SetMarkerSize(markerSize); - ge_r_high->SetMarkerStyle(24); - ge_r_high->SetMarkerColor(kRed); - ge_r_high->SetLineColor(kRed); - - mg_frame->Add(ge_norm_low,"pe"); - mg_frame->Add(ge_norm_mid,"pe"); - mg_frame->Add(ge_norm_high,"pe"); - mg_frame->Add(ge_r_low,"p"); - mg_frame->Add(ge_r_mid,"p"); - mg_frame->Add(ge_r_high,"p"); - mg_frame->Draw("a"); - - // Plot graphs - - mg_frame->GetXaxis()->SetTitle("#phi_{1 MeV-n,eq} [cm^{-2}]"); - mg_frame->GetYaxis()->SetTitle("#it{V}_{depl} [V]"); - mg_frame->GetXaxis()->SetLimits(1e+8,1e+14); - mg_frame->GetXaxis()->SetLabelOffset(-0.005); - mg_frame->GetXaxis()->SetTitleOffset(1.2); - mg_frame->SetMaximum(300.0); - mg_frame->SetMinimum(0.0); - - - - // Draw legend - - double baseY = 0.30; - if (det.EqualTo("IT")) { - baseY = 0.90; - } - - - - TMarker* m_leg = new TMarker(0.25,baseY,20); - TLatex* t_leg = new TLatex(0.27,baseY,"#it{V}_{depl} > 225 V"); - t_leg->SetTextSize(0.04); - t_leg->SetTextAlign(12); - t_leg->SetNDC(); - m_leg->SetNDC(); - - if (det.EqualTo("TT") && (set==3 || set==0)) { - m_leg->SetX(0.25); - m_leg->SetY(baseY); - m_leg->SetMarkerColor(kRed); - m_leg->SetMarkerStyle(20); - m_leg->DrawClone(); - m_leg->SetMarkerStyle(24); - m_leg->SetX(0.55); - if (det.EqualTo("TT")) { - m_leg->DrawClone(); - } - t_leg->DrawLatex(0.30,baseY,Form("#it{V}_{depl} #in [%3.0f,%3.0f] V",v2Bound,vTop)); - if (det.EqualTo("TT")) { - t_leg->DrawLatex(0.60,baseY,"Innermost region"); - } - baseY -= 0.05; - } - - if (det.EqualTo("TT") && (set==2 || set==0)) { - m_leg->SetX(0.25); - m_leg->SetY(baseY); - m_leg->SetMarkerColor(kBlue); - m_leg->SetMarkerStyle(22); - m_leg->DrawClone(); - m_leg->SetMarkerStyle(26); - m_leg->SetX(0.55); - if (det.EqualTo("TT")) { - m_leg->DrawClone(); - } - t_leg->DrawLatex(0.30,baseY,Form("#it{V}_{depl} #in [%3.0f,%3.0f] V",v1Bound,v2Bound)); - if (det.EqualTo("TT")) { - t_leg->DrawLatex(0.60,baseY,"Innermost region"); - } - baseY -= 0.05; - } - - if (det.EqualTo("IT") && (set==1 || set==0)) { - m_leg->SetX(0.25); - m_leg->SetY(baseY); - m_leg->SetMarkerColor(kGreen+3); - m_leg->SetMarkerStyle(23); - m_leg->DrawClone(); - m_leg->SetMarkerStyle(32); - m_leg->SetX(0.55); - if (det.EqualTo("TT")) { - m_leg->DrawClone(); - } - t_leg->DrawLatex(0.30,baseY,Form("#it{V}_{depl} #in [%3.0f,%3.0f] V",vBottom,v1Bound)); - if (det.EqualTo("TT")) { - t_leg->DrawLatex(0.60,baseY,"Innermost region"); - } - baseY -= 0.05; - } - - - - TLine* l_line = new TLine(0.20,baseY,0.25,baseY); - l_line->DrawLineNDC(0.23,baseY,0.27,baseY); - t_leg->DrawLatex(0.30,baseY,Form("Hamburg Model (#it{V}_{depl} = %3.0f V)",startV)); - - mg_frame->GetHistogram()->Draw("axissame"); - - - - Printf("=========================================\n"); - - return EXIT_SUCCESS; -} - - diff --git a/macros/CCEScan/deplVphiN_normalized b/macros/CCEScan/deplVphiN_normalized deleted file mode 100755 index af2da75..0000000 --- a/macros/CCEScan/deplVphiN_normalized +++ /dev/null Binary files differ diff --git a/macros/CCEScan/landauFit b/macros/CCEScan/landauFit deleted file mode 100755 index af96281..0000000 --- a/macros/CCEScan/landauFit +++ /dev/null Binary files differ diff --git a/macros/CCEScan/landauFit.C~ b/macros/CCEScan/landauFit.C~ deleted file mode 100755 index aacb851..0000000 --- a/macros/CCEScan/landauFit.C~ +++ /dev/null @@ -1,521 +0,0 @@ -// -// landauFit.C -// Obtaining the MPV from a fit of the ADC counts distribution -// -// Created by Christian Elsasser on 31.05.13. -// University of Zurich, elsasser@cern.ch -// Copyright only for commercial use -// - -// General include -#include -#include -#include -#include -#include -#include -#include -#include - - -// ROOT include -#include "incBasicROOT.h" -#include "incDrawROOT.h" -#include "incHistROOT.h" -#include "incIOROOT.h" -#include "incRooFit.h" - - -#include "basicROOTTools.h" -#include "basicRooFitTools.h" - - -#include "lhcbstyle.C" - - -#include "deplTool.h" - - -// -d detector (TT/IT) -// -l layer (TTaX,TTaU,TTbV,TTbX, T[1-3]{X1,U,V,X2}) -// -s read-out sector number -// -f fill number -// -c Calibration step (Odin step) [0-65] -// -r consider only hit in a distance closer than r mm from the beam axis -// -v number of summed up strips -// -ns Do not save the value in the root file -// -ss Perform a temporary save - - -int landauFit(int argc, char* argv[]); - - -int main(int argc, char* argv[]){ - - TStyle* style = lhcbStyle(); - - int argc_copy = argc; - - TApplication* theApp = new TApplication("Analysis", &argc, argv); - argc = theApp->Argc(); - argv = theApp->Argv(); - - bool runBatch = false; - - for (int i = 1; iSetBatch(); - } - } - - - - int exit = landauFit(argc,argv); - - - - Printf("End"); - if (!runBatch) { - theApp->Run(!runBatch); - } - - - return exit; -} - - -// Executable method -int landauFit(int argc, char* argv[]){ - - TVectorD v_res(5); - - gStyle->SetPadLeftMargin(0.20); - gStyle->SetTitleOffset(1.5,"Y"); - gStyle->SetTitleOffset(1.0,"X"); - - const char* diskVar = std::getenv ("DISK"); - - int fill = 2083; - TString det("TT"); - TString lay("TTaU"); - int sector = 2634; - int caliStep = 0; - int val = 3; - bool notSave = false; - bool temSave = false; - double radius = -10.0; - - for (int i=1; iGet(tn_sig.Data()); - if (!t_sig) { - Error("landauFit","Signal tree not found!"); - Printf("=========================================\n"); - return EXIT_FAILURE; - } - - TTree* t_bkg = (TTree*)f_data->Get(tn_bkg.Data()); - if (!t_bkg) { - Error("landauFit","Background tree not found!"); - Printf("=========================================\n"); - return EXIT_FAILURE; - } - - TH1F* h_sig = new TH1F("h_sig","Signal histo", - 208,-32.0,176.0); - TH1F* h_bkg = new TH1F("h_bkg","Background histo", - 208,-32.0,176.0); - - // Select calibration and read-out sector - TCut cu_bkg(Form("(odinStep==%d) && (sector==%d)", - caliStep,sector)); - - // Select only first run in fill if you consider calibration step 0 - if (caliStep==0) { - cu_bkg = cu_bkg && TCut(Form("(runNumber==%d)",STTool::zeroRunMap[fill])); - } - - // Select tracks based on their chi^2 and ghost prob - TCut cu_sel = STTool::cu_track_TT; - if (strcmp(det.Data(),"TT")==0) { - cu_sel = STTool::cu_track_IT; - } - - // Get full selection for signal tracks/hits including radius - TCut cu_sig = cu_bkg && cu_sel && TCut(Form("(TMath::Sqrt(xHit*xHit+yHit*yHit)<%f)",radius)); - - // Fill Histograms - t_sig->Draw(Form("val%d>>h_sig",val),cu_sig,"goff"); - t_bkg->Draw(Form("val%d>>h_bkg",val),cu_bkg,"goff"); - - - delete t_sig; - delete t_bkg; - - // Skip histogram if there are too few entries in the histograms - if (h_sig->GetEntries()<10 || h_bkg->GetEntries()<10) { - Error("landauFit","Too few entries in histograms. Skip this step!"); - return EXIT_FAILURE; - } - - // ADC value - RooRealVar* adc = new RooRealVar("adc","Signal height", - h_sig->GetXaxis()->GetXmin(), - h_sig->GetXaxis()->GetXmax(), - "ADC value"); - - // RooFit histograms - RooDataHist* hr_bkg = new RooDataHist("hr_bkg","Histogram of raw noise", - *adc,RooFit::Import(*h_bkg)); - RooDataHist* hr_sig = new RooDataHist("hr_sig","Histogram of raw signal", - *adc,RooFit::Import(*h_sig)); - - // PDF for noise - RooRealVar* nMu = new RooRealVar("nMu","Mean of raw noise distribution",0.0,-10.0,10.0,""); - RooRealVar* nSig1 = new RooRealVar("nSig1","Sigma 1 of raw noise distribution",3.0,0.5,6.0,""); - RooRealVar* nSig2 = new RooRealVar("nSig2","Sigma 2 of raw noise distribution",5.0,0.5,25.0,""); - RooRealVar* nFrac = new RooRealVar("nFrac","Fractition of Gaussian 1 to total noise distribution",0.5,0.0,1.0); - - RooGaussian* nG1 = new RooGaussian("nG1","Gaussian 1 of raw noise distribution",*adc,*nMu,*nSig1); - RooGaussian* nG2 = new RooGaussian("nG2","Gaussian 2 of raw noise distribution",*adc,*nMu,*nSig2); - - RooAddPdf* nPDF = new RooAddPdf("nPDF","Total distribution of raw noise", - RooArgList(*nG1,*nG2), - RooArgList(*nFrac)); - - // Fit and plot noise sample - RooFitResult* nRes = nPDF->fitTo(*hr_bkg,RooFit::Save()); - RooPlot* nFrame = adc->frame(RooFit::Title(Form("Noise distribution for %4.0f V and %4.1f ns",volt,time))); - hr_bkg->plotOn(nFrame,RooFit::DataError(RooAbsData::SumW2),RooFit::DataError(RooAbsData::SumW2)); - nPDF->plotOn(nFrame,RooFit::LineColor(kBlue)); - - - // average noise value - v_res(2) = nSig1->getVal()*nFrac->getVal()+nSig2->getVal()*(1.0-nFrac->getVal()); - // uncertainty on average noise - v_res(3) = TMath::Sqrt(TMath::Power(nSig1->getError()*nFrac->getVal(),2.0)+ - TMath::Power(nSig2->getError()*(1.0-nFrac->getVal()),2.0)+ - TMath::Power((nSig1->getVal()-nSig2->getVal())*nFrac->getError(),2.0)); - - - // PDF for signal - - // Constrain parameters fixed in the noise fit - nMu->setConstant(kTRUE); - nSig1->setConstant(kTRUE); - nSig2->setConstant(kTRUE); - nFrac->setConstant(kTRUE); - - RooGaussian* sG1 = new RooGaussian("sG1","Gaussian 1",*adc,*nMu,*nSig1); - RooGaussian* sG2 = new RooGaussian("sG2","Gaussian 2",*adc,*nMu,*nSig2); - - RooAddPdf* sG = new RooAddPdf("sG","Gaussian smearing", - RooArgList(*sG1,*sG2), - RooArgList(*nFrac)); - - RooAddPdf* nG = new RooAddPdf("nG","Distribution of raw noise", - RooArgList(*nG1,*nG2), - RooArgList(*nFrac)); - - - RooRealVar* mpv = new RooRealVar("mpv","Most probable ADC value", - TMath::Max(10.0,h_sig->GetBinCenter(h_sig->GetMaximumBin())) - ,1.0,60.0,""); - RooRealVar* spread = new RooRealVar("spread","Sigma of Landau distribution" - ,3.0,0.1,10.0,""); - - RooRealVar* ampv = new RooRealVar("ampv","Most probable ADC value", - TMath::Max(10.0,h_sig->GetBinCenter(h_sig->GetMaximumBin())) - ,1.0,60.0,""); - RooRealVar* aspread = new RooRealVar("aspread","Sigma of Landau distribution" - ,3.0,0.1,10.0,""); - - // Variables for the Landau coming from photo conversion - RooFormulaVar* scaledMpv = new RooFormulaVar("scaledMpv","2*mpv",RooArgList(*mpv)); - RooFormulaVar* scaledSpread = new RooFormulaVar("scaledSpread","2*spread",RooArgList(*spread)); - - - RooLandau* landauDist = new RooLandau("landauDist","Landau distribution of ADC values",*adc,*mpv,*spread); - RooLandau* landauDistScaled = new RooLandau("landauDistScaled","Landau distribution of conversion",*adc,*scaledMpv,*scaledSpread); - - RooLandau* altlandauDist = new RooLandau("altlandauDist","AlternativeLandau distribution of ADC values",*adc,*ampv,*aspread); - - RooRealVar* cFrac = new RooRealVar("cFrac","Fraction of photo conversion",0.025); - if (det.EqualTo("IT")) { - cFrac->setVal(0.0); - } - - RooRealVar* sFrac = new RooRealVar("sFrac","Fraction of signal hits",0.8,0.0,1.0); - RooRealVar* aFrac = new RooRealVar("aFrac","Fraction of signal hits",0.8,0.0,1.0); - - adc->setBins(10000,"cache") ; - - RooFFTConvPdf* landauConv = new RooFFTConvPdf("landauConv","Landau (x) gauss",*adc,*landauDist,*sG); - RooFFTConvPdf* landauConvPhoto = new RooFFTConvPdf("landauConvPhoto","Landau Conv (x) gauss",*adc,*landauDistScaled,*sG); - - RooAddPdf* sTot = new RooAddPdf("sTot","Total signal distribution",*landauConvPhoto,*landauConv,*cFrac); - RooAddPdf* sPDF = new RooAddPdf("sPDF","Total distribution",*sTot,*nG,*sFrac); - - // Alternative distribution only done from a single Landau - RooFFTConvPdf* altlandauConv = new RooFFTConvPdf("altlandauConv","Landau Conv (x) gauss",*adc,*altlandauDist,*sG); - RooAddPdf* aPDF = new RooAddPdf("aPDF","Alternative total distribution",*altlandauConv,*nG,*aFrac); - - - // Fit the distribution - RooFitResult* sRes = sPDF->fitTo(*hr_sig,RooFit::Save()); - - - - Info("landauFit","Basline result:"); - Printf(" mpv = %8.4f+/-%8.4f",mpv->getVal(),mpv->getError()); - Printf(" spread = %8.4f+/-%8.4f",spread->getVal(),spread->getError()); - Printf(" conv Frac = %8.4f+/-%8.4f",cFrac->getVal(),cFrac->getError()); - Printf(" noise Frac = %8.4f+/-%8.4f",1.0 - sFrac->getVal(),sFrac->getError()); - - // mpv value: 1st place in result vector - v_res(0) = mpv->getVal(); - // mpv uncertainty: 2nd place in result vector - v_res(1) = mpv->getError(); - - - - - RooPlot* sFrame = adc->frame(RooFit::Title(Form("Signal distribution for %4.0f V and %4.1f ns",volt,time))); - hr_sig->plotOn(sFrame,RooFit::DataError(RooAbsData::SumW2),RooFit::DataError(RooAbsData::SumW2)); - - sPDF->plotOn(sFrame,RooFit::Components(*landauConv),RooFit::LineColor(kRed),RooFit::LineStyle(9)); - - Info("landauFit","Results: "); - Printf(" MPV: %5.2f+/-%5.2f ADC counts",v_res(0),v_res(1)); - Printf(" Noise: %5.2f+/-%5.2f ADC counts",v_res(2),v_res(3)); - - if (nRes->status()==0) { - Info("landauFit","Signal fit status: %d",nRes->status()); - }else{ - Warning("landauFit","Signal fit status: %d",nRes->status()); - } - - if (sRes->status()==0) { - Info("landauFit","Noise fit status: %d",sRes->status()); - }else{ - Warning("landauFit","Noise fit status: %d",sRes->status()); - } - - - // Plot the results - - TCanvas* c_sig = new TCanvas("c_sig","Signal Canvas",800,600); - c_sig->SetLeftMargin(0.20); - SwapParenthesis(sFrame); - ReplaceSubstring(sFrame,"Events","Hits"); - ReplaceSubstring(sFrame,"1 ",""); - sFrame->GetYaxis()->SetTitleOffset(1.5); - sFrame->Draw(); - - TCanvas* c_bkg = new TCanvas("c_bkg","Background Canvas",800,600); - c_bkg->SetLeftMargin(0.20); - nFrame->GetXaxis()->SetLimits(-32.0,32.0); - SwapParenthesis(nFrame); - ReplaceSubstring(nFrame,"Events","Hits"); - ReplaceSubstring(nFrame,"1 ",""); - nFrame->GetYaxis()->SetTitleOffset(1.5); - nFrame->Draw(); - - - // Save data and plots - - TFile* f_output; - - while (!gSystem->Exec(Form("ls %s &> /dev/null",fn_dummy.Data()))) { - Info("landauFit","Waiting until file is closed"); - gSystem->Sleep(1000); - } - gSystem->Exec(Form("touch %s &> /dev/null",fn_dummy.Data())); - - f_output = new TFile(fn_output.Data(),"UPDATE"); - - - TString dn_result = TString::Format("%s/%s/%d/%d", - det.Data(), - lay.Data(), - sector, - fill); - - // Get right directory - bool isThere = f_output->GetListOfKeys()->Contains(det.Data()); - if (!isThere) { - f_output->mkdir(det.Data()); - } - - isThere = f_output->GetListOfKeys()->Contains(Form("%s/%s", - det.Data(), - lay.Data())); - if (!isThere) { - f_output->mkdir(Form("%s/%s", - det.Data(), - lay.Data())); - } - - isThere = f_output->GetListOfKeys()->Contains(Form("%s/%s/%d", - det.Data(), - lay.Data(), - sector)); - if (!isThere) { - f_output->mkdir(Form("%s/%s/%d", - det.Data(), - lay.Data(), - sector)); - } - - isThere = f_output->GetListOfKeys()->Contains(Form("%s/%s/%d/%d", - det.Data(), - lay.Data(), - sector, - fill)); - if (!isThere) { - f_output->mkdir(Form("%s/%s/%d/%d", - det.Data(), - lay.Data(), - sector, - fill)); - } - - f_output->cd(Form("%s/%s/%d/%d", - det.Data(), - lay.Data(), - sector, - fill)); - if(!notSave){ - gDirectory->Delete(Form("%s;*", - vn_result.Data())); - - v_res.Write(Form("%s", - vn_result.Data())); - - Info("landauFit","Result written to %s",fn_output.Data()); - } - - f_output->Close(); - gSystem->Exec(Form("rm -f %s &> /dev/null",fn_dummy.Data())); - - if (!notSave) { - // Save picture - gSystem->Exec(Form("mkdir -p %s",dn_graph.Data())); - gSystem->Exec(Form("rm %s &> /dev/null",fn_graph_s.Data())); - gSystem->Exec(Form("rm %s &> /dev/null",fn_graph_n.Data())); - c_sig->SaveAs(fn_graph_s.Data()); - c_bkg->SaveAs(fn_graph_n.Data()); - } - - - Printf("=========================================\n"); - - return EXIT_SUCCESS; -} - - diff --git a/macros/CCEScan/plotSector b/macros/CCEScan/plotSector deleted file mode 100755 index 8edee82..0000000 --- a/macros/CCEScan/plotSector +++ /dev/null Binary files differ diff --git a/macros/CCEScan/plotSector.C b/macros/CCEScan/plotSector.C index 0105a19..5242271 100755 --- a/macros/CCEScan/plotSector.C +++ b/macros/CCEScan/plotSector.C @@ -336,7 +336,10 @@ double feb2016 = 1455926400.0; double may2016 = 1463040975.0; double sep2016 = 1475243957.0; - mg_frame->GetXaxis()->SetLimits(jan2010, sep2016); + double nov2016 = 1480427833.0; + double feb2018 = 1519776000.0; + double jan2019 = 1546300800.0; + mg_frame->GetXaxis()->SetLimits(jan2010, feb2018); mg_frame->SetMinimum(100.0); mg_frame->SetMaximum(300.0); diff --git a/macros/CCEScan/plotSector_gcc3.4.14 b/macros/CCEScan/plotSector_gcc3.4.14 deleted file mode 100755 index 5149ab5..0000000 --- a/macros/CCEScan/plotSector_gcc3.4.14 +++ /dev/null Binary files differ diff --git a/macros/CCEScan/pulseFit b/macros/CCEScan/pulseFit deleted file mode 100755 index 7b72b68..0000000 --- a/macros/CCEScan/pulseFit +++ /dev/null Binary files differ diff --git a/macros/CCEScan/ratioFit b/macros/CCEScan/ratioFit deleted file mode 100755 index 62660e9..0000000 --- a/macros/CCEScan/ratioFit +++ /dev/null Binary files differ diff --git a/macros/CCEScan/ratioFit_spline b/macros/CCEScan/ratioFit_spline deleted file mode 100755 index ba01189..0000000 --- a/macros/CCEScan/ratioFit_spline +++ /dev/null Binary files differ diff --git a/macros/CCEScan/runAllITForGivenFill.sh~ b/macros/CCEScan/runAllITForGivenFill.sh~ deleted file mode 100755 index 4be5493..0000000 --- a/macros/CCEScan/runAllITForGivenFill.sh~ +++ /dev/null @@ -1,28 +0,0 @@ -#!/bin/bash -# Script to submit all fits to evaluate CCEscans -# run from spinor with: . ./runAllITForGivenFill.sh 3960 -fill=$1 -echo "****************************************" -echo "Running all fills and sectors in IT" -echo "****************************************" -echo "Full sectors" -echo "Analyzing fill "$fill -echo "****************************************" -for sector in `cat ITsectors.dat` -do -. ./runFullFitLoop.sh IT T3X2 $sector $fill -#./runFullFitLoop.sh $sector $fill -done -echo "****************************************" -echo "Inner most" -for sector in `cat ITinner.dat` -do -. ./runFullFitLoop.sh IT T3X2 $sector $fill 175 -#./runFullFitLoop.sh $sector $fill 175 -done - - - -echo "Done!" - -echo "****************************************" \ No newline at end of file diff --git a/macros/CCEScan/runAllTTForGivenFill.sh~ b/macros/CCEScan/runAllTTForGivenFill.sh~ deleted file mode 100755 index 2a1df58..0000000 --- a/macros/CCEScan/runAllTTForGivenFill.sh~ +++ /dev/null @@ -1,32 +0,0 @@ -#!/bin/bash -# Script to submit all fits to evaluate CCEscans -# run from spinor with: . ./runAllTTForGivenFill.sh 3960 -fill=$1 -echo "****************************************" -echo "Running all fills and sectors in TT" -echo "****************************************" -echo "Full sectors" -echo "Analyzing fill "$fill -echo "****************************************" -for sector in `cat TTsectors.dat` -do -. ./runFullFitLoop.sh TT TTaU $sector $fill -#./runFullFitLoop.sh $sector $fill -done -echo "****************************************" -echo "Inner most" -for sector in `cat TTinner.dat` -do -. ./runFullFitLoop.sh TT TTaU $sector $fill 75 -. ./runFullFitLoop.sh TT TTaU $sector $fill 45 -. ./runFullFitLoop.sh TT TTaU $sector $fill 40 -#./runFullFitLoop.sh $sector $fill 75 -#./runFullFitLoop.sh $sector $fill 45 -#./runFullFitLoop.sh $sector $fill 40 -done - - - -echo "Done!" - -echo "****************************************" \ No newline at end of file diff --git a/macros/CCEScan/runFullFitLoop.sh~ b/macros/CCEScan/runFullFitLoop.sh~ deleted file mode 100755 index 33d9466..0000000 --- a/macros/CCEScan/runFullFitLoop.sh~ +++ /dev/null @@ -1,75 +0,0 @@ -#!/bin/bash -# @Author: Elena Graverini -# @Date: 2015-07-14 14:48:41 -# @Last Modified by: Elena Graverini -# @Last Modified time: 2015-08-02 19:06:01 - -echo "========================================" -echo "Submit for:" -echo " Detector $1" -echo " Layer $2" -echo " Sector $3" -echo " Fill $4" -echo " N Strips $5" -echo " Radius $6" - -# VAL=3 -# -# if [ "$1" = "IT" ] -# then -# VAL=7 -# fi - -source /home/hep/egraveri/.bashrc - -restart_root(){ - if [[ -f /etc/SuSE-release ]]; then - CERN=/app/cern - #ROOTSYS=$CERN/root_v5.34.25/root - ROOTSYS=$CERN/root_v6.04.00/root-6.04.00 - PATH=$PATH:$ROOTSYS/bin - LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ROOTSYS/lib - export LD_LIBRARY_PATH=$ROOTSYS/lib:$PYTHONDIR/lib:$LD_LIBRARY_PATH - export PYTHONPATH=$PYTHONPATH:$ROOTSYS/lib - else - ROOTSYS=/cvmfs/lhcb.cern.ch/lib/lcg/app/releases/ROOT/5.34.10/x86_64-slc6-gcc48-opt/root - PATH=$PATH:$ROOTSYS/bin - LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ROOTSYS/lib - cd $ROOTSYS - source $ROOTSYS/bin/thisroot.sh - cd - > /dev/null - fi -} - -restart_root - -#echo "Sector is " $sector " and fill is " $fill - -# Preparation -kinit -l7d egraveri@CERN.ch -k -t /home/hep/egraveri/torqueJobs/private/egraveri.keytab -aklog -#. /cvmfs/lhcb.cern.ch/lib/lhcb/LBSCRIPTS/LBSCRIPTS_v8r4p1/InstallArea/scripts/SetupProject.sh Vetra v15r0 -export LHCBSTYLE=1 -export DISK=/disk/data1/hep/elena - -# Run all analysis -# Run Landau fits -for i in {0..65} -do - echo "Landau fit for calibration step " $i - #./landauFit -b -d $1 -l $2 -s $3 -f $4 -c $i -v $VAL -ss -r $5 - ./landauFit -b -d $1 -l $2 -s $3 -f $4 -c $i -v $5 -ss -r $6 - #./landauFit -f $fill -s $sector -c $i -b -done -# Run pulse fits -for j in {0..10} -do - echo "Pulse fit for voltage step " $j - #./pulseFit -f $fill -s $sector -c $j -b - #./pulseFit -b -d $1 -l $2 -s $3 -f $4 -c $i -v $VAL -ss -r $5 - ./pulseFit -b -d $1 -l $2 -s $3 -f $4 -c $i -v $5 -ss -r $6 -done -# Run voltage fit -#./voltageFit_spline -f $fill -s $sector -b -#./voltageFit -b -d $1 -l $2 -s $3 -f $4 -v $VAL -ss -r $5" -./voltageFit -b -d $1 -l $2 -s $3 -f $4 -v $5 -ss -r $6" diff --git a/macros/CCEScan/tmp_scripts/test.sh~ b/macros/CCEScan/tmp_scripts/test.sh~ deleted file mode 100644 index e69de29..0000000 --- a/macros/CCEScan/tmp_scripts/test.sh~ +++ /dev/null diff --git a/macros/CCEScan/voltageFit b/macros/CCEScan/voltageFit deleted file mode 100755 index 74a0fef..0000000 --- a/macros/CCEScan/voltageFit +++ /dev/null Binary files differ diff --git a/macros/CCEScan/voltageFit_noPulse b/macros/CCEScan/voltageFit_noPulse deleted file mode 100755 index 9acb9f1..0000000 --- a/macros/CCEScan/voltageFit_noPulse +++ /dev/null Binary files differ diff --git a/macros/CCEScan/voltageFit_spline b/macros/CCEScan/voltageFit_spline deleted file mode 100755 index 74a0fef..0000000 --- a/macros/CCEScan/voltageFit_spline +++ /dev/null Binary files differ diff --git a/macros/leakageCurrent/getPredictions.pyc b/macros/leakageCurrent/getPredictions.pyc deleted file mode 100644 index a68f9ab..0000000 --- a/macros/leakageCurrent/getPredictions.pyc +++ /dev/null Binary files differ diff --git a/macros/runCond/deplVandLumiExtrac b/macros/runCond/deplVandLumiExtrac deleted file mode 100755 index 88d2b98..0000000 --- a/macros/runCond/deplVandLumiExtrac +++ /dev/null Binary files differ diff --git a/macros/runCond/fluenceExtract b/macros/runCond/fluenceExtract deleted file mode 100755 index 2478b99..0000000 --- a/macros/runCond/fluenceExtract +++ /dev/null Binary files differ diff --git a/macros/runCond/flukaExtract b/macros/runCond/flukaExtract deleted file mode 100755 index c4cc251..0000000 --- a/macros/runCond/flukaExtract +++ /dev/null Binary files differ diff --git a/macros/runCond/lumiExtract b/macros/runCond/lumiExtract deleted file mode 100755 index 364bd13..0000000 --- a/macros/runCond/lumiExtract +++ /dev/null Binary files differ diff --git a/macros/runCond/tempExtract b/macros/runCond/tempExtract deleted file mode 100755 index ddea8b8..0000000 --- a/macros/runCond/tempExtract +++ /dev/null Binary files differ diff --git a/scripts/checkPulseData.py b/scripts/checkPulseData.py index 7bd16a8..d826906 100644 --- a/scripts/checkPulseData.py +++ b/scripts/checkPulseData.py @@ -54,7 +54,7 @@ return np.array(prune(pulse)) -def monotonic(x): +def is_monotonic(x): if not len(x): return False dx = np.diff(x) @@ -62,9 +62,6 @@ def prune(x): - #for element in x: - # if element is None: - # x.remove(element) return [element for element in x if element is not None] @@ -95,7 +92,24 @@ return self.substituted_str +class Table(dict): + '''Dictionary tweaked for conversion into pandas DataFrame + + Dictionary conversion into DataFrames is faster than + appending rows to a DataFrame, as pointed out here: + http://stackoverflow.com/a/17496530/3324012 + ''' + def __init__(self, *keylist): + super(dict, self).__init__() + self.__keylist = keylist + for key in keylist: + self[key] = [] + + def append(self, *vals): + assert(len(vals) == len(self.__keylist)) + for i, v in enumerate(vals): + self[self.__keylist[i]].append(v) if __name__ == '__main__': @@ -134,14 +148,9 @@ not_found = StringTemplate('{LAY}/{SEC} on fill {FILL}: calibration step {CS} not found for ns={NS}') not_found = not_found.format(LAY=layer) - # Init data frame - data_frame = pd.DataFrame(columns=['Detector', - 'Sector', - 'Fill', - 'N strips', - 'V bias', - 'MPVs']) - + # Init dictionary for data frame + table = Table('Detector', 'Sector', 'Fill', 'N strips', 'V bias', 'MPVs') + # Loop on data for fill in fills: for sector in sectors: @@ -163,17 +172,10 @@ SEC=sector, CS=cstep)) pulse.append(None) - #pulse = build_pulse(pulse) - if monotonic(build_pulse(pulse)): - # Append this occurrence to the data frame - # Probably it would be faster to do as suggested here: - # http://stackoverflow.com/a/17496530/3324012 - data_frame.loc[len(data_frame)] = [detector, - sector, - fill, - ns, - vmap[vstep], - pulse] + pulse = build_pulse(pulse) + # Append this occurrence to the data frame dictionary + table.append(detector, sector, fill, ns, vmap[vstep], pulse) + if is_monotonic(pulse): if VERBOSE: pulse = [round(v,2) for v in pulse] print(warning.format(NS=ns, @@ -182,4 +184,11 @@ V=vmap[vstep], PULSE=pulse)) - # datafile.Close() + datafile.Close() + # Create the pandas DataFrame from the dictionary + data_frame = pd.DataFrame.from_dict(table) + # Select the instances where pulse data are monotonic + monotonics = data_frame[data_frame['MPVs'].apply(is_monotonic)] + + import pickle + pickle.dump(monotonics, open('monotonics.pkl', 'wb'))