Commit 3cbfcd07 authored by Adrien Matta's avatar Adrien Matta
Browse files

* Progress on BDC analysis

parent 8ab31976
......@@ -47,7 +47,7 @@ ClassImp(TSamuraiBDCPhysics)
m_EventPhysics = this ;
//m_Spectra = NULL;
ToTThreshold_L = 0;
ToTThreshold_H = 180;
ToTThreshold_H = 1000;
DriftLowThreshold=0.1 ;
DriftUpThreshold=2.4;
PowerThreshold=5;
......@@ -63,162 +63,164 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
PreTreat();
static unsigned int det,layer,wire,size;
for(auto it = m_DCHit.begin(); it!=m_DCHit.end(); it++){
det = it->first;
size = it->second.size();
for (auto itt = it->second.begin(); itt != it->second.end(); itt++){
double angle = itt->first;
// cout << det << " " << angle << " " << itt->second.size() << endl;
}
}
/*
// RemoveNoise();
static double dl;
// Map[plane angle, vector of spatial information]
static map<double, vector<double> > X ;
static map<double, vector<double> > Z ;
static map<double, vector<double> > R ;
static int det,layer,wire;
X.clear();Z.clear();R.clear();
unsigned int size = Detector.size();
for(unsigned int i = 0 ; i < size ; i++){
if(DriftLength[i] > DriftLowThreshold && DriftLength[i] < DriftUpThreshold){
det = Detector[i];
layer = Layer[i];
wire = Wire[i];
SamuraiDCIndex idx(det,layer,wire);
X[Wire_Angle[idx]].push_back(Wire_X[idx]);
Z[Wire_Angle[idx]].push_back(Wire_Z[idx]);
R[Wire_Angle[idx]].push_back(DriftLength[i]);
}
}
// Reconstruct the vector for each of the plane of each of the detector
static double X0,X100,a,b; // store the BuildTrack2D results
static map<double, TVector3 > VX0 ;
static map<double, TVector3 > VX100 ;
static map<double, double > D ;// the minimum distance
static unsigned int uid; uid=0;
VX0.clear();VX100.clear(),D.clear();
static vector<TVector3> C ;
static vector<double > W ; // weight based on D
unsigned int count = 0 ;
for(auto it = m_DCHit.begin(); it!=m_DCHit.end(); it++){
// Each entry in the map is a detector
det = it->first;
Detector.push_back(det);
PosX.push_back(0);
PosY.push_back(0);
ThetaX.push_back(0);
PhiY.push_back(0);
devX.push_back(0);
devY.push_back(0);
Dir.push_back(TVector3());
PileUp.push_back(0);
X.clear();Z.clear();R.clear();
// Build the necessary X,Z,R vector
for(auto itt = it->second.begin() ; itt!= it->second.end() ; itt++){
dl = (*itt).DriftLength;
if(dl > DriftLowThreshold && dl < DriftUpThreshold){
SamuraiDCIndex idx(det,(*itt).Layer,(*itt).Wire);
X[Wire_Angle[idx]].push_back(Wire_X[idx]);
Z[Wire_Angle[idx]].push_back(Wire_Z[idx]);
R[Wire_Angle[idx]].push_back(dl);
}
}
// Reconstruct the vector for each of the plane of each of the detector
VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
uid=0;
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]);
m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]);
#else
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
#endif
}
#if __cplusplus > 199711L && NPMULTITHREADING
// do all plan at once in parallele, return when all plan are done
m_reconstruction.BuildTrack2D();
uid=0;
#endif
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b);
// do all plan at once in parallele, return when all plan are done
m_reconstruction.BuildTrack2D();
uid=0;
#endif
// Loop over the results
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b);
#endif
// for Debug, write a file of
// { std::ofstream f("distance.txt", std::ios::app);
// f<< D[it->first] << endl;
// f.close();
// }
// very large "a" means track perpendicular to the chamber, what happen when there is pile up
if(abs(a)>5000)
PileUp[count]++;
// Position at z=0
TVector3 P(X0,0,0);
P.RotateZ(it->first);
VX0[it->first]=P;
// Position at z=100
TVector3 P100= TVector3(X100,0,0);
P100.RotateZ(it->first);
VX100[it->first]=P100;
// Reconstruct the central position (z=0) for each detector
C.clear(),W.clear();
for(auto it1 = VX0.begin();it1!=VX0.end();++it1){
for(auto it2 = it1;it2!=VX0.end();++it2){
if(it1!=it2){// different plane
//cout << "BDC" << endl;
m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
// cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
C.push_back(P);
// Mean pos are weighted based on the the sum of distance from track
// to hit obtained during the minimisation
W.push_back(1./sqrt(D[it1->first]*D[it2->first]));
}
}
}
}
// for Debug, write a file of
// { std::ofstream f("distance.txt", std::ios::app);
// f<< D[it->first] << endl;
// f.close();
// }
// very large a means track perpendicular to the chamber, what happen when there is pile up
if(abs(a)>5000)
PileUp++;
Mult+=X[it->first].size();
// Position at z=0
TVector3 P(X0,0,0);
P.RotateZ(it->first);
VX0[it->first]=P;
// Position at z=100
TVector3 P100= TVector3(X100,0,0);
P100.RotateZ(it->first);
VX100[it->first]=P100;
// Reconstruct the position at z=100 for each detector
static vector<TVector3> C100 ;
C100.clear();
for(auto it1 = VX100.begin();it1!=VX100.end();++it1){
for(auto it2 = it1;it2!=VX100.end();++it2){
if(it1!=it2 ){// different plane, same detector
m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
}
// Reconstruct the central position (z=0) for each detector
static vector<TVector3> C ;
static vector<double > W ; // weight based on D
C.clear(),W.clear();
TVector3 P;
for(auto it1 = VX0.begin();it1!=VX0.end();++it1){
for(auto it2 = it1;it2!=VX0.end();++it2){
if(it1!=it2){// different plane
//cout << "BDC" << endl;
m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
// cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
C.push_back(P);
// Mean pos are weighted based on the the sum of distance from track
// to hit obtained during the minimisation
W.push_back(1./sqrt(D[it1->first]*D[it2->first]));
if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
C100.push_back(P);
}
}
}
}
}
// Reconstruct the position at z=100 for each detector
static vector<TVector3> C100 ;
C100.clear();
for(auto it1 = VX100.begin();it1!=VX100.end();++it1){
for(auto it2 = it1;it2!=VX100.end();++it2){
if(it1!=it2 ){// different plane, same detector
m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
C100.push_back(P);
// Build the Reference position by averaging all possible pair
size = C.size();
static double PosX100,PosY100,norm;
if(size){
norm=0;
for(unsigned int i = 0 ; i < size ; i++){
PosX[count]+= C[i].X()*W[i];
PosY[count]+= C[i].Y()*W[i];
PosX100+= C100[i].X()*W[i];
PosY100+= C100[i].Y()*W[i];
norm+=W[i];
}
//MultMean=size;
// Mean position at Z=0
PosX[count]/=norm;
PosY[count]/=norm;
// Mean position at Z=100
PosX100/=norm;
PosY100/=norm;
for(unsigned int i = 0 ; i < size ; i++){
devX[count]+=W[i]*(C[i].X()-PosX[count])*(C[i].X()-PosX[count]);
devY[count]+=W[i]*(C[i].Y()-PosY[count])*(C[i].Y()-PosY[count]);
}
devX[count]=sqrt(devX[count]/((size-1)*norm));
devY[count]=sqrt(devY[count]/((size-1)*norm));
// Compute ThetaX, angle between the Direction vector projection in XZ with
// the Z axis
//ThetaX=atan((PosX100-PosX)/100.);
ThetaX[count] = (PosX100-PosX[count])/100.;
// Compute PhiY, angle between the Direction vector projection in YZ with
// the Z axis
//PhiY=atan((PosY100-PosY)/100.);
PhiY[count]=(PosY100-PosY[count])/100.;
Dir[count]=TVector3(PosX100-PosX[count],PosY100-PosY[count],100).Unit();
}
}
}
// Build the Reference position by averaging all possible pair
size = C.size();
static double PosX100,PosY100,norm;
if(size){
PosX=0;
PosY=0;
PosX100=0;
PosY100=0;
norm=0;
for(unsigned int i = 0 ; i < size ; i++){
PosX+= C[i].X()*W[i];
PosY+= C[i].Y()*W[i];
PosX100+= C100[i].X()*W[i];
PosY100+= C100[i].Y()*W[i];
norm+=W[i];
}
MultMean=size;
// Mean position at Z=0
PosX=PosX/norm;
PosY=PosY/norm;
// Mean position at Z=100
PosX100=PosX100/norm;
PosY100=PosY100/norm;
devX=0;
devY=0;
for(unsigned int i = 0 ; i < size ; i++){
devX+=W[i]*(C[i].X()-PosX)*(C[i].X()-PosX);
devY+=W[i]*(C[i].Y()-PosY)*(C[i].Y()-PosY);
if(PosX[count]==0){
PosX[count]=-10000;
PosY[count]=-10000;
ThetaX[count]=-10000;
PhiY[count]=-10000;
}
devX=sqrt(devX/((size-1)*norm));
devY=sqrt(devY/((size-1)*norm));
// Compute ThetaX, angle between the Direction vector projection in XZ with
// the Z axis
//ThetaX=atan((PosX100-PosX)/100.);
ThetaX = (PosX100-PosX)/100.;
// Compute PhiY, angle between the Direction vector projection in YZ with
// the Z axis
//PhiY=atan((PosY100-PosY)/100.);
PhiY=(PosY100-PosY)/100.;
Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit();
}
*/
count++;
}// detector loop
return;
}
///////////////////////////////////////////////////////////////////////////
......@@ -257,8 +259,7 @@ void TSamuraiBDCPhysics::PreTreat(){
if(etime && time && etime-time>ToTThreshold_L && etime-time<ToTThreshold_H){
channel="SamuraiBDC"+NPL::itoa(det)+"/L" + NPL::itoa(layer);
SamuraiDCIndex idx(det,layer,wire);
double angle = Wire_Angle[idx];
m_DCHit[det][angle].push_back(DCHit(layer,wire,time,etime-time,2.5-Cal->ApplySigmoid(channel,etime)));
m_DCHit[det].push_back(DCHit(det,layer,wire,time,etime-time,2.5-Cal->ApplySigmoid(channel,etime)));
}
}
......@@ -268,25 +269,15 @@ void TSamuraiBDCPhysics::PreTreat(){
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::Clear(){
m_DCHit.clear();
/* DriftLength.clear();
Detector.clear();
Layer.clear();
Wire.clear();
Time.clear();
ToT.clear();
*/
// Computed variable
ParticleDirection.clear();
MiddlePosition.clear();
PosX.clear();
PosY.clear();
ThetaX.clear();
PhiY.clear();
devX.clear();
devY.clear();
Dir.clear();
PileUp.clear();
// Computed variable
PosX.clear();
PosY.clear();
ThetaX.clear();
PhiY.clear();
devX.clear();
devY.clear();
Dir.clear();
PileUp.clear();
}
///////////////////////////////////////////////////////////////////////////
......@@ -320,33 +311,33 @@ void TSamuraiBDCPhysics::ReadConfiguration(NPL::InputParser parser){
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::AddDC(int det, NPL::XmlParser& xml){
std::string name = "SAMURAIBC"+NPL::itoa(det);
std::string name = "SAMURAIBDC"+NPL::itoa(det);
std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName(name);
unsigned int sizeB = b.size();
for(unsigned int i = 0 ; i < sizeB ; i++){
unsigned int layer = b[i]->AsInt("layer");
unsigned int wire = b[i]->AsInt("wireid");
double X = b[i]->AsDouble("wirepos");
double Z = b[i]->AsDouble("wirez");
string sDir = b[i]->AsString("anodedir");
double T=0;
if(sDir=="X")
T= 0*deg;
else if(sDir=="Y")
T= -90*deg;
else if(sDir=="U")
T=-30*deg;
else if(sDir=="V")
T=+30*deg;
else{
cout << "ERROR: Unknown layer orientation for Samurai BDC"<< endl;
exit(1);
}
SamuraiDCIndex idx(det,layer,wire);
Wire_X[idx]=X;
Wire_Z[idx]=Z;
Wire_Angle[idx]=T;
unsigned int sizeB = b.size();
for(unsigned int i = 0 ; i < sizeB ; i++){
unsigned int layer = b[i]->AsInt("layer");
unsigned int wire = b[i]->AsInt("wireid");
double X = b[i]->AsDouble("wirepos");
double Z = b[i]->AsDouble("wirez");
string sDir = b[i]->AsString("anodedir");
double T=0;
if(sDir=="X")
T= 0*deg;
else if(sDir=="Y")
T= -90*deg;
else if(sDir=="U")
T=-30*deg;
else if(sDir=="V")
T=+30*deg;
else{
cout << "ERROR: Unknown layer orientation for Samurai BDC"<< endl;
exit(1);
}
SamuraiDCIndex idx(det,layer,wire);
Wire_X[idx]=X;
Wire_Z[idx]=Z;
Wire_Angle[idx]=T;
}
}
///////////////////////////////////////////////////////////////////////////
......@@ -391,10 +382,10 @@ void TSamuraiBDCPhysics::AddParameterToCalibrationManager(){
// for each det
for( int d = 1 ; d < 3 ; ++d){
// each layer
for( int l = 0 ; l < 8 ; ++l){
Cal->AddParameter("SamuraiBDC"+NPL::itoa(d), "L"+ NPL::itoa(l),"BDC"+NPL::itoa(d)+"_L"+ NPL::itoa(l));
}
// each layer
for( int l = 0 ; l < 8 ; ++l){
Cal->AddParameter("SamuraiBDC"+NPL::itoa(d), "L"+ NPL::itoa(l),"BDC"+NPL::itoa(d)+"_L"+ NPL::itoa(l));
}
}
}
......
......@@ -50,8 +50,8 @@
class DCHit{
public:
DCHit(unsigned int layer , unsigned int wire, double time, double tot, double drift){
DCHit(unsigned int detector, unsigned int layer , unsigned int wire, double time, double tot, double drift){
Detector=detector;
Layer=layer;
Wire=wire;
Time=time;
......@@ -61,7 +61,8 @@ class DCHit{
DCHit(){};
~DCHit(){};
private:
public:
unsigned int Detector;
unsigned int Layer;
unsigned int Wire;
double Time;
......@@ -80,13 +81,11 @@ class TSamuraiBDCPhysics : public TObject, public NPL::VDetector{
void Clear(const Option_t*) {};
public:
// double map of [bdc number][layer angle]=vector of hit
std::map<unsigned int,std::map<double,std::vector<DCHit>>> m_DCHit;
// map of [bdc number]=vector of hit
std::map<unsigned int, std::vector<DCHit>> m_DCHit;//!
// Computed variable
std::map<unsigned int, std::vector<TVector3> > ParticleDirection;
std::map<unsigned int, std::vector<TVector3> > MiddlePosition;
// Computed variable
std::vector<double> PosX;
std::vector<double> PosY;
std::vector<double> ThetaX;
......
BDC1_L0 2.5 0.0787697 1800.25
BDC1_L1 2.5 0.0781409 1800.37
BDC1_L2 2.5 0.0773193 1800.64
BDC1_L3 2.5 0.0780069 1800.26
BDC1_L4 2.5 0.0795404 1798.41
BDC1_L5 2.5 0.0792411 1798.26
BDC1_L6 2.5 0.0791663 1798.37
BDC1_L7 2.5 0.0789406 1798.97
BDC2_L0 2.5 0.0780368 1793.31
BDC2_L1 2.5 0.0774875 1793.08
BDC2_L2 2.5 0.0774375 1792.99
BDC2_L3 2.5 0.0782832 1791.84
BDC2_L4 2.5 0.0766781 1792.56
BDC2_L5 2.5 0.0786876 1790.96
BDC2_L6 2.5 0.078485 1790.88
BDC2_L7 2.5 0.0779132 1792.36
void CalibBDC(){
for(int run12 =824; run12 <= 824 ; run12++){
TChain* c = new TChain("RawTree");
c->Add(Form("../../root/mrdc/run%d/run%d_SAMURAIBDC.root",run12,run12));
c->SetBranchStatus("*",false);
// BDC2
int BDC=1;
int NLayer = 8;
c->SetBranchStatus("SamuraiBDC",true);
ofstream fout(Form("BDC%d_Calib_%d.txt",BDC, run12));
double inter = 2.5;// inter wire distance / 2
double low=1600;
double high=1900;
cout <<"Calibrating BDC" << BDC << " run number : " << run12 << endl;
string cond;
TF1* f = new TF1("sigmoid","[0]/(exp([1]*([2]-x))+1)");
//TF1* f = new TF1("sigmoid","0*(x<[0])+(x-[0])*[1]*(x>[0]&&x<[2])+[3]*(x>[2])");
TH1D* h = new TH1D("h","h",1500,0,3000);
for(unsigned int i = 0 ; i < NLayer; i++){
cond = Form("fBDC_LayerNbr==%d&&fBDC_Edge==0&&fBDC_Time>%f&&fBDC_Time<%f&&fBDC_Time@.size()>10&&fBDC_Time@.size()<100 &&fBDC_DetectorNbr==%d ",i,low,high,BDC);
c->Draw("fBDC_Time>>h",cond.c_str());
TH1D* g = new TH1D(*h);
unsigned int size = h->GetNbinsX();
for(unsigned int i = 0 ; i < size ; i++){
g->SetBinContent(i,h->Integral(0,i));
}
//new TCanvas();
//g->Draw();
cout << g->GetMaximum() << endl;
f->SetParameter(0,g->GetMaximum());
f->SetParameter(1,0.1);
f->SetParameter(2,1700);
f->SetParLimits(1,0,1);
f->SetParLimits(2,1600,1900);
g->Fit(f);
// Renormalize the distribution to the maximum driftlength
f->SetParameter(0,inter);
double p0 = f->GetParameter(0);
double p1 = f->GetParameter(1);
double p2 = f->GetParameter(2);
//new TCanvas();
//c->Draw(Form("%f/(exp(%f*(%f-(Time+ToT)))+1)>>hh(200,0,2.5)",p0,p1,p2),cond.c_str());
fout << Form("BDC%d_L%d ",BDC,i) << p0 << " " << p1 << " " << p2 << endl;
// gPad->Update();
// gPad->WaitPrimitive();
}
fout.close();
/* ////////////////////////////////////////// */
/* // BDC2 */
/* ////////////////////////////////////////// */
/* BDC=2; */
/* NLayer = 14; */
/* c->SetBranchStatus("SamuraiBDC2",true); */
/* ofstream fout2(Form("CalibFiles/BDC2_Calib_%d.txt", run12)); */
/* inter = 10;// inter wire distance */
/* low=1200; */
/* high=1700; */
/* h->Clear(); */
/* for(unsigned int i = 0 ; i < NLayer; i++){ */
/* cond = Form("fBDC%d_LayerNbr==%d&&fBDC%d_Edge==0&&fBDC%d_Time>%f&&fBDC%d_Time<%f&&fBDC%d_Time@.size()>10&&fBDC%d_Time@.size()<100",BDC,i,BDC,BDC,low,BDC,high,BDC,BDC); */
/* c->Draw(Form("fBDC%d_Time>>h",BDC),cond.c_str()); */
/* TH1D* g = new TH1D(*h); */
/* unsigned int size = h->GetNbinsX(); */
/* for(unsigned int i = 0 ; i < size ; i++){ */
/* g->SetBinContent(i,h->Integral(0,i)); */
/* } */
/* //new TCanvas(); */
/* //g->Draw(); */
/* cout << g->GetMaximum() << endl; */
/* f->SetParameter(0,g->GetMaximum()); */
/* f->SetParameter(1,8e-3); */
/* f->SetParameter(2,1500); */
/* g->Fit(f); */
/* // Renormalize the distribution to 10 mm */
/* f->SetParameter(0,inter); */
/* double p0 = f->GetParameter(0); */
/* double p1 = f->GetParameter(1); */
/* double p2 = f->GetParameter(2); */
/* //new TCanvas(); */
/* //c->Draw(Form("%f/(exp(%f*(%f-(Time+ToT)))+1)>>hh(200,0,2.5)",p0,p1,p2),cond.c_str()); */
/* fout2 << Form("BDC%d_L%d ",BDC,i) << p0 << " " << p1 << " " << p2 << endl; */
/* } */
/* fout2.close(); */
}
}
CalibrationFilePath
Calibration/BDC1_Calib.txt
Calibration/BDC/BDC1_Calib_824.txt
Calibration/BDC/BDC2_Calib_824.txt
Calibration/FDC2_Calib.txt
Calibration/FDC0_Calib.txt
Calibration/CalibVDrift/CalibVDrift.txt
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