Commit 442d401a authored by positron cool guys's avatar positron cool guys
Browse files

updated scripts

parent 937252d2
GCC=g++
FLAGS=-std=c++11
INCLUDE := include/
DEP := include/Micromega.cc include/MicromegasPlane.cc include/physicalStrip.cc include/stripsCluster.cc
DEP := include/Micromega.cc include/MicromegasPlane.cc include/physicalStrip.cc include/stripsCluster.cc
LIBS := `root-config --cflags --glibs`
all: map_multiplex ConvertOutputMM analyze_PositionData reconstructAndPlot
......
......@@ -181,7 +181,7 @@ int main(int argc, char* argv[])
if (passedCuts) theNtuples[mmIndex]->Fill(iEntry,theMMvector[mmIndex].GetHitPositionX(),theMMvector[mmIndex].GetHitPositionY(),theMMvector[mmIndex].GetHitTimeRise(),theMMvector[mmIndex].GetHitTimeMax(),theMMvector[mmIndex].GetHitChargeMax(),theMMvector[mmIndex].GetHitChargeMin(),theMMvector[mmIndex].GetHitChargeTot(),float(theMMvector[mmIndex].GetHitNumberOfStrips()));
}
}
if (iEntry%100 == 0) {std::cout << "Analyzed " << iEntry << " events of total " << entries << std::endl;}
if (iEntry%10000 == 0) {std::cout << "Analyzed " << iEntry << " events of total " << entries << std::endl;}
}
......
......@@ -2,7 +2,7 @@
#ifdef __MAKECINT__
#pragma link C++ class vector<vector<short> >+;
#endif
//#pragma link C++ class vector<vector<short> >+;
#include "Dict.h"
#include <iostream>
#include <string>
using namespace std;
......@@ -93,7 +93,7 @@ int main(int argc, char *argv[]){
std::string run_name = file.substr(found_run,found_dotroot-found_run);
std::string path = file.substr(0,found_run);
std::string new_file = path+run_name+"_mapped.root";
std::string new_file = path+"/mapped/" + run_name+"_mapped.root";
......@@ -133,14 +133,14 @@ int main(int argc, char *argv[]){
aux->Branch("raw_q" , &write_raw_q_parse );
std::cout << "multiplexed maps: \n";
// loop on multiplexed map file X
int numOfLines_mapx=320;
double first_mapx[numOfLines_mapx], second_mapx[numOfLines_mapx], third_mapx[numOfLines_mapx], fourth_mapx[numOfLines_mapx], fifth_mapx[numOfLines_mapx];
ifstream InputRead_mapx ("multiplexing_maps/Mapping_new_1.txt");
for (int i=0; i<numOfLines_mapx; i++) {
InputRead_mapx >> first_mapx[i] >> second_mapx[i] >> third_mapx[i] >> fourth_mapx[i] >> fifth_mapx[i];
std::cout<< first_mapx[i] <<" "<< second_mapx[i] <<" "<< third_mapx[i] <<" "<< fourth_mapx[i] <<" "<< fifth_mapx[i] <<" "<<std::endl;
//std::cout<< first_mapx[i] <<" "<< asecond_mapx[i] <<" "<< third_mapx[i] <<" "<< fourth_mapx[i] <<" "<< fifth_mapx[i] <<" "<<std::endl;
} // loop on map file
// loop on multiplexed map file X
......@@ -179,7 +179,6 @@ int main(int argc, char *argv[]){
for (int j=0; j<(int)apv_raw_obj->srsFec->size(); j++) {
if(apv_raw_obj->mmChamber->at(j)=="Chamber_1"){
//X remapping
if(apv_raw_obj->mmReadout->at(j)==49){
for (int p=0; p<numOfLines_mapx; p++) {
......@@ -225,28 +224,11 @@ int main(int argc, char *argv[]){
}
}
}
else {
charge_info_to_be_puched_back_temp.clear();
for(int i_sample=0; i_sample<apv_raw_obj->raw_q->at(j).size(); i_sample++)
{
charge_info_to_be_puched_back_temp.push_back(apv_raw_obj->raw_q->at(j).at(i_sample));
}
write_raw_q_parse.push_back(charge_info_to_be_puched_back_temp);
write_mmStrip.push_back(apv_raw_obj->mmStrip->at(j));
write_srsChan.push_back(apv_raw_obj->srsChan->at(j));
write_srsFec.push_back(apv_raw_obj->srsFec->at(j));
write_srsChip.push_back(apv_raw_obj->srsChip->at(j));
write_mmChamber.push_back(apv_raw_obj->mmChamber->at(j));
write_mmLayer.push_back(apv_raw_obj->mmLayer->at(j));
write_mmReadout.push_back(apv_raw_obj->mmReadout->at(j));
}
}
}
aux->Fill();
}
}
// make second tree
......
......@@ -122,17 +122,17 @@ void Analyze(int argc, char *argv[])
MM_Names.push_back("MM1");
MM_Names.push_back("MM2");
// distance between two detectors
float dist_MM12 = 50.; // 66.5;
// 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 = 100.; //35.75;
float dist_MM1 = 75.; //modify for the setup, preliminary
float size_MM = 80.0; // size in mm
// cut angle in degrees
std::vector<float> maxAngle_degree = {10.0,8.0,6.0,4.0};
std::vector<float> maxAngle_degree = {20.0,17.0,13,10.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());}
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";}
std::cout << " position 1\n";
......@@ -159,6 +159,7 @@ void Analyze(int argc, char *argv[])
std::cout << " position 2 after reading in\n";
TFile* newFile = new TFile(TString(outputFileName),"RECREATE");
// prepare a canvas for the plots
TCanvas* c1 = new TCanvas;
......@@ -166,42 +167,43 @@ void Analyze(int argc, char *argv[])
///// initialize histograms
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
std::vector<TH2*> projectedMaps;
std::vector<TH2D*> projectedMaps;
/// concidences detected
TH2* projectedHeatMap_10deg = new TH2D("projectedHeatMap_10deg","reconstructed Map, angular cut 10 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
TH2* projectedHeatMap_8deg = new TH2D("projectedHeatMap_8deg","reconstructed Map, angular cut 8 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
TH2* projectedHeatMap_6deg = new TH2D("projectedHeatMap_6deg","reconstructed Map, angular cut 6 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
TH2* projectedHeatMap_4deg = new TH2D("projectedHeatMap_4deg","reconstructed Map, angular cut 4 degrees; x [mm] ; z [mm] ",160,-43,43,160,-43,43);
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);
projectedMaps.push_back(projectedHeatMap_8deg);
projectedMaps.push_back(projectedHeatMap_6deg);
projectedMaps.push_back(projectedHeatMap_4deg);
TH2* projectedZoom_10deg = new TH2D("projectedZoom_10deg","reconstructed Map, angular cut 10 degrees; x [mm] ; z [mm] ",110,-14,14,125,-11,20);
TH2* projectedZoom_8deg = new TH2D("projectedZoom_8deg","reconstructed Map, angular cut 8 degrees; x [mm] ; z [mm] ",110,-14,14,125,-11,20);
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);
TH2* singleHeatMap_1 = new TH2D("singleHeatMap_1","histogram for MM2; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
TH2* singleHeatMap_2 = new TH2D("singleHeatMap_2","histogram for MM1; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
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);
TH2* qmaxMap_1 = new TH2D("qmaxMap_1","qmax distribution for MM2; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
TH2* qmaxMap_2 = new TH2D("qmaxMap_2","qmax distribution for MM1; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
TH2* qmaxMap_projected_6deg = new TH2D("qmaxMap_projected_6deg","average qmax for projected, cut at 6 degrees; x [mm] ; z [mm] ",80,-43,43,80,-43,43);
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);
TH1* decayPositions_z = new TH1D("decayPositions_z","reconstructed decay positions in z direction; z [mm] ",100,-32,32);
TH1* decayPositions_z_target = new TH1D("decayPositions_z_target","reconstructed decay positions in z direction for -10<x1,x2,x3<10; z [mm] ",100,-32,32);
TH1* decayPositions_z_10deg = new TH1D("decayPositions_z_10deg","reconstructed decay positions in z direction for -10<x<10, angle<10deg; z [mm] ",120,-32,32);
TH1D* decayPositions_z = new TH1D("decayPositions_z","reconstructed decay positions in z direction; z [mm] ",100,-32,32);
TH1D* decayPositions_z_target = new TH1D("decayPositions_z_target","reconstructed decay positions in z direction for -10<x1,x2,x3<10; z [mm] ",100,-32,32);
TH1D* decayPositions_z_20deg = new TH1D("decayPositions_z_20deg","reconstructed decay positions in z direction for -10<x<10, angle<20deg; z [mm] ",120,-32,32);
TH1* qMax_All = new TH1D("qMax_All","distribution of (average) qmax; qmax ",200,0,1.4e6);
TH1* qMax_double = new TH1D("qMax_double","distribution of average qmax (only double hits); qmax ",200,0,1.4e6);
TH1* qMax_6deg = new TH1D("qMax_6deg","average qmax passed 6 degree cuts; qmax",200,0,1.4e6);
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);
TH1D* qMax_13deg = new TH1D("qMax_13deg","average qmax passed 6 degree cuts; qmax",200,0,1.4e6);
TH1* qMax_double_small = new TH1D("qMax_double_small","distribution of average qmax (only double hits); qmax ",200,0,1.e5);
TH1* qMax_6deg_small = new TH1D("qMax_6deg_small","average qmax passed 6 degree cuts; qmax",200,0,1.e5);
TH1D* qMax_double_small = new TH1D("qMax_double_small","distribution of average qmax (only double hits); qmax ",200,0,1.e5);
TH1D* qMax_13deg_small = new TH1D("qMax_13deg_small","average qmax passed 6 degree cuts; qmax",200,0,1.e5);
TH1D* angle = new TH1D("angle","angle to vertical; angle [rad]",1000,0,1.);
//TH2* timeVSz = new TH2D("timeVSz","reconstructed decay positions over time; timeBin; z [mm]",27,-0.5,26.5,40,-20,20);
//TH2D* timeVSz = new TH2D("timeVSz","reconstructed decay positions over time; timeBin; z [mm]",27,-0.5,26.5,40,-20,20);
std::cout << " position 3 defined histograms\n";
......@@ -300,9 +302,10 @@ void Analyze(int argc, char *argv[])
// average for time
int reconstructed_Time = (hitsVector_MM1[i].t+hitsVector_MM2[i].t)/2.;
float angle_to_Vertical = TMath::ATan( TMath::Sqrt(TMath::Power(Delta_x,2)+TMath::Power(Delta_y,2)) / dist_MM12);
std::cout << angle_to_Vertical << "\n";
//std::cout << angle_to_Vertical << "\n";
angle->Fill(angle_to_Vertical);
// angular cut for 2D histo
for (int l=0;l<maxAngle_rad.size();l++)
......@@ -312,17 +315,17 @@ void Analyze(int argc, char *argv[])
projectedMaps[l]->Fill(reconstructed_x,reconstructed_y);
if (l==0)
{
projectedZoom_10deg->Fill(reconstructed_x,reconstructed_y);
projectedZoom_20deg->Fill(reconstructed_x,reconstructed_y);
}
else if (l==1)
{
projectedZoom_8deg->Fill(reconstructed_x,reconstructed_y);
projectedZoom_17deg->Fill(reconstructed_x,reconstructed_y);
}
else if (l==2)
{
qmaxMap_projected_6deg->Fill(reconstructed_x,reconstructed_y,average_qMax);
qMax_6deg->Fill(average_qMax);
qMax_6deg_small->Fill(average_qMax);
qmaxMap_projected_13deg->Fill(reconstructed_x,reconstructed_y,average_qMax);
qMax_13deg->Fill(average_qMax);
qMax_13deg_small->Fill(average_qMax);
}
}
}
......@@ -344,7 +347,7 @@ void Analyze(int argc, char *argv[])
}
if (TMath::Abs(Delta_y)/dist_MM12 < TMath::Tan(maxAngle_rad[0]))
{
if (std::abs(reconstructed_x)<10) decayPositions_z_10deg->Fill(reconstructed_y);
if (std::abs(reconstructed_x)<10) decayPositions_z_20deg->Fill(reconstructed_y);
}
}
else
......@@ -367,8 +370,8 @@ void Analyze(int argc, char *argv[])
if (singleHeatMap_2->GetBinContent(nBinsX,nBinsY)>0)
{ qmaxMap_2->SetBinContent(nBinsX,nBinsY,qmaxMap_2->GetBinContent(nBinsX,nBinsY)/singleHeatMap_2->GetBinContent(nBinsX,nBinsY));}
if (projectedHeatMap_6deg->GetBinContent(nBinsX,nBinsY)>0)
{ qmaxMap_projected_6deg->SetBinContent(nBinsX,nBinsY,qmaxMap_projected_6deg->GetBinContent(nBinsX,nBinsY)/projectedHeatMap_6deg->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));}
}
}
......@@ -378,43 +381,10 @@ void Analyze(int argc, char *argv[])
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// plot histograms
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
std::cout << " position 7 write to file\n";
TFile* newFile = new TFile(TString(outputFileName),"RECREATE");
newFile->cd();
gDirectory->WriteObject(projectedHeatMap_10deg,"projectedHeatMap_10deg");
gDirectory->WriteObject(projectedHeatMap_8deg,"projectedHeatMap_8deg");
gDirectory->WriteObject(projectedHeatMap_6deg,"projectedHeatMap_6deg");
gDirectory->WriteObject(projectedHeatMap_4deg,"projectedHeatMap_4deg");
gDirectory->WriteObject(projectedZoom_10deg,"projectedZoom_10deg");
gDirectory->WriteObject(projectedZoom_8deg,"projectedZoom_8deg");
gDirectory->WriteObject(qmaxMap_1,"qmaxMap_1");
gDirectory->WriteObject(qmaxMap_2,"qmaxMap_2");
gDirectory->WriteObject(qmaxMap_projected_6deg,"qmaxMap_projected_6deg");
gDirectory->WriteObject(decayPositions_z,"decayPositions_z");
gDirectory->WriteObject(decayPositions_z_target,"decayPositions_z_target");
gDirectory->WriteObject(decayPositions_z_10deg,"decayPositions_z_10deg");
gDirectory->WriteObject(qMax_All,"qMax_All");
gDirectory->WriteObject(qMax_double,"qMax_double");
gDirectory->WriteObject(qMax_6deg,"qMax_6deg");
gDirectory->WriteObject(qMax_double_small,"qMax_double_small");
gDirectory->WriteObject(qMax_6deg_small,"qMax_6deg_small");
gDirectory->WriteObject(singleHeatMap_1,"singleHeatMap_1");
gDirectory->WriteObject(singleHeatMap_2,"singleHeatMap_2");
newFile->Write();
newFile->Close();
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///// plot the histograms
......
Supports Markdown
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