diff --git a/macros/CCEScan/batch/base.sh b/macros/CCEScan/batch/base.sh index 5c3a886..05f6d7a 100644 --- a/macros/CCEScan/batch/base.sh +++ b/macros/CCEScan/batch/base.sh @@ -4,7 +4,6 @@ # export DISK=/disk/data1/hep/elena export DISK=/disk/data12/lhcb/STAgeingData export CCEHOME=/home/hep/egraveri/STAging - # restart_root(){ # if [[ -f /etc/SuSE-release ]]; then # CERN=/app/cern @@ -12,6 +11,7 @@ # #ROOTSYS=$CERN/root_v6.06.00/root-6.06.00 # #ROOTSYS=$CERN/root_v6.02.10/root-6.02.10 # ROOTSYS=$CERN/root_v6.08.06 + # PATH=$PATH:$ROOTSYS/bin # export PATH=$PATH:$ROOTSYS/bin # LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ROOTSYS/lib @@ -26,7 +26,6 @@ # cd - > /dev/null # fi # } - restart_root(){ # if [[ -f /etc/SuSE-release ]]; then CERN=/app/cern @@ -48,5 +47,4 @@ # cd - > /dev/null # fi } - restart_root diff --git a/macros/CCEScan/batch/submitLandau.sh b/macros/CCEScan/batch/submitLandau.sh old mode 100644 new mode 100755 index a1dfb4b..0bbfbf5 --- a/macros/CCEScan/batch/submitLandau.sh +++ b/macros/CCEScan/batch/submitLandau.sh @@ -9,11 +9,11 @@ # TORQUE # to remove erroneous jobs: # for job in `qstat | grep ${USER} | awk '{print $1;}'`; do qdel ${job}@farm-ce; done - swsource=$CCEHOME/macros/CCEScan scripts=$swsource/batch logs=$DISK/logs/ST/CCEScan + #vals=(3 5 7) vals=(5 7) rads=(75 45 40) @@ -25,6 +25,7 @@ #onefill=(5448) #onefill=(4643 4856 5162 5448) + echo "Submitting TT:" for fill in `cat $swsource/Fills.dat` #for fill in "${onefill[@]}" @@ -34,6 +35,7 @@ do echo " Sector "$sector scriptname="landauFit_TT_f${fill}_s${sector}_allSteps_allVal" + mkdir -p $scripts/temp_scripts #Added by Michele cp $scripts/base.sh $scripts/temp_scripts/${scriptname}.slurm chmod 777 $scripts/temp_scripts/${scriptname}.slurm for i in {0..65} @@ -43,7 +45,7 @@ echo "${swsource}/landauFit -d TT -l TTaU -s ${sector} -f ${fill} -c ${i} -v ${val} -b -ss" >> $scripts/temp_scripts/${scriptname}.slurm done done - sbatch -p standard $scripts/temp_scripts/${scriptname}.slurm -o ${logs}/${scriptname}.log -e ${logs}/${scriptname}.log + sbatch -p standard --ntasks-per-node=2 $scripts/temp_scripts/${scriptname}.slurm -o ${logs}/${scriptname}.log -e ${logs}/${scriptname}.log #sleep 2 done done @@ -69,7 +71,7 @@ done done done - sbatch -p standard $scripts/temp_scripts/${scriptname}.slurm -o ${logs}/${scriptname}.log -e ${logs}/${scriptname}.log + sbatch -p standard --ntasks-per-node=2 $scripts/temp_scripts/${scriptname}.slurm -o ${logs}${scriptname}.log -e ${logs}${scriptname}.log #sleep 2 done done @@ -94,7 +96,7 @@ do echo "${swsource}/landauFit -d IT -l T3X2 -s ${sector} -f ${fill} -c ${i} -v ${val} -b -ss" >> $scripts/temp_scripts/${scriptname}.slurm done - sbatch -p standard $scripts/temp_scripts/${scriptname}.slurm -o ${logs}/${scriptname}.log -e ${logs}/${scriptname}.log + sbatch -p standard --ntasks-per-node=2 $scripts/temp_scripts/${scriptname}.slurm -o ${logs}/${scriptname}.log -e ${logs}/${scriptname}.log sleep 1 done done @@ -108,6 +110,7 @@ do echo " Fill "$fill val=7 + scriptname="landauFit_ITinner_s${sector}_f${fill}_allSteps_v${val}_r175" cp $scripts/base.sh $scripts/temp_scripts/${scriptname}.slurm chmod 777 $scripts/temp_scripts/${scriptname}.slurm @@ -127,7 +130,9 @@ done sbatch -p standard $scripts/temp_scripts/${scriptname}.slurm -o ${logs}/${scriptname}.log -e ${logs}/${scriptname}.log sleep 1 + done done +#rm -r $scripts/temp_scripts diff --git a/macros/CCEScan/pulseFit.C b/macros/CCEScan/pulseFit.C index ee91c6e..f10bf61 100755 --- a/macros/CCEScan/pulseFit.C +++ b/macros/CCEScan/pulseFit.C @@ -92,7 +92,8 @@ // Executable method int pulseFit(int argc, char* argv[]){ - TVectorD v_res(11); + //TVectorD v_res(11); + TVectorD v_res(12); //MIchele gStyle->SetPadLeftMargin(0.20); gStyle->SetTitleOffset(1.5,"Y"); @@ -192,7 +193,7 @@ // Directory and file name for the input data TString dn_data(Form("%s/data/ST/Aging" ,diskVar)); - TString fn_data(Form("%s/CCEScan.root", + TString fn_data(Form("%s/CCEScan.root", dn_data.Data())); TString fn_dummy(Form("%s/check",dn_data.Data())); @@ -267,6 +268,7 @@ double chi2ndf = fu_pulse->GetChisquare()/fu_pulse->GetNDF(); Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF()); + // Make PDG-like correction of the uncertainties //Commented out by Michele - 8/31 2017 @@ -281,6 +283,92 @@ ge_pulse = new TGraphErrors(vr_time_val,vr_adc_val,vr_time_err,vr_adc_err); ge_pulse->Fit(fu_pulse,"M0Q","0"); Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF()); + + + + //Evaluation of the systematic error on the integral of the pulse shape using a polinomial + //Define the interval were most likely the zeros should be + //TO DO: Set min and max from the time vector + double low_limit = -50; + double high_limit = 50; + //HERE + //POLY3 + + cout << "================ poly3 ====================" << endl<SetLineColor(1); + ge_pulse->Fit("poly3","M"); + Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly3->GetChisquare(),poly3->GetNDF()); + + + std::vector v_max_p3; + std::vector v_time_srtd; + bool b_max_p3 = true; + + //Copy the vector of interest + for (int i=0; iGetMaximum(v_time_srtd[i], v_time_srtd[i+1]); + + //Checks on the maximum -> do not accept maxima at the boundary + if( (max_tmp <= poly3->Eval(v_time_srtd[i])) || (max_tmp <= poly3->Eval(v_time_srtd[i+1])) ){ + continue; + }else{ + v_max_p3.push_back(max_tmp); + } + } + + if (v_max_p3.size() == 1){ + + //Final check + + cout << "Found the maximum of interest!It is "<SetLineColor(4); + ge_pulse->Fit("poly4","M0Q","0"); + Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly4->GetChisquare(),poly4->GetNDF()); + + double max_poly4 = poly4->GetMaximum(low_limit, high_limit); + + + //POLY5 + cout << "================ poly5 =================" << endl<SetLineColor(6); + ge_pulse->Fit("poly5","M"); + Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly5->GetChisquare(),poly5->GetNDF()); + + double max_poly5 = poly5->GetMaximum(low_limit, high_limit); + */ + //MIchele - End + /* + TF1* fu_pulse2 = STTool::halfGauss2; + ge_pulse->Fit(fu_pulse2,"MQ","0"); + Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse2->GetChisquare(),fu_pulse2->GetNDF()); + */ + TVirtualFitter *fit_volt = TVirtualFitter::GetFitter(); Double_t* cov_volt = fit_volt->GetCovarianceMatrix(); @@ -305,20 +393,44 @@ 2.0*cov_volt[2]/(tau_val*fu_pulse->GetParameter(0)) ); + + //Michele 3/9/2017 + //Estimate the systematics: # if a polynomial with zeros is found take as sys the semidifference + // between its integral and the integral obtained from the fitted theoretical model + // # if no polynomial behaved properly set the systematic error to -999 + //Michele - end + double max_sys = -999; + if (b_max_p3 == true){ + max_sys = TMath::Abs(max_val - v_max_p3[0])/2.; + } + + - mg_pulse->Add(ge_pulse,"pe"); + mg_pulse->Add(ge_pulse,"PE"); Info("pulseFit","pulse Height: %8.4f+/-%8.4f",max_val,max_err); Info("pulseFit"," tau: %8.4f+/-%8.4f ns",tau_val,tau_err); - Info("pulseFit"," charge eq: %8.4f+/-%8.4f",int_val,int_err); - + + //Michele + //TO DO: control + if(b_max_p3 == true){ + Info("pulseFit"," charge eq: %8.4f +/- %8.4f(stat) +/- %8.4f (sys)",max_val,max_err, max_sys); + }else{ + + Info("pulseFit"," charge eq: %8.4f +/- %8.4f(stat) +/- No sys eval.", max_val,max_err); + } + // // Plot pulse TCanvas* c_pul = new TCanvas("c_pul","Pulse Canvas",800,600); mg_pulse->Draw("A"); fu_pulse->Draw("Same"); + //fu_pulse2->Draw("Same"); + poly3->Draw("Same"); + //poly4->Draw("Same"); + //poly5->Draw("Same"); mg_pulse->Draw(""); mg_pulse->GetXaxis()->SetLimits(-30.0,30.0); mg_pulse->GetXaxis()->SetTitle("#delta#it{t} [ns]"); @@ -332,6 +444,8 @@ // 0,1: Maximal signal height v_res(0) = max_val; v_res(1) = max_err; + //v_res(1) = max_sys; + // 2,3: tau (width of pulse) v_res(2) = tau_val; v_res(3) = tau_err; @@ -342,7 +456,7 @@ v_res(8) = cov_volt[1]; v_res(9) = cov_volt[2]; v_res(10)= cov_volt[5]; - + v_res(11) = max_sys; // Write the data and the plots diff --git a/macros/CCEScan/voltageFit_spline.C b/macros/CCEScan/voltageFit_spline.C index 31a1e49..50d3344 100755 --- a/macros/CCEScan/voltageFit_spline.C +++ b/macros/CCEScan/voltageFit_spline.C @@ -130,7 +130,7 @@ } - + std::vector v_volt_val; std::vector v_volt_err; @@ -177,7 +177,7 @@ TString dn_data(Form("%s/data/ST/Aging" ,diskVar)); - TString fn_data(Form("%s/CCEScan.root", + TString fn_data(Form("%s/CCEScan.root",//Michele for TESTING dn_data.Data())); TString fn_dummy(Form("%s/check",dn_data.Data())); @@ -222,9 +222,14 @@ fill)); - - + FILE* myfile; + TString myfile_nm(Form("%s/noSystematic_val%d.txt",dn_graph.Data(),val)); + cout << myfile_nm << endl; + myfile =fopen(myfile_nm,"w"); + fprintf(myfile, "List of the calibration steps with no estimated systematic on pulse shape\n\n"); + for (int i=0; i<11; i++) { + TString vn_input(Form("%s/v_pulse%d_val%d",dn_input.Data(),i,val)); if (radius>0.0) { vn_input +=Form("_r%d",(int)radius); @@ -258,7 +263,10 @@ double height_val = v_input(0); double height_err = v_input(1); - + double height_sys = v_input(11);//Michele + + /* Michele - Use maximum value, not integral + //This is a fake if... if (true || det.EqualTo("IT")) { height_val /= 0.1306; // transform back to the amplitude height_err /= 0.1306; // transform back to the amplitude @@ -267,22 +275,28 @@ TMath::Power(v_input(3)/v_input(2),2)+ 2*v_input(8)/(height_val*v_input(2))); } - + */ - Printf(" %8.2f V: %8.4f+/-%8.4f",volt,height_val,height_err); + Printf(" %8.2f V: %8.4f+/-%8.4f(stat) +/-%8.4f (sys)",volt,height_val,height_err, height_sys);//MIchele // Load the input values in to vector - v_adc_val.push_back(height_val); - v_adc_err.push_back(height_err); - v_volt_val.push_back(volt); - v_volt_err.push_back(0.0); - + + if (height_sys == -999){ + fprintf(myfile," det=%s lay=%s sect=%d fill=%d val=%d rad=%f\n", det.Data(),lay.Data(),sector,fill,val,radius); + + }else{ + + v_adc_val.push_back(height_val); + v_adc_err.push_back(height_err + height_sys); //Putting systematic and statistical togther, temporary fix + v_volt_val.push_back(volt); + v_volt_err.push_back(0.0); + } } - + fclose(myfile); if (v_adc_val.size()<3) { Error("voltageFit_spline","Too few (%d) data points",v_adc_val.size()); @@ -307,7 +321,7 @@ TMultiGraph* mg_volt = new TMultiGraph("mg_volt","Voltage multigraph"); - // Extract the ratio values for the sampling point + // Extract the ratio values for the sampling point TVectorD* vp_ratio = (TVectorD*)f_ratio->Get(Form("v_ratio_%s",lay.Data())); TVectorD* vp_ratio_err = (TVectorD*)f_ratio->Get(Form("v_ratio_err_%s",lay.Data())); TVectorD* vp_secRatio = (TVectorD*)f_ratio->Get(Form("v_sector_%s",lay.Data())); @@ -332,22 +346,22 @@ ratio_err = STTool::TTratio_err; } - /* Elena - // if (false && det.EqualTo("TT")) { - */ - // Using ratio per sector - if (true) { - int ratioInd = getClosestIndex(v_secRatio,(double)sector); - if (!(ratioInd<0) && (int)v_secRatio(ratioInd)==sector && true) { - Info("voltageFit_spline","Using specific ratio for sector %d",sector); - ratio_val = v_ratio(ratioInd); - ratio_err = v_ratio_err(ratioInd); + // Elena - Cancelled by Michele because of what said in l.159 in the Radiation Damage paper + if (false && det.EqualTo("TT")) { + + // Using ratio per sector + if (true) { + int ratioInd = getClosestIndex(v_secRatio,(double)sector);//Gets the index in v_secRatio that has a double that is closest to sector + if (!(ratioInd<0) && (int)v_secRatio(ratioInd)==sector && true) {//if this is not senseless and the sector is correct + Info("voltageFit_spline","Using specific ratio for sector %d",sector); + ratio_val = v_ratio(ratioInd); + ratio_err = v_ratio_err(ratioInd); + } } } - // Define fit range double fitMax = 510.0; double fitMin = 280.0;