Commit 186dc63a authored by Emilio Depero's avatar Emilio Depero
Browse files

analysis update

parent 033dd2ab
......@@ -72,12 +72,13 @@ int main(int argc, char *argv[])
//choose histogram to cast, only 1D supported
vector<TString> histonames;
histonames.push_back("decayPositions_z");
histonames.push_back("decayPositions_z_target");
histonames.push_back("decayPositions_z");
histonames.push_back("decayPositions_z_20deg");
histonames.push_back("decayPositions_z_17deg");
histonames.push_back("decayPositions_z_13deg");
histonames.push_back("decayPositions_z_10deg");
histonames.push_back("decayPositions_z_16deg");
histonames.push_back("decayPositions_z_12deg");
histonames.push_back("decayPositions_z_8deg");
histonames.push_back("decayPositions_z_4deg");
histonames.push_back("decayPositions_z_0deg");
histonames.push_back("angle");
for(unsigned int d(0); d < histonames.size(); ++d)
......
......@@ -8,16 +8,15 @@ export DAQOUTDIR=$MASTERPATH"/run_November2021"
export rotangle="0"
export OUTDIR=$DAOUTDIR"/analysed"
if [ $# -eq 1 ] || [ $# -eq 2 ] || [ $# -eq 3 ]
if [ $# -eq 1 ] || [ $# -eq 2 ]
then
rotangle=$1
distmm1=$2
fi
#optional target output
if [ $# -eq 3 ]
if [ $# -eq 2 ]
then
OUTDIR=$3
OUTDIR=$2
mkdir -p $OUTDIR
fi
......@@ -42,8 +41,18 @@ do
export analfileout=$OUTDIR"/Analysed_"$(basename $file)
#run the reco
echo "running reco for $file"
$ANALEXE $file $analfileout $rotangle $distmm1
$ANALEXE $file $analfileout $rotangle
done
#merge all that is needed
$MASTERPATH/mergerun.sh gate5mum_delay0_en7.5keV.root $OUTDIR/*run{30,31,38,39,40,41,42,43,45,46,47,48,49,50,53,54,55,56,57,58,59}_*root
$MASTERPATH/mergerun.sh gate1mum_delay1_en7.5keV.root $OUTDIR/*run{68,69,75}_*root
$MASTERPATH/mergerun.sh gate1mum_delay4_en5keV.root $OUTDIR/*run{70,73}_*root
$MASTERPATH/mergerun.sh gate1mum_all.root $OUTDIR/*run{68,69}_*root $OUTDIR/*run7?_*root
echo "ALL DONE"
......@@ -2,24 +2,16 @@
if [ $# -lt 2 ]; then
echo "INVALID PARAMETER NUMBER: $#"
echo "usage: ./mergerun <1-run> <2-run> <3-run> ..."
echo "usage: ./mergerun <name-of-output> <1-run> <2-run> <3-run> ..."
exit
fi
#build outname
export outname="Analysed"
for file in $@;
do
filename=$(basename $file)
run=$(echo "$filename" | grep -o -E '[0-9]+')
outname=$outname"-$run"
done
outname=$outname".root"
outname=$1".root"
echo "creating file with name $outname"
hadd -k -j $outname $@
hadd -f -k -j $outname ${@-1}
echo "FINISHED MERGING"
......@@ -36,6 +36,7 @@
#include "TF2.h"
#include "TLegend.h"
#include "TGraph.h"
#include "TParameter.h"
#include "TGraphErrors.h"
#include "TVector.h"
#include "TList.h"
......@@ -122,42 +123,33 @@ void Analyze(int argc, char *argv[])
// define variables to extract data from the g4bl ntuples
Float_t x,y,t,EventID, qmax, rotangle;
if (argc < 3 || argc>5)
if (argc < 3 || argc>4)
{
std::cout << "Usage: ./executable [inputFileName] [outpuFileName] [rotangle (optional, in degree)] [distance MM1 (optional, in mm)]." << std::endl;
std::cout << "Usage: ./executable [inputFileName] [outpuFileName] [rotangle (optional, in degree)]." << std::endl;
return;
}
// distance between two detectors, exactly the same as Jonas run, mechanics was fixed already
float dist_MM12 = 66.5; //mm
// distance of MM1C from center plane
float dist_MM1 = 135.; //measured by mike, 59 mm dist + half of cube
float size_MM = 80.0; // size in mm
// get parameters from command line arguments
TString inputFileName(argv[1]);
TString outputFileName(argv[2]);
TString outputFileName(argv[2]);
int degreeangle;
if ( argc > 4)
{
degreeangle = atoi(argv[3]);
dist_MM1 = atof(argv[4]);
rotangle = degreeangle/180.*TMath::Pi(); //radiant convertion
outputFileName += TString::Format("-rotation%i-distmm1%0.1f.root", degreeangle, dist_MM1);
}
else if(argc > 3)
if(argc > 3)
{
degreeangle = atoi(argv[3]);
rotangle = degreeangle/180.*TMath::Pi(); //radiant convertion
outputFileName += TString::Format("-rotation%i.root", degreeangle);
outputFileName += TString::Format("-rotation%i.root", degreeangle);
}
else
{
degreeangle = 0;
rotangle = 0;
}
}
// distance between two detectors, exactly the same as Jonas run, mechanics was fixed already
float dist_MM12 = 66.5; //mm
// distance of MM1C from center plane
float dist_MM1 = 135.; //measured by mike, 59 mm dist + half of cube
float size_MM = 80.0; // size in mm
std::cout << "PERFORMING ANALYSIS WITH ROTATION ANGLE: " << degreeangle <<" degrees \n";
......@@ -166,7 +158,7 @@ void Analyze(int argc, char *argv[])
MM_Names.push_back("MM2");
// cut angle in degrees
std::vector<float> maxAngle_degree = {20.0,17.0,13,10.0};
std::vector<float> maxAngle_degree = {20.0,17.0,13,10.0,5.0,3.0};
std::vector<float> maxAngle_rad;
for (int i=0;i<maxAngle_degree.size();i++){maxAngle_rad.push_back(maxAngle_degree[i]/180.*TMath::Pi()); std::cout << "angle in [rad]:" << maxAngle_degree[i]/180.*TMath::Pi() << "\n";}
......@@ -204,20 +196,16 @@ void Analyze(int argc, char *argv[])
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
std::vector<TH2D*> projectedMaps;
/// concidences detected
TH2D* projectedHeatMap_20deg = new TH2D("projectedHeatMap_20deg","reconstructed Map, angular cut 10 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
TH2D* projectedHeatMap_17deg = new TH2D("projectedHeatMap_17deg","reconstructed Map, angular cut 8 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
TH2D* projectedHeatMap_13deg = new TH2D("projectedHeatMap_13deg","reconstructed Map, angular cut 6 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
TH2D* projectedHeatMap_10deg = new TH2D("projectedHeatMap_10deg","reconstructed Map, angular cut 4 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
projectedMaps.push_back(projectedHeatMap_20deg);
projectedMaps.push_back(projectedHeatMap_17deg);
projectedMaps.push_back(projectedHeatMap_13deg);
projectedMaps.push_back(projectedHeatMap_10deg);
TH2D* projectedZoom_20deg = new TH2D("projectedZoom_20deg","reconstructed Map, angular cut 10 degrees; x [mm] ; z [mm] ",110,-14,14,125,-11,20);
TH2D* projectedZoom_17deg = new TH2D("projectedZoom_17deg","reconstructed Map, angular cut 8 degrees; x [mm] ; z [mm] ",110,-14,14,125,-11,20);
TString histname, histtitle;
double histrange[2] = {-100, 100};
unsigned int histbin = 200;
for(unsigned int i(0); i < maxAngle_degree.size(); ++i)
{
histname.Form("projectedHeatMap_%ideg", maxAngle_degree[i]);
histtitle.Form("reconstructed Map, angular cut %i degrees; x [mm] ; z [mm] ", maxAngle_degree[i]);
projectedMaps.push_back(new TH2D(histname, histtitle, histbin, histrange[0], histrange[1], histbin, histrange[0], histrange[1]));
}
TH2D* singleHeatMap_1 = new TH2D("singleHeatMap_1","histogram for MM2; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
TH2D* singleHeatMap_2 = new TH2D("singleHeatMap_2","histogram for MM1; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
......@@ -225,14 +213,19 @@ void Analyze(int argc, char *argv[])
TH2D* qmaxMap_1 = new TH2D("qmaxMap_1","qmax distribution for MM2; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
TH2D* qmaxMap_2 = new TH2D("qmaxMap_2","qmax distribution for MM1; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
TH2D* qmaxMap_projected_13deg = new TH2D("qmaxMap_projected_13deg","average qmax for projected, cut at 6 degrees; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
const double decaypost_r[] = {-80, 80};
const double decaypost_r[] = {-100, 100};
const int decaypost_nbin = 200;
TH1D* decayPositions_z = new TH1D("decayPositions_z","reconstructed decay positions in z direction; z [mm] ",decaypost_nbin, decaypost_r[0], decaypost_r[1]);
TH1D* decayPositions_z_target = new TH1D("decayPositions_z_target","reconstructed decay positions in z direction for -10<x1,x2,x3<10; z [mm] ",decaypost_nbin, decaypost_r[0], decaypost_r[1]);
TH1D* decayPositions_z_20deg = new TH1D("decayPositions_z_20deg","reconstructed decay positions in z direction, angle<20deg; z [mm] ",decaypost_nbin, decaypost_r[0], decaypost_r[1]);
TH1D* decayPositions_z_17deg = new TH1D("decayPositions_z_17deg","reconstructed decay positions in z direction, angle<17deg; z [mm] ",decaypost_nbin, decaypost_r[0], decaypost_r[1]);
TH1D* decayPositions_z_13deg = new TH1D("decayPositions_z_13deg","reconstructed decay positions in z direction, angle<13deg; z [mm] ",decaypost_nbin, decaypost_r[0], decaypost_r[1]);
TH1D* decayPositions_z_10deg = new TH1D("decayPositions_z_10deg","reconstructed decay positions in z direction, angle<10deg; z [mm] ",decaypost_nbin, decaypost_r[0], decaypost_r[1]);
//with angle cut
std::vector<TH1D*> decayPositions_vec;
for(unsigned int i(0); i < maxAngle_degree.size(); ++i)
{
histname.Form("decayPositions_z_%ideg", maxAngle_degree[i]);
histtitle.Form("reconstructed decay positions in z direction for angle <%ideg; z [mm]", maxAngle_degree[i]);
decayPositions_vec.push_back(new TH1D(histname, histtitle, decaypost_nbin, decaypost_r[0], decaypost_r[1]));
}
TH1D* qMax_All = new TH1D("qMax_All","distribution of (average) qmax; qmax ",200,0,1.4e6);
TH1D* qMax_double = new TH1D("qMax_double","distribution of average qmax (only double hits); qmax ",200,0,1.4e6);
......@@ -251,7 +244,6 @@ void Analyze(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////////////
//// do reconstruction and histogram filling
////////////////////////////////////////////////////////////////////////////////////
// analysis for one detector pair
......@@ -353,21 +345,11 @@ void Analyze(int argc, char *argv[])
if (angle_to_Vertical<maxAngle_rad[l])
{
projectedMaps[l]->Fill(reconstructed_x,reconstructed_y);
if (l==0)
{
projectedZoom_20deg->Fill(reconstructed_x,reconstructed_y);
}
else if (l==1)
{
projectedZoom_17deg->Fill(reconstructed_x,reconstructed_y);
}
else if (l==2)
}
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[l]))
{
qmaxMap_projected_13deg->Fill(reconstructed_x,reconstructed_y,average_qMax);
qMax_13deg->Fill(average_qMax);
qMax_13deg_small->Fill(average_qMax);
decayPositions_vec[l]->Fill(reconstructed_y);
}
}
}
// now prepare z distributions. Cross check: is this the "correct" z dimension?
......@@ -379,90 +361,34 @@ void Analyze(int argc, char *argv[])
//float reconstructed_z = (z_C - z_F)/dist_MM2*(dist_MM2+dist_MM1);
// accepotance of vertical line is largest in center of micromegas - balance this (only very rough approximation)
float acceptance_Factor = 1-TMath::Sqrt(2.)/size_MM*TMath::Abs(reconstructed_y);
decayPositions_z->Fill(reconstructed_y,acceptance_Factor);
if (std::abs(reconstructed_x)<10 && std::abs(hitsVector_MM2[i].x)<10 && std::abs(hitsVector_MM1[i].x)<10) {decayPositions_z_target->Fill(reconstructed_y);}
//timeVSz->Fill(reconstructed_Time,reconstructed_y,acceptance_Factor);
}
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[0]))
{
decayPositions_z_20deg->Fill(reconstructed_y);
}
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[1]))
{
decayPositions_z_17deg->Fill(reconstructed_y);
}
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[2]))
{
decayPositions_z_13deg->Fill(reconstructed_y);
}
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[3]))
{
decayPositions_z_10deg->Fill(reconstructed_y);
}
}
else
{
//std::cout << "no coinc: EventID_C=" << hitsVector_MM1[i].EventID << " EventID_F=" << hitsVector_MM2[i].EventID << std::endl;
}
qMax_All->Fill(average_qMax);
}
std::cout << "fraction of coincidences relative to triggered events:" << count_Coincidences/float(numberOfEvents_MM1) << std::endl;
for (int nBinsX=0;nBinsX<=singleHeatMap_1->GetNbinsX()+1;nBinsX++)
{
for (int nBinsY=0;nBinsY<=singleHeatMap_1->GetNbinsY()+1;nBinsY++)
{
if (singleHeatMap_1->GetBinContent(nBinsX,nBinsY)>0)
{ qmaxMap_1->SetBinContent(nBinsX,nBinsY,qmaxMap_1->GetBinContent(nBinsX,nBinsY)/singleHeatMap_1->GetBinContent(nBinsX,nBinsY));}
if (singleHeatMap_2->GetBinContent(nBinsX,nBinsY)>0)
{ qmaxMap_2->SetBinContent(nBinsX,nBinsY,qmaxMap_2->GetBinContent(nBinsX,nBinsY)/singleHeatMap_2->GetBinContent(nBinsX,nBinsY));}
if (projectedHeatMap_13deg->GetBinContent(nBinsX,nBinsY)>0)
{ qmaxMap_projected_13deg->SetBinContent(nBinsX,nBinsY,qmaxMap_projected_13deg->GetBinContent(nBinsX,nBinsY)/projectedHeatMap_13deg->GetBinContent(nBinsX,nBinsY));}
}
float acceptance_Factor = 1-TMath::Sqrt(2.)/size_MM*TMath::Abs(reconstructed_y);
decayPositions_z->Fill(reconstructed_y,acceptance_Factor);
if (std::abs(reconstructed_x)<10 && std::abs(hitsVector_MM2[i].x)<10 && std::abs(hitsVector_MM1[i].x)<10) {decayPositions_z_target->Fill(reconstructed_y);}
for (int l=0;l<maxAngle_rad.size();l++)
{
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[0]))
{
decayPositions_vec[l]->Fill(reconstructed_y);
}
}
qMax_All->Fill(average_qMax);
}
//std::cout << "fraction of coincidences relative to triggered events:" << count_Coincidences/float(numberOfEvents_MM1) << std::endl;
}
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// plot histograms
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//paramters
TParameter<Int_t>* triggers = new TParameter<Int_t>("Triggers", numberOfEvents_MM1);
TParameter<Int_t>* coincidences = new TParameter<Int_t>("Coincidences", count_Coincidences);
newFile->cd();
newFile->WriteObject(triggers, triggers->GetName());
newFile->WriteObject(coincidences, coincidences->GetName());
newFile->Write();
newFile->Close();
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///// plot the histograms
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//std::cout << " position 8 Plot\n";
//c1->Divide(2,2);
//TCanvas* c2 = new TCanvas;
//c2->Divide(2,2);
//TCanvas* c3 = new TCanvas;
//c3->Divide(2,2);
//c1->cd(1);
//Sum_coinc_1->SetTitle("A1-B1 coincidence");
//Sum_coinc_1->SetStats(kFALSE);
//Sum_coinc_1->SetLineColor(kBlue);
//Sum_coinc_1->DrawCopy();
//BG_coinc_1->SetLineColor(kRed);
//BG_coinc_1->DrawCopy("same");
newFile->Close();
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment