//
// deplTool.h
// Tools and data for depletion voltage analysis
//
// Created by Christian Elsasser on 01.05.12.
// University of Zurich, elsasser@cern.ch
// Copyright only for commercial use
//
// General include
#include <iostream>
#include <math.h>
#include <map>
#include "TROOT.h"
#include "TSpline.h"
namespace STTool {
// Cut on tracks
TCut cu_track_TT("(TrChi2/TrNDoF<3.0) &&"
"(GhostP<0.01+0.1/1.6*TrChi2/TrNDoF) &&"
"(GhostP<0.01-0.1/3.4*(TrChi2/TrNDoF-5.0))");
TCut cu_track_IT("(GhostP<0.1*TrChi2/TrNDoF) &&"
"(GhostP<0.2-0.1*(TrChi2/TrNDoF-5.0))");
// Define first run of fill
static std::map<int,int> zeroRunOfFill(){
std::map<int,int> runMap;
runMap[1020] = 69606;
runMap[1616] = 87147;
runMap[1944] = 95936;
runMap[2083] = 101357;
runMap[2252] = 104091;
runMap[2472] = 111267;
runMap[2797] = 120006;
runMap[3108] = 129494;
runMap[3478] = 135605;//TO BE UPDATED!!
runMap[3820] = 153583;
runMap[3960] = 156887;
runMap[4518] = 166262;
runMap[4638] = 168092;
runMap[4643] = 168257;
runMap[4856] = 173040;
runMap[5162] = 181247;
runMap[5448] = 185508;
return runMap;
}
std::map<int,int> zeroRunMap = zeroRunOfFill();
// Value of one time step used in sampling time change
double timestep = 0.10417; // ns
// Calibration step residual -> time step (i.e. step%6)
static std::map<int,double> createTimeMap(){
std::map<int,double> timeMap;
timeMap[0] = 0.0;
timeMap[1] = -60.0;
timeMap[2] = -30.0;
timeMap[3] = 30.0;
timeMap[4] = 60.0;
timeMap[5] = 90.0;
return timeMap;
}
std::map<int,double> timeMap = createTimeMap();
//Added by Michele Atzeni to take into account the changes in the timing steps after fill 2252
// Calibration step residual -> time step (i.e. step%6)
static std::map<int,double> createTimeMap_post(){
std::map<int,double> timeMap_post;
timeMap_post[0] = 0.0;
timeMap_post[1] = -90.0;
timeMap_post[2] = -45.0;
timeMap_post[3] = 45.0;
timeMap_post[4] = 90.0;
timeMap_post[5] = 135.0;
return timeMap_post;
}
std::map<int,double> timeMap_post = createTimeMap_post();
// Calibration step residual -> TT voltage in volts (i.e. step/6)
static std::map<int,double> createVoltMapTT(){
std::map<int,double> voltMapTT;
voltMapTT[0] = 400.0;
voltMapTT[1] = 350.0;
voltMapTT[2] = 300.0;
voltMapTT[3] = 250.0;
voltMapTT[4] = 225.0;
voltMapTT[5] = 200.0;
voltMapTT[6] = 175.0;
voltMapTT[7] = 150.0;
voltMapTT[8] = 125.0;
voltMapTT[9] = 100.0;
voltMapTT[10] = 60.0;
voltMapTT[11] = 20.0; // Not used, could remove
voltMapTT[12] = 20.0; // Not used, could remove
return voltMapTT;
}
std::map<int,double> voltMapTT = createVoltMapTT();
// Calibration step residual -> IT voltage in volts (i.e. step/6)
static std::map<int,double> createVoltMapIT(){
std::map<int,double> voltMapIT;
voltMapIT[0] = 300.0;
voltMapIT[1] = 200.0;
voltMapIT[2] = 170.0;
voltMapIT[3] = 140.0;
voltMapIT[4] = 120.0;
voltMapIT[5] = 105.0;
voltMapIT[6] = 90.0;
voltMapIT[7] = 75.0;
voltMapIT[8] = 60.0;
voltMapIT[9] = 40.0;
voltMapIT[10] = 20.0;
voltMapIT[11] = 20.0; // Not used, could remove
voltMapIT[12] = 20.0; // Not used, could remove
return voltMapIT;
}
std::map<int,double> voltMapIT = createVoltMapIT();
// Ratio values
double TTratio_val = 0.953; // with time correction
double TTratio_err = 0.019;
double ITratio_val = 0.956; // with time correction
double ITratio_err = 0.017;
// FLUKA conversion factor
double FLUKAConvFac = 7.2*1e4;
// Returns time step for a given calibration step
//Modified to take into account the change in the time steps used
double GetTimeStep(int step, int fill){
if (step<0) {
Error("GetTimeStep","Calibration step is smaller than zero, returning NaN");
return sqrt(-1);
}
if( fill > 2252){
return timeMap_post[step%6];
}
else{
return timeMap[step%6];
}
}
// Returns time in ns for a given calibration step
//Modified to take into account the change in the time steps used
double GetTime(int step, int fill){
if (step<0) {
Error("GetTime","Calibration step is smaller than zero, returning NaN");
return sqrt(-1);
}
if( fill > 2252){
return timeMap_post[step%6]*timestep;
}
else{
return timeMap[step%6]*timestep;
}
}
// Returns TT voltage in V for a given calibration step
double GetTTVoltage(int step){
if (step<0 || step/6>12) {
Error("GetTTVoltage","Calibration step is out of reach, returning NaN");
return sqrt(-1);
}
return voltMapTT[step/6];
}
// Returns IT voltage in V for a given calibration step
double GetITVoltage(int step){
if (step<0 || step/6>12) {
Error("GetITVoltage","Calibration step is out of reach, returning NaN");
return sqrt(-1);
}
return voltMapIT[step/6];
}
// Return calibration step for a given time step and TT voltage
int GetTTCalibStep(double time, double voltage, int fill){
std::map<int,double> timeMap_sw;
if(fill>2252){
timeMap_sw = timeMap_post;
}
else{
timeMap_sw = timeMap;
}
std::map<int,double>::const_iterator timeEnd = timeMap_sw.end();
std::map<int,double>::const_iterator voltEnd = voltMapTT.end();
for (std::map<int,double>::const_iterator itTime = timeMap_sw.begin();
itTime != timeEnd; itTime++) {
for (std::map<int,double>::const_iterator itVolt = voltMapTT.begin();
itVolt != voltEnd; itVolt++) {
if (itTime->second == time && itVolt->second == voltage) {
return itVolt->first*6+itTime->first;
}
}
}
Warning("GetTTCalibStep","No corresponding calibration step found, returning -1");
return -1;
}
// Return calibration step for a given time step and IT voltage
int GetITCalibStep(double time, double voltage, int fill){
std::map<int,double> timeMap_sw;
if(fill>2252){
timeMap_sw = timeMap_post;
}
else{
timeMap_sw = timeMap;
}
std::map<int,double>::const_iterator timeEnd = timeMap_sw.end();
std::map<int,double>::const_iterator voltEnd = voltMapIT.end();
for (std::map<int,double>::const_iterator itTime = timeMap_sw.begin();
itTime != timeEnd; itTime++) {
for (std::map<int,double>::const_iterator itVolt = voltMapIT.begin();
itVolt != voltEnd; itVolt++) {
if (itTime->second == time && itVolt->second == voltage) {
return itVolt->first*6+itTime->first;
}
}
}
Warning("GetITCalibStep","No corresponding calibration step found, returning -1");
return -1;
}
// Half Gaussian for Pulse
TF1* createHalfGauss(){
TF1* func = new TF1("halfGauss",
"[0]*"
"(TMath::Power((x-[1])/[2],2)/2.0-"
"TMath::Power((x-[1])/[2],3)/6.0)*"
"TMath::Exp(-(x-[1])/[2])",-30.0,30.0);
func->SetParNames("A_{0}","t_{0}","#tau");
func->SetParameters(42,-15.0,10.0);
func->SetParLimits(0, 10.0,800.0);
func->SetParLimits(1,-50.0, 10.0);
func->SetParLimits(2, 1.0, 60.0);
func->SetNpx(1000);
return func;
}
TF1* halfGauss = createHalfGauss();
// Level as V_d
double VdeplLevel = 0.8;
// Create sigmoid function 1 for Voltage
Double_t func_Sig1(Double_t *x, Double_t *par)
// par[0] = V_d, par[1] = V_0, par[2] = A_0, par[3] = per
{
double per = par[3];
double s = TMath::Log(per/(1.0-per));
double exp = par[2]/(1+TMath::Exp(-s*(x[0]-par[1])/(par[0]-par[1])));
double lin = par[2]/2.0*(1+s/2.0*(x[0]-par[1])/(par[0]-par[1]));
return TMath::Min(exp,lin);
}
TF1* createSigmoid1(){
TF1* func = new TF1("Sigmoid1",func_Sig1,-100.0,800.0,5);
func->SetParNames("V_{depl}","V_{0}","A_{0}","level");
func->SetParameters(250.0,160.0,42,VdeplLevel);
func->SetParLimits(0, 0.0, 300.0);
func->SetParLimits(1,-20.0, 200.0);
func->SetParLimits(2, 10.0,800.0);
func->SetParLimits(4, 0.0,1.0);
func->FixParameter(3,VdeplLevel);
func->SetNpx(1000);
return func;
}
TF1* sigmoid1 = createSigmoid1();
// Create sigmoid function 2 for Voltage
Double_t func_Sig2(Double_t *x, Double_t *par)
// par[0] = V_d, par[1] = V_0, par[2] = A_0, par[3] = per
{
double x0 = par[2]/par[0]+par[0]/(4*par[1]);
double x1 = x0 - par[0]/(2*par[1]);
if(x[0]<x1){
return par[0]*x[0];
}
if(x[0]<x0){
return par[2]-(x[0]-x0)*(x[0]-x0)*par[1];
}
return par[2];
}
TF1* createSigmoid2(){
TF1* func = new TF1("Sigmoid2",func_Sig2,-100.0,800.0,3);
func->SetParNames("C_{0}","B_{0}","A_{0}");
func->SetParameters(3.0,0.5,42);
func->SetParLimits(0, 0.0, 10.0);
func->SetParLimits(1, 0.0, 20.0);
func->SetParLimits(2,10.0, 100.0);
func->SetNpx(1000);
return func;
}
TF1* sigmoid2 = createSigmoid2();
// Create sigmoid function 3 for Voltage
Double_t func_Sig3(Double_t *x, Double_t *par)
// par[0] = V_d, par[1] = V_0, par[2] = A_0, par[3] = per
{
if(x[0]<par[0]){
return -TMath::Sqrt(par[1]*(x[0]-par[0])*(x[0]-par[0])+par[2])+par[3]+TMath::Sqrt(par[2]);
}
return par[3];
}
TF1* createSigmoid3(){
TF1* func = new TF1("Sigmoid2",func_Sig3,-100.0,800.0,4);
func->SetParNames("X_{0}","C_{0}","B_{0}","A_{0}");
func->SetParameters(250.0,0.5,2.0,42);
func->SetParLimits(0, 0.0, 500.0);
func->SetParLimits(1, 0.0, 200.0);
func->SetParLimits(2, 0.0, 200.0);
func->SetParLimits(3,10.0, 100.0);
func->SetNpx(1000);
return func;
}
TF1* sigmoid3 = createSigmoid3();
// Create sigmoid function 4 for Voltage
Double_t func_Sig4(Double_t *x, Double_t *par)
// par[0] = V_d, par[1] = A_0, par[2] = per
{
double per = par[2];
double s = TMath::Log(per/(1.0-per));
double V0 = 2.0*par[0]/(2.0+s);
double exp = par[1]/(1+TMath::Exp(-s*(x[0]-V0)/(par[0]-V0)));
double lin = par[1]/2.0*(1+s/2.0*(x[0]-V0)/(par[0]-V0));
return TMath::Min(exp,lin);
}
TF1* createSigmoid4(){
TF1* func = new TF1("Sigmoid4",func_Sig4,-100.0,800.0,3);
func->SetParNames("V_{depl}","A_{0}","level");
func->SetParameters(250.0,160.0,42,VdeplLevel);
func->SetParLimits(0, 0.0, 300.0);
func->SetParLimits(1, 10.0,800.0);
func->FixParameter(2,VdeplLevel);
func->SetNpx(1000);
return func;
}
TF1* sigmoid4 = createSigmoid4();
// Create function describing the ballistic deficit
Double_t func_BD(Double_t *x, Double_t *par){
// par[0] = V_d, par[1]= Q_0, par[2] = lambda_0, par[3] = lambda_1, par[4] = varepsilon_s
// par[5] = d (fixed)
double w = par[5]*TMath::Min(1.0,TMath::Sqrt(x[0]/par[0]));
//Printf("V_d = %f",par[0]);
//Printf("%f",TMath::Exp(-(par[5]-w)/par[2]));
double y = par[1]/par[5]*TMath::Exp(-(par[5]-w)/par[2]);
//Printf("y = %f",y);
int nInt = 1000;
double dt = 1.0*w/nInt;
double t = dt/2.0;
double integral = 0.0;
for (int i=0; i<nInt; i++) {
double eps = 1.0/par[5]*(TMath::Sqrt(x[0]*par[0])-t*par[0]/par[5]);
double lam = par[2]+par[3]*TMath::Min(1.0,eps/par[4]);
//Printf("w = %f, dt = %f, eps = %f",w,dt,eps);
integral += TMath::Exp(-t/lam);
t += dt;
}
return y*integral*dt;
}
TF1* createBD(){
TF1* func = new TF1("BDFunc",func_BD,-100.0,800.0,6);
func->SetParNames("V_{d}","Q_{0}","#lambda_{0}",
"#lambda_{1}","#varepsilon_{s}","d");
func->SetParameters(200.0,40.0,500.0,1000.0,0.8,500.0);
func->SetParLimits(0,20.0, 500.0);
func->SetParLimits(1, 0.0, 200.0);
func->SetParLimits(2,190.0, 210.0);
func->SetParLimits(3,50.0, 3000.0);
func->SetParLimits(4, 1.0, 1.5);
func->FixParameter(5, 500.0);
func->SetNpx(1000);
return func;
}
TF1* DBfunc = createBD();
// Constant function
TF1* createConst(){
TF1* func = new TF1("Constfunc","[0]",-100.0,800.0);
func->SetParNames("A_{0}");
func->SetParameter(0,40.0);
func->SetParLimits(0,10.0, 500.0);;
func->SetNpx(1000);
return func;
}
TF1* Constfunc = createConst();
// Generic get X
double GetX(TObject* o,double y,double xmin, double xmax, int steps=1000){
if (strcmp(o->IsA()->GetName(),"TF1")==0) {
TF1* f = (TF1*)o;
return f->GetX(y,xmin,xmax);
}
double x = xmin;
double step = (xmax-xmin)/steps;
TSpline* sp = (TSpline*)o;
for (int i=0; i<steps; i++) {
if ( (sp->Eval(x)<y && sp->Eval(x+step)>y) || (sp->Eval(x)>y && sp->Eval(x+step)<y) ) {
if (TMath::Abs(sp->Eval(x)-y)==0.0) {
return x;
}else if (TMath::Abs(sp->Eval(x+step)-y)==0.0){
return x+step;
}
return (x/TMath::Abs(sp->Eval(x)-y)+(x+step)/TMath::Abs(sp->Eval(x+step)-y))/(1.0/TMath::Abs(sp->Eval(x)-y)+1.0/TMath::Abs(sp->Eval(x+step)-y));
}
x+=step;
}
Warning("GetX","Has not found a suitable point, returning 0.0");
return 0.0;
}
// List of sectors to be excluded
bool IsExcluded(int sec){
// create List
std::vector<int> v_ex;
// TT
// Outermost regions -> low statistics
v_ex.push_back(2593);
v_ex.push_back(2594);
v_ex.push_back(2595);
v_ex.push_back(2596);
v_ex.push_back(2597);
v_ex.push_back(2598);
v_ex.push_back(2599);
v_ex.push_back(2600);
v_ex.push_back(2601);
v_ex.push_back(2602);
v_ex.push_back(2603);
v_ex.push_back(2604);
v_ex.push_back(2605);
v_ex.push_back(2606);
v_ex.push_back(2607);
v_ex.push_back(2608);
v_ex.push_back(2665);
v_ex.push_back(2666);
v_ex.push_back(2667);
v_ex.push_back(2668);
v_ex.push_back(2669);
v_ex.push_back(2670);
v_ex.push_back(2671);
v_ex.push_back(2672);
v_ex.push_back(2673);
v_ex.push_back(2674);
v_ex.push_back(2675);
v_ex.push_back(2676);
v_ex.push_back(2677);
v_ex.push_back(2678);
v_ex.push_back(2679);
v_ex.push_back(2680);
return std::find(v_ex.begin(), v_ex.end(), sec)!=v_ex.end();
}
// Systematics on depletion voltage
double voltSyst_TT = 13.4; // V
double voltSyst_IT = 8.8; // V
};
