From effdbe761773352033bd2a7f71adb175325c68c2 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS-MAC" Date: Mon, 10 Apr 2023 15:46:24 -0400 Subject: [PATCH] modified simulation to adopt 2nd array --- Cleopatra/Check_Simulation.C | 28 +- Cleopatra/HELIOS_LIB.h | 425 ++++-------------- Cleopatra/Simulation_Helper.C | 107 +++-- armory/AnalysisLib.h | 183 +++++--- working/{DWBA => DWBA2} | 0 working/{Ex.txt => Ex1.txt} | 0 working/Ex2.txt | 6 + working/Monitor.C | 3 +- working/detectorGeo.txt | 2 +- ...reactionConfig.txt => reactionConfig1.txt} | 14 +- working/reactionConfig2.txt | 25 ++ 11 files changed, 348 insertions(+), 445 deletions(-) rename working/{DWBA => DWBA2} (100%) rename working/{Ex.txt => Ex1.txt} (100%) create mode 100644 working/Ex2.txt rename working/{reactionConfig.txt => reactionConfig1.txt} (72%) create mode 100644 working/reactionConfig2.txt diff --git a/Cleopatra/Check_Simulation.C b/Cleopatra/Check_Simulation.C index c6bdaee..f36aca4 100644 --- a/Cleopatra/Check_Simulation.C +++ b/Cleopatra/Check_Simulation.C @@ -47,7 +47,7 @@ enum plotID { pEZ, /// 0 pElum1RThetaCM, /// 16 pEmpty }; /// 17 plotID StringToPlotID(TString str); -void Check_Simulation(TString filename = "transfer.root", +void Check_Simulation(TString filename = "transfer1.root", TString configFile = "../working/Check_Simulation_Config.txt", Int_t padSize = 500, bool outputCanvas = false){ @@ -139,25 +139,33 @@ void Check_Simulation(TString filename = "transfer.root", AnalysisLib::DetGeo detGeo = AnalysisLib::LoadDetectorGeo(detGeoTxt); + AnalysisLib::Array array; + if( detGeo.use2ndArray){ + array = detGeo.array2; + }else{ + array = detGeo.array1; + } + double field = detGeo.Bfield; TString fdmsg = field > 0 ? "out of plan" : "into plan"; TString msg2; msg2.Form("field = %.2f T, %s", field, fdmsg.Data()); - double prepDist = detGeo.array1.detPerpDist; - double length = detGeo.array1.detLength; - double posRecoil = detGeo.recoilPos; + double prepDist = array.detPerpDist; + double length = array.detLength; + + double posRecoil = detGeo.recoilPos; double rhoRecoilIn = detGeo.recoilInnerRadius; double rhoRecoilOut = detGeo.recoilOuterRadius; double posRecoil1 = detGeo.recoilPos1; double posRecoil2 = detGeo.recoilPos2; - vector pos = detGeo.array1.detPos; + vector pos = array.detPos; - float firstPos = detGeo.array1.firstPos; - int rDet = detGeo.array1.nDet; - int cDet = detGeo.array1.mDet; + float firstPos = array.firstPos; + int rDet = array.nDet; + int cDet = array.mDet; double elum1 = detGeo.elumPos1; @@ -176,8 +184,8 @@ void Check_Simulation(TString filename = "transfer.root", int numDet = rDet * cDet; - zRange[1] = detGeo.zMin - 50; - zRange[2] = detGeo.zMax + 50; + zRange[1] = array.zMin - 50; + zRange[2] = array.zMax + 50; printf(" zRange : %f - %f \n", zRange[1], zRange[2]); printf("=================================\n"); diff --git a/Cleopatra/HELIOS_LIB.h b/Cleopatra/HELIOS_LIB.h index c93aa3f..101e4ec 100644 --- a/Cleopatra/HELIOS_LIB.h +++ b/Cleopatra/HELIOS_LIB.h @@ -444,22 +444,22 @@ public: bool SetDetectorGeometry(string filename); void SetBeamPosition(double x, double y) { xOff = x; yOff = y;} - void OverrideMagneticField(double BField){ this->Bfield = BField; this->sign = BField > 0 ? 1: -1;} - void OverrideMagneticFieldDirection(double BfieldThetaInDeg){ this->BfieldTheta = BfieldThetaInDeg;} + void OverrideMagneticField(double BField){ this->detGeo.Bfield = BField; this->detGeo.BfieldSign = BField > 0 ? 1: -1;} + void OverrideMagneticFieldDirection(double BfieldThetaInDeg){ this->detGeo.BfieldTheta = BfieldThetaInDeg;} void OverrideFirstPos(double firstPos){ overrideFirstPos = true; printf("------ Overriding FirstPosition to : %8.2f mm \n", firstPos); - this->firstPos = firstPos; + this->array.firstPos = firstPos; } void OverrideDetectorDistance(double perpDist){ overrideDetDistance = true; printf("------ Overriding Detector Distance to : %8.2f mm \n", perpDist); - this->perpDist = perpDist; + this->array.detPerpDist = perpDist; } void SetDetectorOutside(bool isOutside){ - this->isFromOutSide = isOutside; - printf(" Detectors are facing %s\n", isFromOutSide ? "outside": "inside" ); + this->array.detFaceOut = isOutside; + printf(" Detectors are facing %s\n", array.detFaceOut ? "outside": "inside" ); } int DetAcceptance(); @@ -469,7 +469,7 @@ public: void CalTrajectoryPara(TLorentzVector P, int Z, int id); // id = 0 for Pb, id = 1 for PB. - int GetNumberOfDetectorsInSamePos(){return mDet;} + int GetNumberOfDetectorsInSamePos(){return array.mDet;} double GetEnergy(){return e;} double GetDetX(){return detX;} // position in each detector, range from -1, 1 @@ -486,17 +486,17 @@ public: return sqrt(x*x+y*y); } - double GetXPos(double ZPos){ return XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, sign); } - double GetYPos(double ZPos){ return YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, sign); } - double GetR(double ZPos) { return RPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, sign); } + double GetXPos(double ZPos){ return XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } + double GetYPos(double ZPos){ return YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } + double GetR(double ZPos) { return RPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } double GetRecoilEnergy(){return eB;} - double GetRecoilXPos(double ZPos){ return XPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, sign); } - double GetRecoilYPos(double ZPos){ return YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, sign); } - double GetRecoilR(double ZPos) { return RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, sign); } + double GetRecoilXPos(double ZPos){ return XPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } + double GetRecoilYPos(double ZPos){ return YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } + double GetRecoilR(double ZPos) { return RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } - double GetBField() {return Bfield;} - double GetDetRadius() {return perpDist;} + double GetBField() {return detGeo.Bfield;} + double GetDetRadius() {return array.detPerpDist;} trajectory GetTrajectory_b() {return orbitb;} trajectory GetTrajectory_B() {return orbitB;} @@ -523,6 +523,7 @@ private: } AnalysisLib::DetGeo detGeo; + AnalysisLib::Array array; trajectory orbitb, orbitB; @@ -536,22 +537,21 @@ private: double xOff, yOff; // beam position //detetcor Geometry - double Bfield; // T - int sign ; // sign of B-field - double BfieldTheta; // rad, 0 = z-axis, pi/2 = y axis, pi = -z axis - double bore; // bore , mm - double perpDist; // distance from axis - double width; // width - double posRecoil; // recoil, downstream - double rhoRecoilin; // radius recoil inner - double rhoRecoilout; // radius recoil outter - double length; // length - double blocker; - double firstPos; // m - vector pos; // near position in m - int nDet, mDet; // nDet = number of different pos, mDet, number of same pos - - bool isFromOutSide; + //int sign ; // sign of B-field + // double Bfield; // T + // double BfieldTheta; // rad, 0 = z-axis, pi/2 = y axis, pi = -z axis + // double bore; // bore , mm + // double perpDist; // distance from axis + // double width; // width + // double posRecoil; // recoil, downstream + // double rhoRecoilin; // radius recoil inner + // double rhoRecoilout; // radius recoil outter + // double length; // length + // double blocker; + // double firstPos; // m + // vector pos; // near position in m + // int nDet, mDet; // nDet = number of different pos, mDet, number of same pos + // bool isFromOutSide; bool overrideDetDistance; bool overrideFirstPos; @@ -575,24 +575,6 @@ HELIOS::HELIOS(){ isDetReady = false; - Bfield = 0; - BfieldTheta = 0; - sign = 1; - bore = 1000; - perpDist = 10; - width = 10; - posRecoil = 0; - rhoRecoilin = 0; - rhoRecoilout = 0; - length = 0; - blocker = 0; - firstPos = 0; - pos.clear(); - nDet = 0; - mDet = 0; - - isFromOutSide = true; //default is facing outside - overrideDetDistance = false; overrideFirstPos = false; isCoincidentWithRecoil = false; @@ -609,25 +591,13 @@ bool HELIOS::SetDetectorGeometry(string filename){ if( haha->ReadFile(filename.c_str()) > 0 ) { detGeo = AnalysisLib::LoadDetectorGeo(haha); - PrintDetGeo(detGeo); + if( detGeo.use2ndArray ){ + array = detGeo.array2; + }else{ + array = detGeo.array1; + } - Bfield = detGeo.Bfield; - BfieldTheta = detGeo.BfieldTheta; - sign = detGeo.BfieldSign; - bore = detGeo.bore; - posRecoil = detGeo.recoilPos; - rhoRecoilin = detGeo.recoilInnerRadius; - rhoRecoilout = detGeo.recoilOuterRadius; - - length = detGeo.array1.detLength; - width = detGeo.array1.detWidth; - perpDist = detGeo.array1.detPerpDist; - blocker = detGeo.array1.blocker; - firstPos = detGeo.array1.firstPos; - pos = detGeo.array1.detPos; - nDet = detGeo.array1.nDet; - mDet = detGeo.array1.mDet; - isFromOutSide = detGeo.array1.detFaceOut; + PrintDetGeo(detGeo, false); isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil; @@ -647,45 +617,45 @@ int HELIOS::DetAcceptance(){ if( isDetReady == false ) return 0; // -1 ========= when recoil direction is not same side of array - if( firstPos < 0 && orbitb.z > 0 ) return -1; - if( firstPos > 0 && orbitb.z < 0 ) return -1; + if( array.firstPos < 0 && orbitb.z > 0 ) return -1; + if( array.firstPos > 0 && orbitb.z < 0 ) return -1; // -11 ======== rho is too small - if( 2 * orbitb.rho < perpDist ) return -11; + if( 2 * orbitb.rho < array.detPerpDist ) return -11; // -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist) if( orbitb.detRowID == -1 ) return -15; // -10 =========== when rho is too big . - if( bore < 2 * orbitb.rho) return -10; + if( detGeo.bore < 2 * orbitb.rho) return -10; // -14 ========== check particle-B hit radius on recoil dectector - if( isCoincidentWithRecoil && orbitB.R > rhoRecoilout ) return -14; + if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) return -14; //if( isCoincidentWithRecoil && (orbitB.R > rhoRecoilout || orbitB.R < rhoRecoilin) ) return -14; // -12 ========= check is particle-b was blocked by recoil detector - rhoHit = GetR(posRecoil); - if( orbitb.z > 0 && posRecoil > 0 && orbitb.z > posRecoil && rhoHit < rhoRecoilout ) return -12; - if( orbitb.z < 0 && posRecoil < 0 && orbitb.z < posRecoil && rhoHit < rhoRecoilout ) return -12; + rhoHit = GetR(detGeo.recoilPos); + if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) return -12; + if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) return -12; // -13 ========= not more than 3 loops if( orbitb.loop > 3 ) return -13; // -2 ========= calculate the "y"-distance from detector center - if( sqrt(orbitb.R*orbitb.R - perpDist*perpDist)> width/2 ) return -2; + if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) return -2; // -3 ==== when zPos further the range of whole array, more loop would not save - if( firstPos < 0 && orbitb.z < pos[0] - length ) return -3; - if( firstPos > 0 && orbitb.z > pos[nDet-1] + length ) return -3; + if( array.firstPos < 0 && orbitb.z < array.detPos[0] - array.detLength ) return -3; + if( array.firstPos > 0 && orbitb.z > array.detPos[array.nDet-1] + array.detLength ) return -3; // -4 ======== Hit on blacker - if( blocker != 0 && firstPos > 0 && pos[0] - blocker < orbitb.z && orbitb.z < pos[0] ) return -4; - if( blocker != 0 && firstPos < 0 && pos[nDet-1] < orbitb.z && orbitb.z < pos[nDet-1] + blocker ) return -4; + if( array.blocker != 0 && array.firstPos > 0 && array.detPos[0] - array.blocker < orbitb.z && orbitb.z < array.detPos[0] ) return -4; + if( array.blocker != 0 && array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z && orbitb.z < array.detPos[array.nDet-1] + array.blocker ) return -4; // 2 ====== when zPos less then the nearest position, more loop may hit int increaseLoopFlag = 0; - if( firstPos < 0 && pos[nDet-1] < orbitb.z ) increaseLoopFlag = 2; - if( firstPos > 0 && pos[0] > orbitb.z ) increaseLoopFlag = 2; + if( array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z ) increaseLoopFlag = 2; + if( array.firstPos > 0 && array.detPos[0] > orbitb.z ) increaseLoopFlag = 2; if (increaseLoopFlag == 2 ) { orbitb.z += orbitb.z0; orbitb.effLoop += 1.0; @@ -695,33 +665,34 @@ int HELIOS::DetAcceptance(){ } // 1 ======= check hit array z- position - if( firstPos < 0 ){ - for( int i = 0; i < nDet; i++){ - if( pos[i]-length <= orbitb.z && orbitb.z <= pos[i]) { + if( array.firstPos < 0 ){ + for( int i = 0; i < array.nDet; i++){ + if( array.detPos[i] - array.detLength <= orbitb.z && orbitb.z <= array.detPos[i]) { orbitb.detID = i; - detX = ( orbitb.z - (pos[i] + length/2 ))/ length*2 ;// range from -1 , 1 + detX = ( orbitb.z - (array.detPos[i] + array.detLength/2 ))/ array.detLength * 2 ;// range from -1 , 1 return 1; } } }else{ - for( int i = 0; i < nDet ; i++){ - if( pos[i] <= orbitb.z && orbitb.z <= pos[i] + length) { - ///printf(" %d | %f < z = %f < %f \n", i, pos[i], orbitb.z, pos[i]+length); + for( int i = 0; i < array.nDet ; i++){ + if( array.detPos[i] <= orbitb.z && orbitb.z <= array.detPos[i] + array.detLength) { + ///printf(" %d | %f < z = %f < %f \n", i, array.detPos[i], orbitb.z, array.detPos[i]+length); orbitb.detID = i; - detX = ( orbitb.z - (pos[i] - length/2 ))/ length*2 ;// range from -1 , 1 + detX = ( orbitb.z - (array.detPos[i] - array.detLength/2 ))/ array.detLength*2 ;// range from -1 , 1 return 1; } } } + // -5 ======== check hit array gap - if( firstPos < 0 ){ - for( int i = 0; i < nDet-1 ; i++){ - if( pos[i] < orbitb.z && orbitb.z < pos[i+1] - length ) return -5; //increaseLoopFlag = 3; + if( array.firstPos < 0 ){ + for( int i = 0; i < array.nDet-1 ; i++){ + if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) return -5; //increaseLoopFlag = 3; } }else{ - for( int i = 0; i < nDet-1 ; i++){ - if( pos[i] + length < orbitb.z && orbitb.z < pos[i+1] ) return -5; //increaseLoopFlag = 3; + for( int i = 0; i < array.nDet-1 ; i++){ + if( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) return -5; //increaseLoopFlag = 3; } } if (increaseLoopFlag == 3 ) { @@ -741,7 +712,7 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){ if( id == 0 ){ orbitb.theta = P.Theta(); orbitb.phi = P.Phi(); - orbitb.rho = P.Pt() / abs(Bfield) / Z / c * 1000; //mm + orbitb.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm orbitb.vt = P.Beta() * TMath::Sin(P.Theta()) * c ; // mm / nano-second orbitb.vp = P.Beta() * TMath::Cos(P.Theta()) * c ; // mm / nano-second orbitb.t0 = TMath::TwoPi() * orbitb.rho / orbitb.vt; // nano-second @@ -755,7 +726,7 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){ if( id == 1 ){ orbitB.theta = P.Theta(); orbitB.phi = P.Phi(); - orbitB.rho = P.Pt() / abs(Bfield) / Z / c * 1000; //mm + orbitB.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm orbitB.vt = P.Beta() * TMath::Sin(P.Theta()) * c ; // mm / nano-second orbitB.vp = P.Beta() * TMath::Cos(P.Theta()) * c ; // mm / nano-second orbitB.t0 = TMath::TwoPi() * orbitB.rho / orbitB.vt; // nano-second @@ -776,7 +747,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ CalTrajectoryPara(Pb, Zb, 0); int targetLoop = 1; - int inOut = isFromOutSide == true ? 1: 0; //1 = from Outside, 0 = from inside + int inOut = array.detFaceOut == true ? 1: 0; //1 = from Outside, 0 = from inside bool debug = false; @@ -792,18 +763,18 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ vector zPossible; vector dID; //detRowID - int iStart = ( sign == 1 ? 0 : -mDet ); - int iEnd = ( sign == 1 ? 2*mDet : mDet ); + int iStart = ( detGeo.BfieldSign == 1 ? 0 : -array.mDet ); + int iEnd = ( detGeo.BfieldSign == 1 ? 2 * array.mDet : array.mDet ); for( int i = iStart; i < iEnd ; i++){ - double phiD = TMath::TwoPi()/mDet * i ; + double phiD = TMath::TwoPi()/array.mDet * i ; double dphi = orbitb.phi - phiD; - double aEff = perpDist - (xOff * TMath::Cos(phiD) + yOff * TMath::Sin(phiD)); - double hahaha = asin( aEff/ orbitb.rho - sign * sin(dphi)); + double aEff = array.detPerpDist - (xOff * TMath::Cos(phiD) + yOff * TMath::Sin(phiD)); + double hahaha = asin( aEff/ orbitb.rho - detGeo.BfieldTheta * sin(dphi)); int n = 2*targetLoop + inOut; - double zP = orbitb.z0 /TMath::TwoPi() * ( sign * dphi + n * TMath::Pi() + pow(-1, n) * hahaha ); + double zP = orbitb.z0 /TMath::TwoPi() * ( detGeo.BfieldSign * dphi + n * TMath::Pi() + pow(-1, n) * hahaha ); if( debug ) { double xP = GetXPos(zP) ; @@ -848,15 +819,15 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ orbitb.detRowID = (12+dID[i])%4; orbitb.t = orbitb.t0 * orbitb.effLoop; - double phiD = TMath::TwoPi()/mDet * dID[i] ; + double phiD = TMath::TwoPi()/array.mDet * dID[i] ; double dphi = orbitb.phi - phiD ; if( debug ) { // Check is in or out - double hitDir = cos( orbitb.z/orbitb.z0 * TMath::TwoPi() - sign * dphi ); + double hitDir = cos( orbitb.z/orbitb.z0 * TMath::TwoPi() - detGeo.BfieldSign * dphi ); printf(" hitDir : %4.1f ", hitDir); if( ( inOut == 1 && hitDir > 0 ) || (inOut == 0 && hitDir < 0 ) ) { - printf(" != %f ", perpDist); + printf(" != %f ", array.detPerpDist); orbitb.z = TMath::QuietNaN(); orbitb.loop = -1; orbitb.detRowID = -1; @@ -868,8 +839,8 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ double yPos = GetYPos(orbitb.z ) ; double a = xPos * cos(phiD) + yPos * sin(phiD); printf(" a : %f ", a); - if( abs(a - perpDist) > 0.01) { - printf(" != %f ", perpDist); + if( abs(a - array.detPerpDist) > 0.01) { + printf(" != %f ", array.detPerpDist); orbitb.z = TMath::QuietNaN(); orbitb.loop = -1; orbitb.detRowID = -1; @@ -892,244 +863,16 @@ int HELIOS::CalRecoilHit(TLorentzVector PB, int ZB){ CalTrajectoryPara(PB, ZB, 1); - orbitB.z = posRecoil; - orbitB.x = GetRecoilXPos(posRecoil) ; - orbitB.y = GetRecoilYPos(posRecoil) ; - orbitB.R = GetRecoilR(posRecoil); + orbitB.z = detGeo.recoilPos; + orbitB.x = GetRecoilXPos(detGeo.recoilPos) ; + orbitB.y = GetRecoilYPos(detGeo.recoilPos) ; + orbitB.R = GetRecoilR(detGeo.recoilPos); orbitB.effLoop = orbitB.z/orbitB.z0; orbitB.t = orbitB.t0 * orbitB.effLoop ; return 1; } -/* -int HELIOS::CalHit(TLorentzVector Pb, int Zb, TLorentzVector PB, int ZB, double xOff, double yOff){ - - //initialization - int hit = 0; - const double c = 299.792458; // mm/ns, standard speed for MeV conversion - theta = TMath::QuietNaN(); - e = TMath::QuietNaN(); - z = TMath::QuietNaN(); - detX = TMath::QuietNaN(); - t = TMath::QuietNaN(); - rho = TMath::QuietNaN(); - dphi = TMath::QuietNaN(); - detID = -1; - loop = -1; - detRowID = -1; - - eB = TMath::QuietNaN(); - zB = TMath::QuietNaN(); - tB = TMath::QuietNaN(); - rhoB = TMath::QuietNaN(); - rxHit = TMath::QuietNaN(); - ryHit = TMath::QuietNaN(); - phiB = TMath::QuietNaN(); - - //range of detector azimuth angle, for beam at center - double azimu = TMath::Pi()/ mDet; - double azimuDet = TMath::ATan2(width/2., perpDist); - - //rotate Pb and PB to B-Field - //Pb.RotateX(BfieldTheta); - //PB.RotateX(BfieldTheta); - - //====================== X-Y plane, light particle - rho = Pb.Pt() / abs(Bfield) / Zb / c * 1000; //mm - theta = Pb.Theta(); - phi = Pb.Phi(); - - //====================== recoil detector - thetaB = PB.Theta(); - phiB = PB.Phi(); - rhoB = PB.Pt() / abs(Bfield) / ZB / c * 1000; //mm - vt0B = PB.Beta() * TMath::Sin(thetaB) * c ; // mm / nano-second - vp0B = PB.Beta() * TMath::Cos(thetaB) * c ; // mm / nano-second - //tB = TMath::TwoPi() * rhoB / vt0B; // nano-second - tB = posRecoil / vp0B; // nano-second - eB = PB.E() - PB.M(); - zB = vp0B * tB; - rhoBHit = GetRecoilR(posRecoil); - rxHit = GetRecoilXPos(posRecoil); - ryHit = GetRecoilYPos(posRecoil); - - if( isDetReady == false ) { - //====================== infinite small detector - vt0 = Pb.Beta() * TMath::Sin(theta) * c ; // mm / nano-second - vp0 = Pb.Beta() * TMath::Cos(theta) * c ; // mm / nano-second - t = TMath::TwoPi() * rho / vt0; // nano-second - z = vp0 * t; // mm - e = Pb.E() - Pb.M(); - dphi = TMath::TwoPi(); - - return 0; - } - - if( bore > 2 * rho && 2*rho > perpDist && ((firstPos > 0 && theta < TMath::PiOver2()) || (firstPos < 0 && theta > TMath::PiOver2())) ){ - //====================== infinite small detector - vt0 = Pb.Beta() * TMath::Sin(theta) * c ; // mm / nano-second - vp0 = Pb.Beta() * TMath::Cos(theta) * c ; // mm / nano-second - t0 = TMath::TwoPi() * rho / vt0; // nano-second - z0 = vp0 * t0; // mm - - //========== check particle-B hit radius on recoil dectector - if(isCoincidentWithRecoil && rhoBHit > rhoRecoilout ) return -2; // when particle-B miss the recoil detector - - //========= check is particle-b was blocked by recoil detector - rhoHit = GetR(posRecoil) ;// radius of light particle b at recoil detector - if( z0 > 0 && posRecoil > 0 && z0 > posRecoil && rhoHit < rhoRecoilout) { - return -1 ; // when particle-b blocked by recoil detector - } - if( z0 < 0 && posRecoil < 0 && z0 < posRecoil && rhoHit < rhoRecoilout) { - return -1 ; - } - - //================ Calulate z hit - double zHit = TMath::QuietNaN(); - bool isHit = false; - bool isHitFromOutside = false; - bool isReachArrayCoverage = false; - - loop = 0; - int startJ = (int) fmod(TMath::Ceil(mDet*phi/TMath::TwoPi() - 0.5) ,mDet) ; - - //printf("======================= T : %5.2f, theta : %7.2f , phi : %7.2f \n", Pb.E() - Pb.M(), theta*TMath::RadToDeg(), phi*TMath::RadToDeg()); - - // loop until reach the detector position covrage. - do{ - loop += 1; - //int n = 2*loop + sign; - int n = 2*loop -1; - - if( blocker != 0.0 && abs(firstPos/blocker) < loop ) return -6; - - //======= check Block By blocker using z0, speed thing up - if( firstPos > 0 ){ - if( pos[0] - blocker < z0 * loop && z0 * loop < pos[0] ) return -6; // blocked by blocker - }else{ - if( pos[nDet-1] < z0 * loop && z0 * loop < pos[nDet-1] + blocker ) return -6; - } - - if( loop > 10 ) { - return -3; // when loop > 10 - break; // maximum 10 loops - } - - for( int j = startJ ; j < startJ + mDet; j++){ - - double phiDet = TMath::TwoPi() / mDet * (j); // detector plane angle - isReachArrayCoverage = false; - isHitFromOutside = false; - - //========== calculate zHit by solving the cycltron orbit and the detector planes. - double aEff = perpDist - (xOff * TMath::Cos(phiDet) + yOff * TMath::Sin(phiDet)); - double dphi = phi - phiDet; - double z0 = TMath::TwoPi() * rho / TMath::Tan(theta); // the cycle - - //with sign is the correct formula, but somehow, probably the ASin? give incorrect result for sign = 1; - //zHit = z0 / TMath::TwoPi() * ( sign * dphi + TMath::Power(-1, n) * TMath::ASin(aEff/rho - sign * TMath::Sin(dphi)) + TMath::Pi() * n ); - zHit = z0 / TMath::TwoPi() * ( -1 * dphi + TMath::Power(-1, n) * TMath::ASin(aEff/rho + TMath::Sin(dphi)) + TMath::Pi() * n ); - e = Pb.E() - Pb.M(); - z = zHit; - - //======= calculate the distance from middle of detector - double xHit = GetXPos(zHit) + xOff; // resotre the beam to be at the center - double yHit = GetYPos(zHit) + yOff; - double sHit = TMath::Sqrt(xHit*xHit + yHit*yHit - perpDist*perpDist); // is it hit on the detector - - ///if( zHit > z0 ) continue; - ///printf("==== %7.2f, %7.2f | n : %d, row : %2d, phiD : %4.0f, rho : %9.4f, z0 : %9.4f, zHit : %9.4f, xHit : %9.4f, yHit : %9.4f \n", - /// theta*TMath::RadToDeg(), phi*TMath::RadToDeg(), n, j, phiDet*TMath::RadToDeg(), rho, z0, zHit, xHit, yHit); - - - if( sHit > width /2.) continue; // if the sHit is large, it does not hit on detector, go to next mDet - - //======== this is the particel direction (normalized) dot normal vector of the detector plane - //double dir = TMath::Cos(zHit/z0 * TMath::TwoPi() - sign * dphi); - double dir = TMath::Cos(zHit/z0 * TMath::TwoPi() + dphi); - if( dir < 0) {// when dir == 0, no solution - isHitFromOutside = true; - }else{ - return -5 ; // hit from inside. - } - - //### after this, the particle is hit on detector, from outside. - - //======= check Block By blocker using z0 and zHit; z0 has to be inside the nearst array, but zHit is in the blocker. - if( firstPos > 0 && blocker > 0.0){ - if( pos[0] < z0*loop && pos[0] - blocker < zHit && zHit < pos[0] ) return -6; - }else{ - if( z0*loop < pos[nDet-1] && pos[nDet-1] < zHit && zHit < pos[nDet-1] + blocker ) return -6; - } - - //### after this, the particle is hit on detector, from outside, and not blocked by blocker - - //======= check within the detector range - if( firstPos > 0 ){ - if( zHit < pos[0] ) continue; // goto next mDet, after progress of all side, still not OK, then next loop - if( zHit > pos[nDet-1] + length) return -4; // since the zHit is mono-increse, when zHit shoot over the detector - }else{ - if( zHit < pos[0] - length ) return -4; - if( zHit > pos[nDet-1]) continue; - } - - //====== check hit - if( !isReachArrayCoverage && isHitFromOutside && sHit < width/2.){ - isHit = true; - isReachArrayCoverage = true; - detRowID = (j+mDet) % mDet; - break; // if isHit, break, don't calculate for the rest of the detector - }else{ - isReachArrayCoverage = false; - } - } - }while(!isReachArrayCoverage); - - if( !isHit ) return -7; // zHit in the gap of detector - - //===== final calculation for light particle - e = Pb.E() - Pb.M(); - z = zHit; - t = zHit / vp0; - dphi = t * vt0 / rho; - - //sometimes, dphi = 2pi + something, i.e. loop a bit more than a single loop. - if( dphi / TMath::TwoPi() > loop ) return -8; - double xHit = GetXPos(zHit) + xOff; // resotre the beam to be at the center - double yHit = GetYPos(zHit) + yOff; - //======= Check inside the detector - double eta = TMath::Abs(TMath::ATan2(yHit,xHit)); - double etaMod = TMath::Abs(fmod(eta + azimu, 2* azimu) - azimu); - if ( etaMod > azimuDet ) return -8; - - //=========== check hit on detector gap - for( int i = 0 ; i < nDet ; i++){ - if( firstPos > 0 ){ - if( pos[i] < z && z < pos[i] + length ){ - detID = i; - detX = ( z - (pos[i] + length/2 ))/ length*2 ;// range from -1 , 1 - } - }else{ - if( pos[i] - length < z && z < pos[i] ){ - detID = i; - detX = ( z - (pos[i] - length/2 ))/ length*2 ;// range from -1 , 1 - } - } - } - - if( detID >=0 ){ - hit = 1; - }else{ - hit = -9; // particle-b hit the gap - } - }else{ - hit = -10; // hit helio wall or (detector upstream, particle go downstream, via versa) - } - - return hit; -} -*/ //======================================================= //####################################################### // Class for multi-scattering of the beam inside target diff --git a/Cleopatra/Simulation_Helper.C b/Cleopatra/Simulation_Helper.C index 78775f4..41b20b7 100644 --- a/Cleopatra/Simulation_Helper.C +++ b/Cleopatra/Simulation_Helper.C @@ -17,6 +17,7 @@ #include "../Cleopatra/ExtractXSec.h" #include "../Cleopatra/PlotTGraphTObjArray.h" #include "../armory/AutoFit.C" +#include "../armory/AnalysisLib.h" #include "../Cleopatra/Check_Simulation.C" #include @@ -60,7 +61,8 @@ private: TGTextEntry * txtName ; TGTextEntry * txtEx ; - + + bool isUse2ndArray; public: MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h); @@ -69,6 +71,8 @@ public: void OpenFile(int); void GetData(); bool IsFileExist(TString filename); + + void CheckIsUse2ndArray(); }; @@ -85,7 +89,7 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) { TGVerticalFrame *hframe2 = new TGVerticalFrame(fMain,600,800 ); hframe->AddFrame(hframe2,new TGLayoutHints( kLHintsExpandX | kLHintsExpandY, 2,2,2,2)); - fileName = "../working/reactionConfig.txt"; + fileName = "../working/detectorGeo.txt"; TGHorizontalFrame *hframe00 = new TGHorizontalFrame(hframe2,600,600 ); hframe2->AddFrame(hframe00, new TGLayoutHints(kLHintsCenterX | kLHintsExpandX , 2,2,2,2)); @@ -127,13 +131,6 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) { {//================= Simulation group TGGroupFrame * simFrame = new TGGroupFrame(hframe1, "Kinematics Simulation", kVerticalFrame); hframe1->AddFrame(simFrame, new TGLayoutHints(kLHintsCenterX, 5,5,3,4)); - - TGTextButton *openRec = new TGTextButton(simFrame, "reaction Config"); - openRec->SetWidth(150); - openRec->SetHeight(20); - openRec->ChangeOptions( openRec->GetOptions() | kFixedSize ); - openRec->Connect("Clicked()","MyMainFrame",this, "OpenFile(=1)"); - simFrame->AddFrame(openRec,new TGLayoutHints(kLHintsRight, 5,5,3,4)); TGTextButton *openDet = new TGTextButton(simFrame, "detector Geo."); openDet->SetWidth(150); @@ -142,6 +139,13 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) { openDet->Connect("Clicked()","MyMainFrame",this, "OpenFile(=0)"); simFrame->AddFrame(openDet,new TGLayoutHints(kLHintsRight, 5,5,3,4)); + TGTextButton *openRec = new TGTextButton(simFrame, "reaction Config"); + openRec->SetWidth(150); + openRec->SetHeight(20); + openRec->ChangeOptions( openRec->GetOptions() | kFixedSize ); + openRec->Connect("Clicked()","MyMainFrame",this, "OpenFile(=1)"); + simFrame->AddFrame(openRec,new TGLayoutHints(kLHintsRight, 5,5,3,4)); + TGTextButton *openEx = new TGTextButton(simFrame, "Ex List"); openEx->SetWidth(150); openEx->SetHeight(20); @@ -168,16 +172,16 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) { openSimChk->Connect("Clicked()","MyMainFrame",this, "OpenFile(=4)"); simFrame->AddFrame(openSimChk,new TGLayoutHints(kLHintsRight, 5,5,3,4)); - TGTextButton *SimChk = new TGTextButton(simFrame,"Plot Simulation"); + TGTextButton *SimChk = new TGTextButton(simFrame,"Re-Plot Simulation"); SimChk->SetWidth(150); - SimChk->SetHeight(40); + //SimChk->SetHeight(40); SimChk->ChangeOptions( SimChk->GetOptions() | kFixedSize ); SimChk->Connect("Clicked()","MyMainFrame",this,"Command(=2)"); simFrame->AddFrame(SimChk, new TGLayoutHints(kLHintsRight,5,5,3,4)); TGTextButton *autoFit = new TGTextButton(simFrame,"AutoFit ExCal"); autoFit->SetWidth(150); - autoFit->SetHeight(40); + //autoFit->SetHeight(40); autoFit->ChangeOptions( autoFit->GetOptions() | kFixedSize ); autoFit->Connect("Clicked()","MyMainFrame",this,"Command(=5)"); simFrame->AddFrame(autoFit, new TGLayoutHints(kLHintsRight,5,5,3,4)); @@ -352,20 +356,40 @@ bool MyMainFrame::IsFileExist(TString filename){ return file.is_open(); } +void MyMainFrame::CheckIsUse2ndArray(){ + + TMacro * haha = new TMacro("../working/detectorGeo.txt"); + AnalysisLib::DetGeo detGeo = AnalysisLib::LoadDetectorGeo(haha); + delete haha; + isUse2ndArray = detGeo.use2ndArray; + +} + void MyMainFrame::OpenFile(int ID){ editor->SaveFile(fileName); TString oldFileName = fileName; - + if ( ID == 0 ) fileName = "../working/detectorGeo.txt"; - if ( ID == 1 ) fileName = "../working/reactionConfig.txt"; - if ( ID == 2 ) fileName = "../working/Ex.txt"; - if ( ID == 3 ) fileName = "../working/DWBA"; + + CheckIsUse2ndArray(); + if( isUse2ndArray){ + if ( ID == 1 ) fileName = "../working/reactionConfig2.txt"; + if ( ID == 2 ) fileName = "../working/Ex2.txt"; + if ( ID == 3 ) fileName = "../working/DWBA2"; + if ( ID == 5 ) fileName = "../working/DWBA2.in"; + if ( ID == 6 ) fileName = "../working/DWBA2.out"; + if ( ID == 7 ) fileName = "../working/DWBA2.Xsec.txt"; + }else{ + if ( ID == 1 ) fileName = "../working/reactionConfig1.txt"; + if ( ID == 2 ) fileName = "../working/Ex1.txt"; + if ( ID == 3 ) fileName = "../working/DWBA1"; + if ( ID == 5 ) fileName = "../working/DWBA1.in"; + if ( ID == 6 ) fileName = "../working/DWBA1.out"; + if ( ID == 7 ) fileName = "../working/DWBA1.Xsec.txt"; + } if ( ID == 4 ) fileName = "../working/Check_Simulation_Config.txt"; - if ( ID == 5 ) fileName = "../working/DWBA.in"; - if ( ID == 6 ) fileName = "../working/DWBA.out"; - if ( ID == 7 ) fileName = "../working/DWBA.Xsec.txt"; if ( ID == 8 ) fileName = isoFileName; //test if file exist @@ -417,6 +441,8 @@ void MyMainFrame::GetData(){ void MyMainFrame::Command(int ID) { editor->SaveFile(fileName); + + CheckIsUse2ndArray(); if( ID == 0 ){ @@ -468,29 +494,54 @@ void MyMainFrame::Command(int ID) { } if( ID == 1 ){ - string basicConfig = "reactionConfig.txt"; + + string basicConfig = "reactionConfig1.txt"; string heliosDetGeoFile = "detectorGeo.txt"; - string excitationFile = "Ex.txt"; //when no file, only ground state + string excitationFile = "Ex1.txt"; //when no file, only ground state TString ptolemyRoot = ""; // when no file, use isotropic distribution of thetaCM - TString saveFileName = "transfer.root"; - TString filename = "reaction.dat"; //when no file, no output + TString saveFileName = "transfer1.root"; + TString filename = "reaction1.dat"; //when no file, no output if( withDWBA->GetState() ) { - ptolemyRoot = "DWBA.root"; - excitationFile = "DWBA.Ex.txt"; + ptolemyRoot = "DWBA1.root"; + excitationFile = "DWBA1.Ex.txt"; } + + if( isUse2ndArray ){ + basicConfig = "reactionConfig2.txt"; + heliosDetGeoFile = "detectorGeo.txt"; + excitationFile = "Ex2.txt"; //when no file, only ground state + ptolemyRoot = ""; // when no file, use isotropic distribution of thetaCM + saveFileName = "transfer2.root"; + filename = "reaction2.dat"; //when no file, no output + + if( withDWBA->GetState() ) { + ptolemyRoot = "DWBA2.root"; + excitationFile = "DWBA2.Ex.txt"; + } + } + statusLabel->SetText("Running simulation......."); Transfer( basicConfig, heliosDetGeoFile, excitationFile, ptolemyRoot, saveFileName, filename); statusLabel->SetText("Plotting simulation......."); - Check_Simulation(); + + if( isUse2ndArray ){ + Check_Simulation("transfer2.root"); + }else{ + Check_Simulation("transfer1.root"); + } statusLabel->SetText("Plotted Simulation result"); } if( ID == 2 ){ - Check_Simulation(); - statusLabel->SetText(" Run Simulation first."); + if( isUse2ndArray ){ + Check_Simulation("transfer2.root"); + }else{ + Check_Simulation("transfer1.root"); + } + statusLabel->SetText(" Run Simulation first."); } if( ID == 3 ){ diff --git a/armory/AnalysisLib.h b/armory/AnalysisLib.h index aaf10f5..6fdef81 100644 --- a/armory/AnalysisLib.h +++ b/armory/AnalysisLib.h @@ -47,7 +47,7 @@ std::vector SplitStr(std::string tempLine, std::string splitter, in }while(pos != std::string::npos ); return output; -} +}; struct Array{ @@ -117,7 +117,7 @@ struct DetGeo{ Array array1; //==================2nd array - bool is2ndArrayExist; + bool use2ndArray; Array array2; double zMin, zMax; /// range of detectors @@ -169,7 +169,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){ detGeo.array1.pos.clear(); detGeo.array2.pos.clear(); - detGeo.is2ndArrayExist = false; + detGeo.use2ndArray = false; int detFlag = 0; int detLine = 0; @@ -219,7 +219,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){ } if( detFlag == 2){ - if ( detLine == 0 ) detGeo.is2ndArrayExist = str[0] == "false" ? false: true; + if ( detLine == 0 ) detGeo.use2ndArray = str[0] == "true" ? true : false; if ( detLine == 1 ) detGeo.array2.detPerpDist = atof(str[0].c_str()); if ( detLine == 2 ) detGeo.array2.detWidth = atof(str[0].c_str()); if ( detLine == 3 ) detGeo.array2.detLength = atof(str[0].c_str()); @@ -240,7 +240,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){ detGeo.zMin = detGeo.array1.zMin; detGeo.zMax = detGeo.array1.zMax; - if( detGeo.is2ndArrayExist) { + if( detGeo.use2ndArray) { detGeo.array2.DeduceAbsolutePos(); detGeo.zMin = TMath::Min(detGeo.array1.zMin, detGeo.array2.zMin); @@ -250,7 +250,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){ return detGeo; } -void PrintDetGeo(DetGeo detGeo){ +void PrintDetGeo(DetGeo detGeo, bool printAll = true){ printf("=====================================================\n"); printf(" B-field: %8.2f T, Theta : %6.2f deg \n", detGeo.Bfield, detGeo.BfieldTheta); @@ -258,15 +258,24 @@ void PrintDetGeo(DetGeo detGeo){ printf(" +---- field angle != 0 is not supported!!! \n"); } printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", detGeo.recoilPos, detGeo.recoilInnerRadius, detGeo.recoilOuterRadius); - printf("------------------------------------- Detector Position \n"); - detGeo.array1.PrintArray(); + if( printAll ){ + printf("------------------------------------- Detector Position \n"); + detGeo.array1.PrintArray(); - if( detGeo.is2ndArrayExist){ - printf("--------------------------------- 2nd Detector Position \n"); - detGeo.array2.PrintArray(); + if( detGeo.use2ndArray){ + printf("--------------------------------- 2nd Detector Position \n"); + detGeo.array2.PrintArray(); + } + }else{ + if( detGeo.use2ndArray){ + printf("--------------------------------- 2nd Detector Position \n"); + detGeo.array2.PrintArray(); + }else{ + printf("------------------------------------- Detector Position \n"); + detGeo.array1.PrintArray(); + } } - if( detGeo.elumPos1 != 0 || detGeo.elumPos2 != 0 || detGeo.recoilPos1 != 0 || detGeo.recoilPos2 != 0){ printf("=================================== Auxillary/Imaginary Detectors\n"); } @@ -334,7 +343,7 @@ ReactionConfig LoadReactionConfig(TMacro * macro){ reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA; reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ; - + return reaction; } @@ -378,9 +387,12 @@ void PrintReactionConfig(ReactionConfig reaction){ } DetGeo detGeo; -ReactionConfig reactionConfig; +ReactionConfig reactionConfig1; +ReactionConfig reactionConfig2; -void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt", std::string reactionConfigFileName = "reactionConfig.txt"){ +void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt", + std::string reactionConfigFileName1 = "reactionConfig1.txt", + std::string reactionConfigFileName2 = "reactionConfig2.txt"){ printf("=====================================================\n"); printf(" loading detector geometery : %s.", detGeoFileName.c_str()); TMacro * haha = new TMacro(); @@ -391,21 +403,34 @@ void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.t }else{ printf("... fail\n"); } + delete haha; printf("=====================================================\n"); - printf(" loading reaction config : %s.", reactionConfigFileName.c_str()); + printf(" loading reaction1 config : %s.", reactionConfigFileName1.c_str()); TMacro * kaka = new TMacro(); - if( kaka->ReadFile(reactionConfigFileName.c_str()) > 0 ) { - reactionConfig = AnalysisLib::LoadReactionConfig(kaka); + if( kaka->ReadFile(reactionConfigFileName1.c_str()) > 0 ) { + reactionConfig1 = AnalysisLib::LoadReactionConfig(kaka); printf("..... done.\n"); - AnalysisLib::PrintReactionConfig(reactionConfig); + AnalysisLib::PrintReactionConfig(reactionConfig1); }else{ printf("..... fail\n"); } - - delete haha; delete kaka; + if( detGeo.use2ndArray){ + printf("=====================================================\n"); + printf(" loading reaction2 config : %s.", reactionConfigFileName2.c_str()); + TMacro * jaja = new TMacro(); + if( jaja->ReadFile(reactionConfigFileName2.c_str()) > 0 ) { + reactionConfig2 = AnalysisLib::LoadReactionConfig(kaka); + printf("..... done.\n"); + AnalysisLib::PrintReactionConfig(reactionConfig2); + }else{ + printf("..... fail\n"); + } + delete jaja; + } + } //************************************** Correction parameters; @@ -502,11 +527,28 @@ void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.d if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]); } -double q, alpha, Et, betRel, gamm, G, massB, mass; //variables for Ex calculation -bool hasReactionPara = false; + +struct ReactionParas{ + + double Et; // total energy in CM frame + double beta; // Lorentz beta from Lab to CM + double gamma; // Lorentz gamma from Lab to CM + double alpha; // E-Z slope / beta + double G; //The G-coefficient.... + double massB; // heavy mass + double q; // charge of light particle + double mass; //light mass + bool hasReactionPara; + + double detPrepDist; + +}; + +ReactionParas reactParas1; +ReactionParas reactParas2; //~========================================= reaction parameters -void LoadReactionParasForArray1(bool verbose = false){ +void LoadReactionParas(int arrayID, bool verbose = false){ //check is the transfer.root is using the latest reactionConfig.txt //sicne reaction.dat is generated as a by-product of transfer.root @@ -527,59 +569,86 @@ void LoadReactionParasForArray1(bool verbose = false){ // system("../Cleopatra/Transfer"); // printf("########################## transfer.root updated\n"); //} + + ReactionParas * reactParas = nullptr; + + std::string fileName; + + if( arrayID == 1){ + reactParas = &AnalysisLib::reactParas1; + fileName = "reaction.dat"; + }else if( arrayID == 2){ + reactParas = &AnalysisLib::reactParas2; + fileName = "reaction2.dat"; + }else{ + printf("arrayID must be either 1 or 2. Abort.\n"); + return; + } + reactParas->detPrepDist = AnalysisLib::detGeo.array1.detPerpDist; + printf(" loading reaction parameters"); std::ifstream file; - file.open("reaction.dat"); - hasReactionPara = false; + file.open(fileName.c_str()); + reactParas->hasReactionPara = false; if( file.is_open() ){ std::string x; int i = 0; while( file >> x ){ if( x.substr(0,2) == "//" ) continue; - if( i == 0 ) mass = atof(x.c_str()); - if( i == 1 ) q = atof(x.c_str()); - if( i == 2 ) betRel = atof(x.c_str()); - if( i == 3 ) Et = atof(x.c_str()); - if( i == 4 ) massB = atof(x.c_str()); + if( i == 0 ) reactParas->mass = atof(x.c_str()); + if( i == 1 ) reactParas->q = atof(x.c_str()); + if( i == 2 ) reactParas->beta = atof(x.c_str()); + if( i == 3 ) reactParas->Et = atof(x.c_str()); + if( i == 4 ) reactParas->massB = atof(x.c_str()); i = i + 1; } printf("........ done.\n"); - hasReactionPara = true; - alpha = 299.792458 * abs(detGeo.Bfield) * q / TMath::TwoPi()/1000.; //MeV/mm - gamm = 1./TMath::Sqrt(1-betRel*betRel); - G = alpha * gamm * betRel * detGeo.array1.detPerpDist ; + reactParas->hasReactionPara = true; + reactParas->alpha = 299.792458 * abs(detGeo.Bfield) * reactParas->q / TMath::TwoPi()/1000.; //MeV/mm + reactParas->gamma = 1./TMath::Sqrt(1-reactParas->beta * reactParas->beta); + reactParas->G = reactParas->alpha * reactParas->gamma * reactParas->beta * reactParas->detPrepDist ; if( verbose ){ - printf("\tmass-b : %f MeV/c2 \n", mass); - printf("\tcharge-b : %f \n", q); - printf("\tE-total : %f MeV \n", Et); - printf("\tmass-B : %f MeV/c2 \n", massB); - printf("\tbeta : %f \n", betRel); + printf("\tmass-b : %f MeV/c2 \n", reactParas->mass); + printf("\tcharge-b : %f \n", reactParas->q); + printf("\tE-total : %f MeV \n", reactParas->Et); + printf("\tmass-B : %f MeV/c2 \n", reactParas->massB); + printf("\tbeta : %f \n", reactParas->beta); printf("\tB-field : %f T \n", detGeo.Bfield); - printf("\tslope : %f MeV/mm \n", alpha * betRel); - printf("\tdet radius: %f mm \n", detGeo.array1.detPerpDist); - printf("\tG-coeff : %f MeV \n", G); + printf("\tslope : %f MeV/mm \n", reactParas->alpha * reactParas->beta); + printf("\tdet radius: %f mm \n", reactParas->detPrepDist); + printf("\tG-coeff : %f MeV \n", reactParas->G); printf("=====================================================\n"); } }else{ printf("........ fail.\n"); - hasReactionPara = false; } file.close(); } std::vector CalExTheta(double e, double z){ - if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()}; + + ReactionParas * reactParas = nullptr; + + if( detGeo.array1.zMin <= z && z <= detGeo.array1.zMax ){ + reactParas = &reactParas1; + if( !reactParas->hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()}; + } + + if( detGeo.array2.zMin <= z && z <= detGeo.array2.zMax ){ + reactParas = &reactParas2; + if( !reactParas->hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()}; + } double Ex = TMath::QuietNaN(); double thetaCM = TMath::QuietNaN(); - double y = e + mass; // to give the KE + mass of proton; - double Z = alpha * gamm * betRel * z; - double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ; + double y = e + reactParas->mass; // to give the KE + mass of proton; + double Z = reactParas->alpha * reactParas->gamma * reactParas->beta * z; + double H = TMath::Sqrt(TMath::Power(reactParas->gamma * reactParas->beta,2) * (y*y - reactParas->mass * reactParas->mass) ) ; if( TMath::Abs(Z) < H ) { ///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi) @@ -590,26 +659,26 @@ std::vector CalExTheta(double e, double z){ int iter = 0; do{ phi = nPhi; - nPhi = phi - (H * TMath::Sin(phi) - G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - G /TMath::Power( TMath::Cos(phi), 2)); + nPhi = phi - (H * TMath::Sin(phi) - reactParas->G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - reactParas->G /TMath::Power( TMath::Cos(phi), 2)); iter ++; if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break; }while( TMath::Abs(phi - nPhi ) > tolerrence); phi = nPhi; /// check f'(phi) > 0 - double Df = H * TMath::Cos(phi) - G / TMath::Power( TMath::Cos(phi),2); + double Df = H * TMath::Cos(phi) - reactParas->G / TMath::Power( TMath::Cos(phi),2); if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){ double K = H * TMath::Sin(phi); - double x = TMath::ACos( mass / ( y * gamm - K)); - double momt = mass * TMath::Tan(x); /// momentum of particel b or B in CM frame - double EB = TMath::Sqrt(mass*mass + Et*Et - 2*Et*TMath::Sqrt(momt*momt + mass * mass)); - Ex = EB - massB; + double x = TMath::ACos( reactParas->mass / ( y * reactParas->gamma - K)); + double momt = reactParas->mass * TMath::Tan(x); /// momentum of particel b or B in CM frame + double EB = TMath::Sqrt(reactParas->mass * reactParas->mass + reactParas->Et * reactParas->Et - 2 * reactParas->Et * TMath::Sqrt(momt*momt + reactParas->mass * reactParas->mass)); + Ex = EB - reactParas->massB; - double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y; - double hahaha2 = gamm* betRel * momt; + double hahaha1 = reactParas->gamma * TMath::Sqrt(reactParas->mass * reactParas->mass + momt * momt) - y; + double hahaha2 = reactParas->gamma * reactParas->beta * momt; thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg(); - } + } } return {Ex, thetaCM}; diff --git a/working/DWBA b/working/DWBA2 similarity index 100% rename from working/DWBA rename to working/DWBA2 diff --git a/working/Ex.txt b/working/Ex1.txt similarity index 100% rename from working/Ex.txt rename to working/Ex1.txt diff --git a/working/Ex2.txt b/working/Ex2.txt new file mode 100644 index 0000000..160e6c1 --- /dev/null +++ b/working/Ex2.txt @@ -0,0 +1,6 @@ +//Ex relative_xsec SF sigma_in_MeV +//<--- use "//" for line comment +0.000 1.0 1.0 0.0100 +//4.400 1.0 1.0 0.0100 +//4.600 1.0 1.0 0.0100 +#============_End_of_file diff --git a/working/Monitor.C b/working/Monitor.C index 5c88ad3..a7d4b62 100644 --- a/working/Monitor.C +++ b/working/Monitor.C @@ -183,7 +183,8 @@ void Monitor::Begin(TTree *tree){ AnalysisLib::LoadXScaleCorr(); AnalysisLib::LoadECorr(); AnalysisLib::LoadRDTCorr(); - AnalysisLib::LoadReactionParasForArray1(true); + AnalysisLib::LoadReactionParas(1, true); + if( AnalysisLib::detGeo.use2ndArray ) AnalysisLib::LoadReactionParas(2, true); if( (int) AnalysisLib::xnCorr.size() < mapping::NARRAY ) { isXNCorrOK = false; printf(" !!!!!!!! size of xnCorr < NARRAY .\n"); } if( (int) AnalysisLib::xfxneCorr.size() < mapping::NARRAY ) { isXFXNCorrOK = false; printf(" !!!!!!!! size of xfxneCorr < NARRAY .\n"); } diff --git a/working/detectorGeo.txt b/working/detectorGeo.txt index d302ecc..d4453f6 100644 --- a/working/detectorGeo.txt +++ b/working/detectorGeo.txt @@ -26,7 +26,7 @@ Out //detector_facing_Out_or_In 235.8 //5th_det 290.0 #===============2nd_Array -true //is_2nd_detctor_exist_false=false_other=true +false //is_2nd_array_exist_is_use_for_Simulation 11.5 //distance_from_axis_[mm] 10.0 //width_of_detector_[mm] 50 //length_of_detector_[mm] diff --git a/working/reactionConfig.txt b/working/reactionConfig1.txt similarity index 72% rename from working/reactionConfig.txt rename to working/reactionConfig1.txt index e21ef48..5cee6fe 100644 --- a/working/reactionConfig.txt +++ b/working/reactionConfig1.txt @@ -1,16 +1,16 @@ -32 //beam_A -14 //beam_Z +32 //beam_A +14 //beam_Z 2 //target_A 1 //target_Z 1 //recoil_light_A 1 //recoil-light_Z -8.8 //beam-energy_in_MeV/u +8.8 //beam-energy_in_MeV/u 0.000 //beam-energy_sigma_in_MeV/u 0.000 //beam-angle_in_mrad 0.000 //beam-emittance_in_mrad 0.00 //x_offset_of_Beam_in_mm 0.00 //y_offset_of_Beam_in_mm -100000 //number_of_Event_being_generated +10000 //number_of_Event_being_generated false //isTargetScattering 0.913 //target_density_in_g/cm3 2.2e-4 //targetThickness_in_cm @@ -18,8 +18,8 @@ false //isTargetScattering ../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil ../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil false //isDacay -32 //decayNucleus_A -14 //decayNucleus_Z +32 //decayNucleus_A +14 //decayNucleus_Z false //isReDo -0.0 //excitation_energy_of_A[MeV] +0.0 //List_of_excitation_energy_of_A[MeV] #===== end of file diff --git a/working/reactionConfig2.txt b/working/reactionConfig2.txt new file mode 100644 index 0000000..801cd8c --- /dev/null +++ b/working/reactionConfig2.txt @@ -0,0 +1,25 @@ +16 //beam_A +8 //beam_Z +2 //target_A +1 //target_Z +1 //recoil_light_A +1 //recoil-light_Z +10.0 //beam-energy_in_MeV/u +0.000 //beam-energy_sigma_in_MeV/u +0.000 //beam-angle_in_mrad +0.000 //beam-emittance_in_mrad +0.00 //x_offset_of_Beam_in_mm +0.00 //y_offset_of_Beam_in_mm +100000 //number_of_Event_being_generated +false //isTargetScattering +0.913 //target_density_in_g/cm3 +2.2e-4 //targetThickness_in_cm +../SRIM/20F_in_CD2.txt //stopping_power_for_beam +../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil +../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil +false //isDacay +32 //decayNucleus_A +14 //decayNucleus_Z +false //isReDo +0.0 //List_of_excitation_energy_of_A[MeV] +#===== end of file