diff --git a/NPLib/Detectors/Strasse/TStrassePhysics.cxx b/NPLib/Detectors/Strasse/TStrassePhysics.cxx index 08b8d7dd7f4837a6cd1147e426a9ed8c7f96fa85..1fee4f854b5ddd334c412170ab0f727e1b9176d9 100644 --- a/NPLib/Detectors/Strasse/TStrassePhysics.cxx +++ b/NPLib/Detectors/Strasse/TStrassePhysics.cxx @@ -75,8 +75,10 @@ ClassImp(TStrassePhysics) Inner_PCB_BevelAngle= 60*deg; Inner_PCB_UpstreamWidth=1*cm; Inner_PCB_DownstreamWidth=2*mm; - Inner_PCB_MidWidth=2*mm; + Inner_PCB_MidWidth=1*mm; Inner_PCB_Thickness=3*mm; + Inner_PCB_Ledge = 1*mm ; + Inner_PCB_Step = 2*mm ; Inner_Wafer_TransverseStrips= 128; Inner_Wafer_LongitudinalStrips= 128; @@ -98,8 +100,10 @@ ClassImp(TStrassePhysics) Outer_PCB_BevelAngle= 60*deg; Outer_PCB_UpstreamWidth=1*cm; Outer_PCB_DownstreamWidth=2*mm; - Outer_PCB_MidWidth=2*mm; + Outer_PCB_MidWidth=1*mm; Outer_PCB_Thickness=3*mm; + Outer_PCB_Ledge = 1*mm ; + Outer_PCB_Step = 2*mm ; Outer_Wafer_TransverseStrips= 128; Outer_Wafer_LongitudinalStrips= 128; @@ -116,7 +120,8 @@ void TStrassePhysics::AddInnerDetector(double R, double Z, double Phi, double Sh //cout << ActiveWidth << " " << ActiveLength << " " << LongitudinalPitch << " " << TransversePitch << endl; // Vector C position of detector face center - TVector3 C(Shift,R,Z);// center of the whole detector, including PCB + double Recess = (Inner_PCB_Thickness-Inner_PCB_Step-Inner_Wafer_Thickness); + TVector3 C(Shift,R+Recess,Z);// center of the whole detector, including PCB C.RotateZ(-Phi); // Vector W normal to detector face (pointing to the back) @@ -212,7 +217,8 @@ void TStrassePhysics::AddOuterDetector(double R, double Z, double Phi, double Sh //cout << ActiveWidth << " " << ActiveLength << " " << LongitudinalPitch << " " << TransversePitch << endl; // Vector C position of detector face center - TVector3 C(Shift,R,Z);// center of the whole detector, including PCB + double Recess = (Inner_PCB_Thickness-Inner_PCB_Step-Inner_Wafer_Thickness); + TVector3 C(Shift,R+Recess,Z);// center of the whole detector, including PCB C.RotateZ(-Phi); // Vector W normal to detector face (pointing to the back) @@ -659,6 +665,8 @@ void TStrassePhysics::ReadConfiguration(NPL::InputParser parser) { "Inner_PCB_DownstreamWidth", "Inner_PCB_MidWidth", "Inner_PCB_Thickness", + "Inner_PCB_Ledge", + "Inner_PCB_Step", "Inner_Wafer_TransverseStrips", "Inner_Wafer_LongitudinalStrips", "Outer_Wafer_Length", @@ -675,8 +683,10 @@ void TStrassePhysics::ReadConfiguration(NPL::InputParser parser) { "Outer_PCB_DownstreamWidth", "Outer_PCB_MidWidth", "Outer_PCB_Thickness", + "Outer_PCB_Ledge", + "Outer_PCB_Step", "Outer_Wafer_TransverseStrips", - "Outer_Wafer_LongitudinalStrips" + "Outer_Wafer_LongitudinalStrips", }; if(blocks_info[0]->HasTokenList(info)){ @@ -697,6 +707,8 @@ void TStrassePhysics::ReadConfiguration(NPL::InputParser parser) { Inner_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Inner_PCB_DownstreamWidth","mm"); Inner_PCB_MidWidth = blocks_info[0]->GetDouble("Inner_PCB_MidWidth","mm"); Inner_PCB_Thickness = blocks_info[0]->GetDouble("Inner_PCB_Thickness","mm"); + Inner_PCB_Ledge = blocks_info[0]->GetDouble("Inner_PCB_Ledge","mm"); + Inner_PCB_Step = blocks_info[0]->GetDouble("Inner_PCB_Step","mm"); Outer_Wafer_Length = blocks_info[0]->GetDouble("Outer_Wafer_Length","mm"); Outer_Wafer_Width = blocks_info[0]->GetDouble("Outer_Wafer_Width","mm"); Outer_Wafer_Thickness = blocks_info[0]->GetDouble("Outer_Wafer_Thickness","mm"); @@ -713,7 +725,8 @@ void TStrassePhysics::ReadConfiguration(NPL::InputParser parser) { Outer_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Outer_PCB_DownstreamWidth","mm"); Outer_PCB_MidWidth = blocks_info[0]->GetDouble("Outer_PCB_MidWidth","mm"); Outer_PCB_Thickness = blocks_info[0]->GetDouble("Outer_PCB_Thickness","mm"); - + Outer_PCB_Ledge = blocks_info[0]->GetDouble("Outer_PCB_Ledge","mm"); + Outer_PCB_Step = blocks_info[0]->GetDouble("Outer_PCB_Step","mm"); } else{ diff --git a/NPLib/Detectors/Strasse/TStrassePhysics.h b/NPLib/Detectors/Strasse/TStrassePhysics.h index a50c24a59b4f0f1d49a566f9ed90c778ce19e033..49e035e051230fd4a9dad00c1cf3c7f5975cb15b 100644 --- a/NPLib/Detectors/Strasse/TStrassePhysics.h +++ b/NPLib/Detectors/Strasse/TStrassePhysics.h @@ -259,6 +259,8 @@ class TStrassePhysics : public TObject, public NPL::VDetector { double Inner_PCB_DownstreamWidth; double Inner_PCB_MidWidth; double Inner_PCB_Thickness; + double Inner_PCB_Ledge; + double Inner_PCB_Step; double Inner_Wafer_TransverseStrips; double Inner_Wafer_LongitudinalStrips; @@ -282,6 +284,8 @@ class TStrassePhysics : public TObject, public NPL::VDetector { double Outer_PCB_DownstreamWidth; double Outer_PCB_MidWidth; double Outer_PCB_Thickness; + double Outer_PCB_Ledge; + double Outer_PCB_Step; double Outer_Wafer_TransverseStrips; double Outer_Wafer_LongitudinalStrips; diff --git a/NPSimulation/Detectors/Strasse/Strasse.cc b/NPSimulation/Detectors/Strasse/Strasse.cc index 0ee0ab3daf49125fde7313d92332e070399c7085..b4d7a1c6df5d09dbb4547b957c8b4c5fdb11655f 100644 --- a/NPSimulation/Detectors/Strasse/Strasse.cc +++ b/NPSimulation/Detectors/Strasse/Strasse.cc @@ -369,7 +369,8 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){ -0.5*Inner_PCB_Thickness -offsetPCB2-0.05,0); - G4ThreeVector WaferShiftZ = G4ThreeVector(0,0,-Inner_PCB_UpstreamWidth-Inner_PCB_MidWidth); + //G4ThreeVector WaferShiftZ = G4ThreeVector(0,0,-Inner_PCB_UpstreamWidth-Inner_PCB_MidWidth); + G4ThreeVector WaferShiftZ = G4ThreeVector(0,0,CentralZOffset); new G4PVPlacement(new G4RotationMatrix(0,0,0), WaferShiftY+WaferShiftZ // Shift along Y @@ -504,6 +505,7 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){ PCB2CrossSection.push_back(G4TwoVector(-Outer_PCB2_Width*0.5,Outer_PCB2_Thickness*0.5)); PCB2CrossSection.push_back(G4TwoVector(-Outer_PCB2_Width*0.5,-Outer_PCB2_Thickness*0.5)); + double Outer_PCB2_Length= 2*Outer_Wafer_Length + Outer_PCB_MidWidth; G4ExtrudedSolid* PCB2Full = @@ -513,14 +515,9 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){ G4TwoVector(0,0),1,// offset, scale G4TwoVector(0,0),1);// offset, scale - // Offset along beam axis between PCB middle and (2*Wafer+MiddleBar) volume center - double CentralZOffset = - Outer_PCB_Length*0.5 - + Outer_PCB_UpstreamWidth - + Outer_Wafer_Length - + Outer_PCB_MidWidth*0.5; double Length_Shift21 = -0.5*Outer_PCB_Length // Flush to border - + 0.5*(Outer_PCB_UpstreamWidth+Outer_PCB_DownstreamWidth) // add Upstream side shift + + 0.5*(Outer_PCB_UpstreamWidth+Outer_PCB_DownstreamWidth) + 0.5*Outer_Wafer_Length; double Length_Shift22 = Length_Shift21 // overlap detector 1 + Outer_Wafer_Length // at opposing edge @@ -546,7 +543,6 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){ Outer_PCB2_Thickness, Outer_PCB2_MidWidth*0.5); - // Substracting the hole Shape from the Stock PCB G4SubtractionSolid* PCB2_3 = new G4SubtractionSolid("PCB2_3", PCB2_2, HoleShapeCenterBar, new G4RotationMatrix,HoleCenterBar); @@ -555,13 +551,19 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){ new G4LogicalVolume(PCB2_3,m_MaterialPCB,"logicPCB2", 0, 0, 0); logicPCB2->SetVisAttributes(PADVisAtt); + // Offset along beam axis between PCB middle and (2*Wafer+MiddleBar) volume center + double CentralZOffset = - Outer_PCB_Length*0.5 + + Outer_PCB_UpstreamWidth + + Outer_Wafer_Length + + Outer_PCB_MidWidth*0.5; + new G4PVPlacement(new G4RotationMatrix(0,0,0), - //G4ThreeVector(0,-0.5*offsetPCB2,Outer_PCB_UpstreamWidth-Outer_PCB_MidWidth), G4ThreeVector(0,-0.5*offsetPCB2,CentralZOffset), logicPCB2,"Strasse_Outer_PCB2",m_OuterDetector, false,0); ////////////////////////////////////////////////////////////////// + // Build the Wafer // Sub volume Wafer G4Box* WaferShape = new G4Box("WaferShape", Outer_Wafer_Width*0.5, @@ -576,19 +578,24 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){ new G4LogicalVolume(WaferShape,m_MaterialSilicon,"logicWafer2", 0, 0, 0); logicWafer2->SetVisAttributes(GuardRingVisAtt); + // Shift along Y to flush the wafer to the pcb ledge on one side + G4ThreeVector WaferShiftY = G4ThreeVector(0,-0.5*Outer_Wafer_Thickness + -Outer_Wafer_AlThickness + -0.5*Outer_PCB_Thickness + -offsetPCB2-0.05,0); + + //G4ThreeVector WaferShiftZ = G4ThreeVector(0,0,-Outer_PCB_UpstreamWidth-Outer_PCB_MidWidth); + G4ThreeVector WaferShiftZ = G4ThreeVector(0,0,CentralZOffset); + new G4PVPlacement(new G4RotationMatrix(0,0,0), - G4ThreeVector(0,0.5*Outer_Wafer_Thickness - +Outer_Wafer_AlThickness - -0.5*Outer_PCB_Thickness,0)// flush the wafer to the pcb on one side - +HoleShift1, // Shift wafer in the hole + WaferShiftY+WaferShiftZ + +HoleShift21, // Shift along Z to put wafer in the 1st hole logicWafer1,"Strasse_Outer_Wafer1",m_OuterDetector, false,0); new G4PVPlacement(new G4RotationMatrix(0,0,0), - G4ThreeVector(0,0.5*Outer_Wafer_Thickness - +Outer_Wafer_AlThickness - -0.5*Outer_PCB_Thickness,0)// flush the wafer to the pcb on one side - +HoleShift2, // Shift wafer in the hole + WaferShiftY+WaferShiftZ + +HoleShift22, // Shift along Z to put wafer in the 1st hole logicWafer2,"Strasse_Outer_Wafer2",m_OuterDetector, false,0); @@ -608,14 +615,13 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){ logicActiveWafer2->SetVisAttributes(SiliconVisAtt); logicActiveWafer2->SetSensitiveDetector(m_OuterScorer2); - new G4PVPlacement(new G4RotationMatrix(0,0,0), - G4ThreeVector(0,0,0.5*(Outer_Wafer_PADExternal-Outer_Wafer_PADInternal)), // assymetric pading for bounding + G4ThreeVector(0,0,0.5*(Outer_Wafer_PADExternal-Outer_Wafer_PADInternal)), logicActiveWafer1,"Strasse_Outer_ActiveWafer1",logicWafer1, false,1); new G4PVPlacement(new G4RotationMatrix(0,0,0), - G4ThreeVector(0,0,-0.5*(Outer_Wafer_PADExternal-Outer_Wafer_PADInternal)), // assymetric pading for bounding + G4ThreeVector(0,0,-0.5*(Outer_Wafer_PADExternal-Outer_Wafer_PADInternal)), logicActiveWafer2,"Strasse_Outer_ActiveWafer2",logicWafer2, false,2); diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx index fbbd5e06dce0bde9e0e5772dd8646efc8cf5a350..61182b5f30ee430848c6c40732b61c3d7db759f8 100644 --- a/Projects/Strasse/Analysis.cxx +++ b/Projects/Strasse/Analysis.cxx @@ -28,6 +28,8 @@ using namespace std; #include"NPFunction.h" #include"NPTrackingUtility.h" #include"NPPhysicalConstants.h" +#include"TRandom3.h" + //////////////////////////////////////////////////////////////////////////////// Analysis::Analysis(){ } @@ -45,7 +47,7 @@ void Analysis::Init(){ InitInputBranch(); Strasse = (TStrassePhysics*) m_DetectorManager -> GetDetector("Strasse"); - Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana"); +// Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana"); // reaction properties m_QFS = new NPL::QFS(); m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile()); @@ -63,14 +65,27 @@ void Analysis::Init(){ protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100); protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100); LV_T.SetVectM(TVector3(0,0,0),NPUNITS::proton_mass_c2); + + rand= new TRandom3(); } //////////////////////////////////////////////////////////////////////////////// void Analysis::TreatEvent(){ // Reinitiate calculated variable ReInitValue(); + + bool perfectBeamEnergy=false; + bool perfectBeamTracking=false; + bool perfectProtonTracking=false; + bool perfectProtonEnergy=false; + bool CATANAEnergyReco=false; + unsigned int size = Strasse->GetEventMultiplicity(); if(size==2){ // 2 proton detected + + //////////////////////////////////////////////// + ////////////Protons (E,P) Reconstruction /////// + //////////////////////////////////////////////// // Proton 1 TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0); TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0); @@ -86,11 +101,6 @@ void Analysis::TreatEvent(){ sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg; Theta12 = Proton1.Angle(Proton2)/deg; - // reject event that make no physical sense - /*if(deltaPhi<170 && sumTheta<80){ - return; - } - */ // computing minimum distance of the two lines TVector3 Vertex; TVector3 delta; @@ -101,12 +111,79 @@ void Analysis::TreatEvent(){ deltaX=delta.X(); deltaY=delta.Y(); deltaZ=delta.Z(); - + + if(perfectProtonTracking){ + // From Reaction Conditions + VertexX=RC->GetVertexPositionX(); + VertexY=RC->GetVertexPositionY(); + VertexZ=RC->GetVertexPositionZ(); + Proton2=RC->GetParticleMomentum(0).Unit(); + Proton1=RC->GetParticleMomentum(1).Unit(); + deltaX=0; + deltaY=0; + deltaZ=0; + } + if(perfectProtonEnergy){ + E1 = RC->GetKineticEnergy(1); + E2 = RC->GetKineticEnergy(0); + } + else { + //General CATANA resolution for proton (1.4 % FWHM) + // simply applied to E, crude approximation (no CATANA reco.) + E1 = ApplyCATANAResoProton(RC->GetKineticEnergy(1)); + E2 = ApplyCATANAResoProton(RC->GetKineticEnergy(0)); + } + if(CATANAEnergyReco){ + //if true, bypass previous proton energy calc. and use full CATANA Reco. + /////////// + // TO DO // + /////////// + } + double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2)); + double P2= sqrt(E2*(E2+2*NPUNITS::proton_mass_c2)); + + /////////////////////////////////////// + //////Beam (E,P) Reconstruction /////// + /////////////////////////////////////// + + TVector3 BeamDirection; + if(perfectBeamTracking) BeamDirection = RC->GetBeamDirection(); + else BeamDirection = TVector3(0,0,1); + //Beam Energy at vertex + TA = RC->GetBeamEnergy(); + if(perfectBeamEnergy) TAcalc = TA; + else TAcalc = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ+75), BeamDirection.Angle(TVector3(0,0,1))); + double beam_mom=sqrt(TAcalc*(TAcalc+2*m_QFS->GetParticleA()->Mass())); + TVector3 PA=beam_mom*BeamDirection.Unit(); + + ////////////////////////////////////////////// + ///// Missing mass using Lorentz Vector ////// + ////////////////////////////////////////////// + + LV_A.SetVectM(PA,m_QFS->GetParticleA()->Mass()); + LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2); + LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2); + // computing Ex from Missing Mass + LV_B = LV_A + LV_T - LV_p1 - LV_p2; + //LV_B = RC->GetParticleMomentum(2); + Ex = LV_B.M() - m_QFS->GetParticleB()->Mass(); + + }//endif size==2 +} + +/* /////////////////////////////////// + * //OLD Tentative CATANA treatment/// + * /////////////////////////////////// + * // Look for associated Catana event double d1,d2; unsigned int i1,i2; i1 = Catana->FindClosestHitToLine(InnerPos1,OuterPos1,d1); i2 = Catana->FindClosestHitToLine(InnerPos2,OuterPos2,d2); +cout<<"------------------------"<<endl; +cout<<"d1="<<d1<<endl; +cout<<"d2="<<d2<<endl; +cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl; //if(i1!=i2){ if(true){ double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),RC->GetBeamDirection().Angle(TVector3(0,0,1))); @@ -118,6 +195,7 @@ void Analysis::TreatEvent(){ TVector3 Proton1s=RC->GetParticleMomentum(0).Unit(); TVector3 Proton2s=RC->GetParticleMomentum(1).Unit(); // Matching the right energy with the right proton + if((Proton1s-Proton1).Mag()<(Proton1s-Proton2).Mag()){ E1=E1s; E2=E2s; @@ -145,43 +223,28 @@ void Analysis::TreatEvent(){ Theta2s=Proton1s.Theta()/deg; Phi2s=Proton1s.Phi()/deg; } + // From detectors E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]); E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]); + //E1 = Catana->Energy[i1]; + //E2 = Catana->Energy[i2]; + cout<<"------------------------"<<endl; - // Vertex = RC->GetVertexPosition(); - //TA = RC->GetBeamEnergy(); - //////////////////////////////////// - - // setting up Lorentz Vector from measured trajectories and energies - // TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetParticleA()->Mass()))); // for like there is no BDC - double beam_mom=sqrt(TA*(TA+2*m_QFS->GetParticleA()->Mass())); - // TVector3 PA(0,0,beam_mom); // for like there is no BDC - TVector3 PA=beam_mom*RC->GetBeamDirection().Unit(); - - LV_A.SetVectM(PA,m_QFS->GetParticleA()->Mass()); - double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2)); - double P2= sqrt(E2*(E2+2*NPUNITS::proton_mass_c2)); - - LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2); - LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2); - // computing Ex from Missing Mass - LV_B = LV_A + LV_T - LV_p1 - LV_p2; - //LV_B = RC->GetParticleMomentum(2); - Ex = LV_B.M() - m_QFS->GetParticleB()->Mass(); - } - } -} +*/ //////////////////////////////////////////////////////////////////////////////// double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir,const double& Ecatana){ + TVector3 Normal = TVector3(0,1,0); Normal.SetPhi(dir.Phi()); double Theta = dir.Angle(Normal); // Catana Al housing double E = protonAl.EvaluateInitialEnergy(Ecatana,0.5*mm,Theta); + cout<<"Ecatana="<<Ecatana<<endl; + cout<<"Erec0="<<E<<endl; // Strasse Chamber - E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta); + //E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta); // Outer Barrel E = protonSi.EvaluateInitialEnergy(E,300*micrometer,Theta); // Inner Barrel @@ -193,11 +256,40 @@ double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir TVector3 T2(0,15,1); T1.SetPhi(dir.Phi()); T2.SetPhi(dir.Phi()); - double d = NPL::MinimumDistancePointLine(T1,T2,x0); + // double d = NPL::MinimumDistancePointLine(T1,T2,x0); + double d = 0; + + cout<<"d="<<d<<endl; + cout<<"Theta="<<Theta<<endl; E = protonTarget.EvaluateInitialEnergy(E,d,Theta); - } + cout<<"Erec="<<E<<endl; + return E; + + //return 0; +} +//////////////////////////////////////////////////////////////////////////////// +double Analysis::ApplyCATANAResoProton(double Ein){ + // Function applying overall CATANA crystal resolution + // For protons (1.4% FWHM) + double FWHM = 1.4/100 * Ein; + double sigma = FWHM/2.35; + double Eout = rand->Gaus(Ein,sigma); + return Eout; +} +//////////////////////////////////////////////////////////////////////////////// +double Analysis::ApplyCATANAResoGamma(double Ein){ + // Function applying overall CATANA crystal resolution + // For gammas defined in smsimulator package + double a = 0.686569; + double b = 0.564352; + // Ein from MeV to keV + Ein = Ein * 1000; + double SigmaE = a * pow(Ein,b); + double Eout = rand->Gaus(Ein,SigmaE); + return Eout/1000.; +} //////////////////////////////////////////////////////////////////////////////// TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj){ TVector3 Vproj(-999,-999,-999); @@ -207,7 +299,6 @@ TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj){ Vproj.SetXYZ(Xproj,Yproj,Zproj); return Vproj; } - //////////////////////////////////////////////////////////////////////////////// void Analysis::End(){ } @@ -237,6 +328,8 @@ void Analysis::InitOutputBranch() { RootOutput::getInstance()->GetTree()->Branch("Phi1s",&Phi1s,"Phi1s/D"); RootOutput::getInstance()->GetTree()->Branch("Theta2s",&Theta2s,"Theta2s/D"); RootOutput::getInstance()->GetTree()->Branch("Phi2s",&Phi2s,"Phi2s/D"); + RootOutput::getInstance()->GetTree()->Branch("TA",&TA,"TA/D"); + RootOutput::getInstance()->GetTree()->Branch("TAcalc",&TAcalc,"TAcalc/D"); RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D"); @@ -274,6 +367,8 @@ void Analysis::ReInitValue(){ Phi1s=-1000; Theta2s=-1000; Phi2s=-1000; + TA=-1000; + TAcalc=-1000; } diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h index 9754678afa6b9f134ce44a6707941336bd9291b8..841ff583ff7b18c190b1533d20e8ce7dc0961c0f 100644 --- a/Projects/Strasse/Analysis.h +++ b/Projects/Strasse/Analysis.h @@ -50,11 +50,15 @@ class Analysis: public NPL::VAnalysis{ void ReInitValue(); static NPL::VAnalysis* Construct(); TVector3 InterpolateInPlaneZ(TVector3,TVector3,double); + double ApplyCATANAResoGamma(double); + double ApplyCATANAResoProton(double); private: double Ex; double E1; double E2; + double E1s; + double E2s; double Theta12; double ThetaCM; double VertexX; @@ -75,8 +79,12 @@ class Analysis: public NPL::VAnalysis{ double Phi1s=-1000; double Theta2s=-1000; double Phi2s=-1000; + double TA; + double TAcalc; + TRandom3 *rand; + TLorentzVector LV_A; TLorentzVector LV_T; TLorentzVector LV_B; diff --git a/Projects/Strasse/ana.sh b/Projects/Strasse/ana.sh new file mode 100755 index 0000000000000000000000000000000000000000..924f3f65240aa504f7ca8f291fb156cf2d3c1e60 --- /dev/null +++ b/Projects/Strasse/ana.sh @@ -0,0 +1 @@ +npanalysis -E ./reaction/17Fp2p.reaction -D ./geometry/strasse_July2021.detector -R RunToTreat.txt -O configJuly2021_ana.root diff --git a/Projects/Strasse/geometry/strasse_CAD.detector b/Projects/Strasse/geometry/strasse_CAD.detector index 1acb1346fd57bcc8452680ae3d0d678122e51f2d..28aeff76d5271140d06d0f94644b3379f047dbc0 100644 --- a/Projects/Strasse/geometry/strasse_CAD.detector +++ b/Projects/Strasse/geometry/strasse_CAD.detector @@ -7,7 +7,7 @@ Target X= 0 mm Y= 0 mm Z= 0 mm - NbSlices= 10 + NbSlices= 150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 Strasse Info @@ -28,7 +28,7 @@ Strasse Info Inner_PCB_MidWidth= 1 mm Inner_PCB_Thickness= 2.4 mm Inner_PCB_Ledge= 1 mm - Inner_PCB_Step= 2 mm + Inner_PCB_Step= 2.2 mm Outer_Wafer_Length= 123 mm Outer_Wafer_Width= 64.6 mm Outer_Wafer_Thickness= 300 micrometer @@ -44,7 +44,7 @@ Strasse Info Outer_PCB_MidWidth= 1 mm Outer_PCB_Thickness= 2.4 mm Outer_PCB_Ledge= 1 mm - Outer_PCB_Step= 2 mm + Outer_PCB_Step= 2.1 mm Outer_Wafer_TransverseStrips= 605 Outer_Wafer_LongitudinalStrips= 313 % unused if using CAD file (.stl) chamber @@ -87,11 +87,11 @@ Strasse Outer Ref= 0 0 0 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Strasse InactiveMaterial - Chamber= ./geometry/STRASSE_Chamber.stl - Stars= ./geometry/STRASSE_StarSupports.stl - Base= ./geometry/STRASSE_Base.stl - Blades= ./geometry/STRASSE_Blades.stl +%Strasse InactiveMaterial +% Chamber= ./geometry/STRASSE_Chamber.stl +% Stars= ./geometry/STRASSE_StarSupports.stl +% Base= ./geometry/STRASSE_Base.stl +% Blades= ./geometry/STRASSE_Blades.stl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 %Catana CSV % Path= geometry/Catana.csv diff --git a/Projects/Strasse/geometry/strasse_July2021.detector b/Projects/Strasse/geometry/strasse_July2021.detector new file mode 100644 index 0000000000000000000000000000000000000000..aa02386b2bb1003f5c2b8cf8b27c428ec3679a74 --- /dev/null +++ b/Projects/Strasse/geometry/strasse_July2021.detector @@ -0,0 +1,101 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Target + THICKNESS= 150 mm + ANGLE= 0 deg + RADIUS= 15 mm + MATERIAL= LH2 + X= 0 mm + Y= 0 mm + Z= 0 mm + NbSlices= 150 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +Strasse Info + Inner_Wafer_Length= 124 mm + Inner_Wafer_Width= 32 mm + Inner_Wafer_Thickness= 200 micrometer + Inner_Wafer_AlThickness= 0.4 micrometer + Inner_Wafer_PADExternal= 0 mm + Inner_Wafer_PADInternal= 0 mm + Inner_Wafer_GuardRing= 1.0 mm + Inner_Wafer_TransverseStrips= 610 + Inner_Wafer_LongitudinalStrips= 150 + Inner_PCB_PortWidth= 1 mm + Inner_PCB_StarboardWidth= 1 mm + Inner_PCB_BevelAngle= 90 deg + Inner_PCB_UpstreamWidth= 12 mm + Inner_PCB_DownstreamWidth= 38 mm + Inner_PCB_MidWidth= 1 mm + Inner_PCB_Thickness= 2.4 mm + Inner_PCB_Ledge= 1 mm + Inner_PCB_Step= 2.2 mm + Outer_Wafer_Length= 123 mm + Outer_Wafer_Width= 64.6 mm + Outer_Wafer_Thickness= 300 micrometer + Outer_Wafer_AlThickness= 0.4 micrometer + Outer_Wafer_PADExternal= 0 mm + Outer_Wafer_PADInternal= 0 mm + Outer_Wafer_GuardRing= 1.0 mm + Outer_PCB_PortWidth= 1 mm + Outer_PCB_StarboardWidth= 1 mm + Outer_PCB_BevelAngle= 90 deg + Outer_PCB_UpstreamWidth= 40 mm + Outer_PCB_DownstreamWidth= 12 mm + Outer_PCB_MidWidth= 1 mm + Outer_PCB_Thickness= 2.4 mm + Outer_PCB_Ledge= 1 mm + Outer_PCB_Step= 2.1 mm + Outer_Wafer_TransverseStrips= 605 + Outer_Wafer_LongitudinalStrips= 313 + % unused if using CAD file (.stl) chamber + Chamber_Thickness= 3 mm + Chamber_Cylinder_Length= 360 mm + Chamber_Radius= 180 mm + Chamber_ExitTube_Radius= 79.5 mm + Chamber_ExitTube_Length= 100 mm + Chamber_Flange_Inner_Radius= 50 mm + Chamber_Sphere_Radius= 220 mm + Chamber_Sphere_Shift= 60 mm + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Alias InnerPhi + Action= Copy + %Value= 0 + Value= -5.5 54.5 114.5 174.5 234.5 294.5 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Strasse Inner + Radius= 28.7 mm + Z= 76.5 mm + Phi= @InnerPhi deg + Shift= 3.5 mm + Ref= 0 0 0 mm + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Alias OuterPhi + Action= Copy + %Value= 0 + Value= 0 60 120 180 240 300 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Strasse Outer + Radius= 60.7 mm + Z= 76.5 mm + Phi= @OuterPhi deg + Shift= 0 mm + Ref= 0 0 0 mm + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Strasse InactiveMaterial +% Chamber= ./geometry/STRASSE_Chamber.stl +% Stars= ./geometry/STRASSE_StarSupports.stl +% Base= ./geometry/STRASSE_Base.stl +% Blades= ./geometry/STRASSE_Blades.stl +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +%Catana CSV +% Path= geometry/Catana.csv +% Pos= 0 0 100 mm +% Rshift= 100 micrometer + + diff --git a/Projects/Strasse/geometry/strasse_ref.detector b/Projects/Strasse/geometry/strasse_ref.detector index c542cea82b01d10d2d9cb363efe2cbae423ed81a..34317b25dda2d3743ab83addeb05cd8b970bc051 100644 --- a/Projects/Strasse/geometry/strasse_ref.detector +++ b/Projects/Strasse/geometry/strasse_ref.detector @@ -7,12 +7,12 @@ Target X= 0 mm Y= 0 mm Z= 0 mm - NbSlices= 10 + NbSlices= 150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 Strasse Info Inner_Wafer_Length= 127 mm - Inner_Wafer_Width= 33 mm + Inner_Wafer_Width= 30.6 mm Inner_Wafer_Thickness= 200 micrometer Inner_Wafer_AlThickness= 0.4 micrometer Inner_Wafer_PADExternal= 0 mm @@ -22,13 +22,13 @@ Strasse Info Inner_Wafer_LongitudinalStrips= 150 Inner_PCB_PortWidth= 1 mm Inner_PCB_StarboardWidth= 1 mm - Inner_PCB_BevelAngle= 45 deg + Inner_PCB_BevelAngle= 90 deg Inner_PCB_UpstreamWidth= 1 mm Inner_PCB_DownstreamWidth= 1 mm Inner_PCB_MidWidth= 1 mm Inner_PCB_Thickness= 1.6 mm Outer_Wafer_Length= 124 mm - Outer_Wafer_Width= 68 mm + Outer_Wafer_Width= 65.6 mm Outer_Wafer_Thickness= 300 micrometer Outer_Wafer_AlThickness= 0.4 micrometer Outer_Wafer_PADExternal= 0 mm @@ -43,6 +43,7 @@ Strasse Info Outer_PCB_Thickness= 1.6 mm Outer_Wafer_TransverseStrips= 605 Outer_Wafer_LongitudinalStrips= 325 + % all radius are external, internal deduced from thickness Chamber_Thickness= 3 mm Chamber_Cylinder_Length= 360 mm Chamber_Radius= 180 mm @@ -51,7 +52,6 @@ Strasse Info Chamber_Flange_Inner_Radius= 50 mm Chamber_Sphere_Radius= 220 mm Chamber_Sphere_Shift= 60 mm - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Alias InnerPhi Action= Copy @@ -63,6 +63,7 @@ Strasse Inner Z= 66.0 mm Phi= @InnerPhi deg Shift= 0 mm + Ref= 0 0 0 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Alias OuterPhi @@ -75,15 +76,15 @@ Strasse Outer Z= 78.0 mm Phi= @OuterPhi deg Shift= 0 mm - + Ref= 0 0 0 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Strasse Chamber - Z= -30 mm +%Strasse Chamber +% Z= -30 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Catana CSV - Path= geometry/Catana.csv - Pos= 0 0 100 mm - Rshift= 100 micrometer +%Catana CSV +% Path= geometry/Catana.csv +% Pos= 0 0 100 mm +% Rshift= 100 micrometer diff --git a/Projects/Strasse/macro/CalculateDetectorOffset.cxx b/Projects/Strasse/macro/CalculateDetectorOffset.cxx index 97dab461602c0cde1e36d5f7fcd9b6080679b92a..bef189e72f89ebf42768c79884c3ab4e9ca01e19 100644 --- a/Projects/Strasse/macro/CalculateDetectorOffset.cxx +++ b/Projects/Strasse/macro/CalculateDetectorOffset.cxx @@ -12,7 +12,7 @@ using namespace std; using namespace NPL; -void CalculateDetectorOffset(const char * fname = "strasse_optimized"){ +void CalculateDetectorOffset(const char * fname = "./geometry/strasse_July2021"){ // Open output ROOT file from NPTool simulation run string path = ""; @@ -214,9 +214,9 @@ double InnerTotalLength = double d_TargetCenter_InnerCenter = - TargetThickness/2. + d_TargetFront_InnerActive - + InnerTotalLength/2. - Inner_Wafer_GuardRing - - Inner_PCB_UpstreamWidth; + - Inner_PCB_UpstreamWidth + + InnerTotalLength/2.; double OuterTotalLength = Outer_PCB_UpstreamWidth diff --git a/Projects/Strasse/macro/CheckAna.cxx b/Projects/Strasse/macro/CheckAna.cxx new file mode 100644 index 0000000000000000000000000000000000000000..6e32e08d62582008b20605942d9b6c755090f479 --- /dev/null +++ b/Projects/Strasse/macro/CheckAna.cxx @@ -0,0 +1,30 @@ +{ +gStyle->SetPalette(1); +//TFile *file= new TFile("../../Outputs/Analysis/PhysicsTree.root"); +TFile *file= new TFile("../../Outputs/Analysis/configJuly2021_ana.root"); +TTree *tree = (TTree*)file->Get("PhysicsTree"); + +TCanvas *c1 = new TCanvas("c1","c1",1000,1000); +c1->Divide(3,3); +c1->cd(1); +tree->Draw("Ex>>h1(200,-10,10)","",""); +c1->cd(2); +tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(400,-150,250,200,-100,100)","","colz"); +c1->cd(3); +tree->Draw("VertexX:VertexZ>>h4(400,-150,250,200,-100,100)","","colz"); +c1->cd(4); +tree->Draw("fRC_Vertex_Position_X-VertexX>>h5(100,-2,2)","",""); +c1->cd(5); +tree->Draw("fRC_Vertex_Position_Y-VertexY>>h6(100,-2,2)","",""); +c1->cd(6); +tree->Draw("fRC_Vertex_Position_Z-VertexZ>>h7(100,-2,2)","",""); +c1->cd(7); +//tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[0]>>h8(100,-10,10)","Catana.GetEventMultiplicity()==2",""); +tree->Draw("fRC_Kinetic_Energy[1]-E1>>h8(100,-50,50)","",""); +c1->cd(8); +//tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[1]>>h9(100,-10,10)","Catana.GetEventMultiplicity()==2",""); +tree->Draw("fRC_Kinetic_Energy[0]-E2>>h9(100,-50,50)","",""); +c1->cd(9); +tree->Draw("fRC_Kinetic_Energy[0]:E2>>h10(300,0,300,300,0,300)","","colz"); + +} diff --git a/Projects/Strasse/macro/CheckSim.cxx b/Projects/Strasse/macro/CheckSim.cxx index 5b0f07c253c6f4e082524baacebab80aa888f3a7..c7e2a762f830f0799da52eae67f30aa19867125d 100644 --- a/Projects/Strasse/macro/CheckSim.cxx +++ b/Projects/Strasse/macro/CheckSim.cxx @@ -1,27 +1,29 @@ { gStyle->SetPalette(1); -//TFile *file= new TFile("../../Outputs/Simulation/test_ref.root"); -TFile *file= new TFile("../../Outputs/Simulation/test_optimized.root"); +//TFile *file= new TFile("../../Outputs/Simulation/SimulatedTree.root"); +TFile *file= new TFile("../../Outputs/Simulation/configJuly2021.root"); TTree *tree = (TTree*)file->Get("SimulatedTree"); TCanvas *c1 = new TCanvas("c1","c1",1000,1000); c1->Divide(4,4); c1->cd(1); -tree->Draw("fRC_Beam_Reaction_Energy:fRC_Vertex_Position_Z>>h1(200,-3.1,3.1,220,,)","","colz"); +tree->Draw("fRC_Beam_Reaction_Energy:fRC_Vertex_Position_Z>>h1(200,-100,100,650,3200,4500)","","colz"); c1->cd(2); -tree->Draw("fRC_Vertex_Position_Y:fRC_Vertex_Position_Z>>h2","",""); +tree->Draw("fRC_Vertex_Position_Y:fRC_Vertex_Position_Z>>h2(200,-100,100,120,-60,60)","","colz"); c1->cd(3); -tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3","",""); +tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(200,-100,100,120,-60,60)","","colz"); c1->cd(4); tree->Draw("fDetected_Position_Y:fDetected_Position_X>>h4","","colz"); c1->cd(5); -tree->Draw("fDetected_Position_Y:fDetected_Position_Z>>h5(1200,-70,230,280,-70,70)","","colz"); +tree->Draw("fDetected_Position_Y:fDetected_Position_Z>>h5(900,-70,230,700,-70,70)","","colz"); c1->cd(6); tree->Draw("fInner_TE_StripNbr:fDetected_Position_Z>>h6","","colz"); c1->cd(7); tree->Draw("fInner_LE_StripNbr:fDetected_Position_X>>h7","","colz"); c1->cd(8); tree->Draw("Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()>>h8(6,0,6)","",""); +cout<<"MULT4="<<h8->GetBinContent(5)<<endl; +cout<<"Efficiency2p="<<h8->GetBinContent(5)/50000*100<<endl; c1->cd(9); tree->Draw("fInner_LE_StripNbr>>h9(1250,0,1250)","",""); c1->cd(10); @@ -45,5 +47,34 @@ tree->Draw("fRC_Theta[0]>>h18","Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMu h18->SetLineColor(2); h18->Scale(1./4);; +/* +TCanvas *c2 = new TCanvas("c2","c2",1000,1000); +c2->Divide(2,2); +c2->cd(1); +tree->Draw("Strasse.GetInnerMultTEnergy()>>hMultInnerT(4,0,4)","",""); +c2->cd(2); +tree->Draw("Strasse.GetInnerMultLEnergy()>>hMultInnerL(4,0,4)","",""); +c2->cd(3); +tree->Draw("Strasse.GetOuterMultTEnergy()>>hMultOuterT(4,0,4)","",""); +c2->cd(4); +tree->Draw("Strasse.GetOuterMultLEnergy()>>hMultOuterL(4,0,4)","",""); + +TCanvas *c3 = new TCanvas("c3","c3",1000,1000); +c3->Divide(3,2); + +c3->cd(1); +tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta(900,0,90,1000,0,350)","","colz"); +c3->cd(2); +tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerT(900,0,90,1000,0,350)","Strasse.GetInnerMultTEnergy()==1","colz"); +c3->cd(3); +tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerL(900,0,90,1000,0,350)","Strasse.GetInnerMultLEnergy()==1","colz"); +c3->cd(4); +tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterT(900,0,90,1000,0,350)","Strasse.GetOuterMultTEnergy()==1","colz"); +c3->cd(5); +tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterL(900,0,90,1000,0,350)","Strasse.GetOuterMultLEnergy()==1","colz"); +c3->cd(6); +tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultAll(900,0,90,1000,0,350)","Strasse.GetInnerMultTEnergy()==1 && Strasse.GetInnerMultLEnergy()==1 && Strasse.GetOuterMultLEnergy()==1 && Strasse.GetOuterMultTEnergy()==1","colz"); + +*/ } diff --git a/Projects/Strasse/sim.sh b/Projects/Strasse/sim.sh new file mode 100755 index 0000000000000000000000000000000000000000..fcbe616643eb6f91d80c6e9888f5db67c2cdff02 --- /dev/null +++ b/Projects/Strasse/sim.sh @@ -0,0 +1,2 @@ +npsimulation -E ./reaction/17Fp2p.reaction -D ./geometry/strasse_July2021.detector -O configJuly2021.root -B run.mac +