From 017cd1188f902d5fa79e5d817bbd077f26aa8a82 Mon Sep 17 00:00:00 2001 From: Adrien Matta <matta@lpccaen.in2p3.fr> Date: Wed, 2 Jun 2021 16:21:24 +0200 Subject: [PATCH] * Progress on s034 analysis - BDC alignment on MINOS/SAMURAI --- NPLib/Detectors/Samurai/SamuraiFieldMap.cxx | 55 ++++++------ Projects/S034/Analysis.cxx | 43 ++++++---- Projects/S034/Calibration/Pos/FDC.cxx | 94 +++++++++++++-------- Projects/S034/macro/rigz.cxx | 31 ++++--- Projects/S034/s034.detector | 18 ++-- 5 files changed, 140 insertions(+), 101 deletions(-) diff --git a/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx b/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx index 4498dcdd3..ccb5f57c1 100644 --- a/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx +++ b/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx @@ -58,7 +58,6 @@ double SamuraiFieldMap::Delta(const double* parameter){ //////////////////////////////////////////////////////////////////////////////// double SamuraiFieldMap::FindBrho(TVector3 p_fdc0,TVector3 d_fdc0,TVector3 p_fdc2,TVector3 d_fdc2){ - m_FitPosFDC0=p_fdc0; m_FitDirFDC0=d_fdc0; m_FitPosFDC2=p_fdc2; @@ -123,7 +122,8 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV N.SetBrho(Brho); // track result - std::vector< TVector3 > track; + static std::vector< TVector3 > track; + track.clear(); // starting point of the track if(store){ @@ -132,11 +132,13 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV pos.RotateY(m_angle); } dir=dir.Unit(); - double r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z()); + static double r; + r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z()); // number of step taken - unsigned int count = 0; + static unsigned int count,limit; + count = 0; // maximum number of state before giving up - unsigned int limit = 1000; + limit = 1000; // First propagate to r_max with one line while(r>m_Rmax && count<limit){ @@ -156,17 +158,18 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV return track; } - TVector3 xk1,xk2,xk3,xk4; // position - TVector3 pk1,pk2,pk3,pk4; // impulsion - - double K = N.GetEnergy(); // kinetic energy - double m = N.Mass(); // mc2 - double P = sqrt(K*K+2*K*m)/c_light; // P - double px = P*dir.X();//px - double py = P*dir.Y();//py - double pz = P*dir.Z();//pz - TVector3 imp = P*dir; - double h = 1*nanosecond; + static TVector3 xk1,xk2,xk3,xk4; // position + static TVector3 pk1,pk2,pk3,pk4; // impulsion + static TVector3 imp; + static double K,m,P,px,py,pz; + K = N.GetEnergy(); // kinetic energy + m = N.Mass(); // mc2 + P = sqrt(K*K+2*K*m)/c_light; // P + px = P*dir.X();//px + py = P*dir.Y();//py + pz = P*dir.Z();//pz + imp = P*dir; + static double h = 1*nanosecond; while(r<=m_Rmax && count < limit){ func(N, pos , imp , xk1, pk1); func(N, pos+xk1*(h/2.), imp+pk1*(h/2.) , xk2, pk2); @@ -193,26 +196,26 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV //////////////////////////////////////////////////////////////////////////////// void SamuraiFieldMap::func(NPL::Particle& N, TVector3 pos, TVector3 imp, TVector3& new_pos, TVector3& new_imp){ - double px,py,pz; + static double px,py,pz,vx,vy,vz,Bx,By,Bz,q,P2,D,m2c4; + static vector<double> B; px=imp.X(); py=imp.Y(); pz=imp.Z(); - double P2,D,m2c4; P2=imp.Mag2(); // P2 m2c4 = N.Mass()*N.Mass(); D=sqrt(m2c4+P2*c_squared); // sqrt(m2c4+P2c2) - double vx=px*c_squared/D;// pxc * c / D = pxc2/D - double vy=py*c_squared/D; - double vz=pz*c_squared/D; + vx=px*c_squared/D;// pxc * c / D = pxc2/D + vy=py*c_squared/D; + vz=pz*c_squared/D; new_pos.SetX(vx); new_pos.SetY(vy); new_pos.SetZ(vz); - vector<double> B = InterpolateB(pos); - double Bx= B[0]; - double By= B[1]; - double Bz= B[2]; - double q = N.GetZ()*eplus; // issue with the tesla/coulomb definition + B = InterpolateB(pos); + Bx= B[0]; + By= B[1]; + Bz= B[2]; + q = N.GetZ()*eplus; // issue with the tesla/coulomb definition new_imp.SetX(q*(vy*Bz-vz*By));// q*pyc2*Bz/D -q*pzc2*By/D new_imp.SetY(q*(vz*Bx-vx*Bz)); new_imp.SetZ(q*(vx*By-vy*Bx)); diff --git a/Projects/S034/Analysis.cxx b/Projects/S034/Analysis.cxx index a78dd53f8..544f76824 100644 --- a/Projects/S034/Analysis.cxx +++ b/Projects/S034/Analysis.cxx @@ -49,7 +49,7 @@ void Analysis::Init(){ InitOutputBranch(); InitInputBranch(); // for fdc/bdc alignement - file.open("Calibration/Pos/fdc.txt"); + //file.open("Calibration/Pos/bdc.txt"); } //////////////////////////////////////////////////////////////////////////////// @@ -64,13 +64,9 @@ void Analysis::TreatEvent(){ && FDC2->PosY>-500 && FDC2->PosY<500 && FDC0->PosX>-80 && FDC0->PosX<80 && FDC0->PosY>-80 && FDC0->PosY<80 // both FDC ok - && (Minos->Tracks_P0.size()==1)|| (Minos->Tracks_P0.size()==2)) { // p,pn or p,2p + && (Minos->Tracks_P0.size()>1)) { // p,pn or p,2p // Compute ThetaX and PhiY using Minos vertex and FDC0 X - double FDC0_ThetaX = FDC0->ThetaX; - double FDC0_PhiY = FDC0->PhiY; - TVector3 FDC0_Dir = FDC0->Dir; // Check if both BDC are reconstructed - // TVector3 BDC1=BDC->GetPos(1); TVector3 BDC2=BDC->GetPos(2); @@ -81,7 +77,7 @@ void Analysis::TreatEvent(){ MinimumDistanceTwoLines(BDC1,BDC2, Minos->Tracks_P0[0], P1, Vertex, delta) ; - FDC0_Dir= FDC0->GetPos()-Vertex; + TVector3 FDC0_Dir= FDC0->GetPos()-Vertex; FDC0_Dir=FDC0_Dir.Unit(); TVector3 BDCDir=BDC2-BDC1; BDCDir=BDCDir.Unit(); @@ -93,23 +89,36 @@ void Analysis::TreatEvent(){ Z=Vertex.Z()+4657.39; //double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX}; - if(FDC0->GetPos().X()>-10000 && FDC0_Dir.Z()>0.8){ - FDC0_ThetaX = atan((FDC0->PosX-Vertex.X())/(1283.7-Z)); - FDC0_PhiY = atan((FDC0->PosY-Vertex.Y())/(1283.7-Z)); - double brho_param[6]={FDC0->PosX-1.77, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX+252.55, FDC2->ThetaX}; + if(FDC0_Dir.Z()>0.6){ + double FDC0_ThetaX = atan((FDC0->PosX-Vertex.X())/(1254.39-Z)); + double FDC0_PhiY = atan((FDC0->PosY-Vertex.Y())/(1254.39-Z)); + double brho_param[6]={FDC0->PosX, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX+252.416, FDC2->ThetaX}; BrhoP=r_fit(brho_param); Brho=m_field.FindBrho(FDC0->GetPos(),FDC0_Dir,FDC2->GetPos(),TVector3(0,0,1)); // cout << Brho-BrhoP << endl; } - /* // Calib////////////////////////////////////////////////////////////////// - static int count=0; - if(FDC2->PosX-252.55>0&&FDC0->GetPos().X()>-10000 && FDC0->Dir.Z()>0.8){ - file << FDC0->GetPos().X() <<" " << FDC0->GetPos().Y() << " " << FDC0->GetPos().Z() <<" " << FDC0_Dir.X() <<" " << FDC0_Dir.Y() << " " << FDC0_Dir.Z()<< " " ; - file << FDC2->GetPos().X() <<" " << FDC2->GetPos().Y() << " " << FDC2->GetPos().Z() <<" " << FDC2->Dir.X() <<" " << FDC2->Dir.Y() << " " << FDC2->Dir.Z()<< endl; + // Calib////////////////////////////////////////////////////////////////// + /* static int count=0; + if(Minos->Delta_Vertex < 5 && FDC2->PosX-252.55>0&&FDC0->GetPos().X()>-10000 && FDC0->Dir.Z()>0.9 && Minos->Z_Vertex>0&& sqrt(Minos->X_Vertex*Minos->X_Vertex+Minos->Y_Vertex*Minos->Y_Vertex)<15){ + file << FDC0->GetPos().X() <<" " << FDC0->GetPos().Y() << " " << FDC0->GetPos().Z() <<" " ; + file << Minos->X_Vertex <<" " << Minos->Y_Vertex << " " << Minos->Z_Vertex << " " ; + file << FDC2->GetPos().X() <<" " << FDC2->GetPos().Y() << " " << FDC2->GetPos().Z() <<" " << FDC2->Dir.X() <<" " << FDC2->Dir.Y() << " " << FDC2->Dir.Z()<< endl; + count ++; + } + if(count>1000) + exit(1); + */ + /* static int count=0; + if(Minos->Delta_Vertex < 5 && sqrt(Minos->X_Vertex*Minos->X_Vertex+Minos->Y_Vertex*Minos->Y_Vertex)<15 && Minos->Z_Vertex>-4650){ + file << BDC1.X() <<" " << BDC1.Y()<< " " << BDC1.Z() <<" " ; + file << BDC2.X() <<" " << BDC2.Y()<< " " << BDC2.Z() <<" " ; + file << Minos->X_Vertex <<" " << Minos->Y_Vertex << " " << Minos->Z_Vertex << endl ; count ++; } if(count>10000) - exit(1);*/ + exit(1); + */ + } } } diff --git a/Projects/S034/Calibration/Pos/FDC.cxx b/Projects/S034/Calibration/Pos/FDC.cxx index 35131be2c..0bf889b1c 100644 --- a/Projects/S034/Calibration/Pos/FDC.cxx +++ b/Projects/S034/Calibration/Pos/FDC.cxx @@ -1,7 +1,7 @@ vector<TVector3*> pos0; -vector<TVector3*> dir0; +vector<TVector3*> posM; vector<TVector3*> pos2; vector<TVector3*> dir2; SamuraiFieldMap field; @@ -10,23 +10,28 @@ SamuraiFieldMap field; double devB(const double* parameter){ // Return the standard deviation in Brho unsigned int size = pos0.size(); - TVector3 o0(parameter[0],parameter[1],parameter[2]); - TVector3 o2(parameter[3],parameter[4],parameter[5]); - TVector3 p0,p2; - field.SetFDC2R(parameter[5]); - double Brho; + TVector3 oM(parameter[1],parameter[2],parameter[3]); + TVector3 o0(parameter[4],parameter[5],parameter[6]); + TVector3 o2(parameter[7],parameter[8],parameter[9]); + TVector3 pM,p0,d0,p2; + oM.Print(); + o0.Print(); + o2.Print(); + field.SetFDC2R(parameter[6]); static auto h = new TH1D("h","h", 1000,0,10); h->Reset(); - for(unsigned int i = 0 ; i < 1000 ; i++){ + for(unsigned int i = 0 ; i < size ; i++){ + pM=*(posM[i]); + pM.RotateZ(parameter[0]*M_PI/180.); + pM+=oM; p0=*(pos0[i])+o0; p2=*(pos2[i])+o2; - - if(dir0[i]->Z()>0.8) - h->Fill(field.FindBrho(p0,*dir0[i],p2,*dir2[i])); + d0=(p0-pM).Unit(); +// p2.Print(); + if(d0.Z()>0.9) + h->Fill(field.FindBrho(p0,d0,p2,*dir2[i])); } - o0.Print(); - o2.Print(); - cout << h->GetStdDev() << endl; + cout << h->GetStdDev() << " " << parameter[0] << endl; return h->GetStdDev(); } //////////////////////////////////////////////////////////////////////////////// @@ -41,18 +46,18 @@ void LoadFile(){ cout << "Success opening file" << endl; } + double xm,ym,zm; double x0,y0,z0; - double dx0,dy0,dz0; double x2,y2,z2; double dx2,dy2,dz2; - while(file >> x0 >> y0 >> z0 >> dx0 >> dy0 >> dz0 >> x2 >> y2 >> z2 >> dx2 >> dy2 >> dz2){ + while(file >> x0 >> y0 >> z0 >> xm >> ym >> zm >> x2 >> y2 >> z2 >> dx2 >> dy2 >> dz2){ auto p0 = new TVector3(x0,y0,z0); - auto d0 = new TVector3(dx0,dy0,dz0); + auto pM = new TVector3(xm,ym,zm); auto p2 = new TVector3(x2,y2,z2); auto d2 = new TVector3(dx2,dy2,dz2); pos0.push_back(p0); - dir0.push_back(d0); + posM.push_back(pM); pos2.push_back(p2); dir2.push_back(d2); } @@ -69,28 +74,47 @@ void FDC(){ field.LoadMap(30*deg,"field_map/180702-2,40T-3000.table.bin",10); field.SetFDC2Angle((59.930-90.0)*deg); - double parameter[6]={0,0,-3456.52,-252.55,0,4123.47}; - cout << devB(parameter) << endl; - // Minimizer auto min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"); - auto func=ROOT::Math::Functor(&devB,6); + auto func=ROOT::Math::Functor(&devB,10); min->SetFunction(func); min->SetPrintLevel(0); - min->SetPrecision(1e-3); - min->SetLimitedVariable(0,"X0",0,1,-10,10); - min->SetLimitedVariable(1,"Y0",0,1,-10,10); - min->SetLimitedVariable(2,"Z0",-3372.07,1,-3500,-3400); - min->SetLimitedVariable(3,"X2",-252.55,1,-260,-240); - min->SetLimitedVariable(4,"Y2",0,1,-10,10); - min->SetLimitedVariable(5,"Z2",4123.47,1,4000,4200); + min->SetPrecision(1e-6); + double parameter[10]={40.6,0,0,-4657.39, 0, 0,-3372.07,-252.55,0,4123.47}; + devB(parameter); + + // min->SetLimitedVariable(0,"AM",parameter[0],1,-90,90); + min->SetFixedVariable(0,"AM",parameter[0]); + min->SetLimitedVariable(1,"XM",parameter[1],1,-10,10); + min->SetLimitedVariable(2,"YM",parameter[2],1,-10,10); + min->SetLimitedVariable(3,"ZM",parameter[3],1,-4650,-4660); + min->SetLimitedVariable(4,"X0",parameter[4],1,-10,10); + min->SetLimitedVariable(5,"Y0",parameter[5],1,-10,10); + min->SetLimitedVariable(6,"Z0",parameter[6],1,-3370,-3375); + min->SetLimitedVariable(7,"X2",parameter[7],1,-260,-240); + min->SetLimitedVariable(8,"Y2",parameter[8],1,-10,10); + min->SetLimitedVariable(9,"Z2",parameter[9],1,4120,4125); + /*min->SetFixedVariable(1,"XM",parameter[1]); + min->SetFixedVariable(2,"YM",parameter[2]); + min->SetFixedVariable(3,"ZM",parameter[3]); + min->SetFixedVariable(4,"X0",parameter[4]); + min->SetFixedVariable(5,"Y0",parameter[5]); + min->SetFixedVariable(6,"Z0",parameter[6]); + min->SetFixedVariable(7,"X2",parameter[7]); + min->SetFixedVariable(8,"Y2",parameter[8]); + min->SetFixedVariable(9,"Z2",parameter[9]); + */ min->Minimize(); const double* x = min->X(); - cout << "X0 =" << x[0]<<endl; - cout << "Y0 =" << x[1]<<endl; - cout << "Z0 =" << x[2]<<endl; - cout << "X2 =" << x[3]<<endl; - cout << "Y2 =" << x[4]<<endl; - cout << "Z2 =" << x[5]<<endl; - + cout << "AM =" << x[0]<<endl; + cout << "XM =" << x[1]<<endl; + cout << "YM =" << x[2]<<endl; + cout << "ZM =" << x[3]<<endl; + cout << "X0 =" << x[4]<<endl; + cout << "Y0 =" << x[5]<<endl; + cout << "Z0 =" << x[6]<<endl; + cout << "X2 =" << x[7]<<endl; + cout << "Y2 =" << x[8]<<endl; + cout << "Z2 =" << x[9]<<endl; + cout << "Minimum: " << devB(x) << endl; } diff --git a/Projects/S034/macro/rigz.cxx b/Projects/S034/macro/rigz.cxx index 19e4cba55..7da954ecf 100644 --- a/Projects/S034/macro/rigz.cxx +++ b/Projects/S034/macro/rigz.cxx @@ -4,9 +4,9 @@ void rigz(){ auto fz = new TFile("root/zaihong/run0582_RIG20210424_6He.root"); auto tz = (TTree*) fz->FindObjectAny("rig"); - auto fl = new TFile("root/analysis/Results582.root"); - auto tl = (TTree*) fl->FindObjectAny("ResultTree"); - tl->AddFriend(tz); + auto fl = new TFile("root/analysis/Result582.root"); + auto tl = (TTree*) fl->FindObjectAny("PhysicsTree"); + tl->AddFriend(tz); double FDC0_X,FDC0_Y,FDC2_X,FDC2_ThetaX,beta; int FragID; tz->SetBranchAddress("FDC0_X",&FDC0_X); @@ -24,11 +24,11 @@ void rigz(){ auto b = new TH1D("rig","rig",1000,2,8); unsigned int entries = tz->GetEntries(); - for(unsigned int i = 0 ; i < entries ; i++){ + for(unsigned int i = 0 ; i < entries; i++){ tz->GetEntry(i); double brho_param[6]={FDC0_X/*+1.77*/, FDC0_Y, 0, 0, FDC2_X/*-252.55*/, FDC2_ThetaX}; double Brho=r_fit(brho_param); - if(Brho>2&& Brho<8 && FragID>0 && FragID<27){ + if(FragID>0 && FragID<27){ h->Fill(Brho); // compute Brho based on beta and FragID double rig ; @@ -48,18 +48,25 @@ void rigz(){ He6.SetBeta(beta); rig = He6.GetBrho(); } - b->Fill(rig); + if(rig>2 && rig<8) + b->Fill(rig); } } - h->Scale(1./h->Integral()); - h->Draw(); h->SetLineColor(kBlack); - b->Draw("same"); +// h->Scale(1./h->Integral()); +// h->Draw(); h->SetLineColor(kBlack); + b->Draw(""); b->Scale(1./b->Integral()); - b->SetLineColor(kOrange+7);b->SetLineWidth(4); - cout << tl->Draw("Brho>>g","Brho>2 && Brho<8 &&FragID>0 && FragID<27","same") << endl; + b->SetLineColor(kOrange+7);b->SetLineWidth(2); + cout << "ratio Z: " << b->Integral(b->FindBin(3),b->FindBin(4.3))/ b->Integral(b->FindBin(4.5),b->FindBin(6.5)); + cout << tl->Draw("Brho>>g(1000,2,8)","Brho>2 && Brho<8 &&FragID>0 && FragID<27","same") << endl; + cout << tl->Draw("BrhoP>>hp(1000,2,8)","Brho>2 && Brho<8 &&FragID>0 && FragID<27","same") << endl; auto g = (TH1*) gDirectory->FindObjectAny("g"); - g->SetLineColor(kAzure+7);g->SetLineWidth(4); + auto hp = (TH1*) gDirectory->FindObjectAny("hp"); + g->SetLineColor(kAzure+7);g->SetLineWidth(2); g->Scale(1./g->Integral()); + cout << "ratio N: " << g->Integral(g->FindBin(3),g->FindBin(4.3))/ g->Integral(g->FindBin(4.5),g->FindBin(6.5)); + hp->Scale(1./hp->Integral()); + hp->SetLineColor(kBlack);b->SetLineWidth(2); auto l = new TLine(3.62,0,3.62,h->GetMaximum());l->Draw(); l = new TLine(5.53,0,5.53,h->GetMaximum());l->Draw(); l = new TLine(5.48,0,5.48,h->GetMaximum());l->Draw(); diff --git a/Projects/S034/s034.detector b/Projects/S034/s034.detector index f4763b912..f0c1a6b7d 100644 --- a/Projects/S034/s034.detector +++ b/Projects/S034/s034.detector @@ -1,8 +1,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SAMURAIBDC 1 XML= db/SAMURAIBDC1.xml -% Offset= -0.69 +1.1 -6876.34 mm - Offset= 0 0 -6876.34 mm + Offset= 0.959479 -0.17903 -6876.34 mm InvertX= 0 InvertY= 0 InvertD= 1 @@ -11,15 +10,16 @@ SAMURAIBDC 1 SAMURAIBDC 2 XML= db/SAMURAIBDC2.xml % Offset= +0.25 +1.6 -5876.7 mm - Offset= 0 0 -5876.7 mm +% Offset= 0 0 -5876.7 mm + Offset= -0.0686464 0.294288 -5876.7 mm InvertX= 0 InvertY= 0 InvertD= 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Minos - Position= 0 0 -4657.39 mm - ZRotation= 35 deg + Position= 0.345399 1.02061 -4650 mm + ZRotation= 40.6 deg TargetLength= 151.72 mm TargetMaterial= LH2 CellMaterial= Mylar @@ -33,10 +33,7 @@ Minos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SAMURAIFDC0 XML= db/SAMURAIFDC0_20200109.xml -%Offset= 3.005 1.864 0 mm -% Offset= 1.77 0 -3372.07 mm - Offset= 4.52872 -9.97899 -3403.03 mm - Offset= 0 0 0 mm + Offset= -0.00666226 0.102191 -3370.01 mm InvertX= 1 InvertY= 0 InvertD= 1 @@ -44,8 +41,7 @@ SAMURAIFDC0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SAMURAIFDC2 XML= db/SAMURAIFDC2.xml - Offset= -255.13 9.86883 4168.46 mm -% Offset= -252.55 0 4123.47 mm + Offset= -252.416 -0.228477 4122.57 mm InvertX= 0 InvertY= 1 InvertD= 1 -- GitLab