diff --git a/Armory/ClassDetGeo.h b/Armory/ClassDetGeo.h index a50c76f..c401a96 100644 --- a/Armory/ClassDetGeo.h +++ b/Armory/ClassDetGeo.h @@ -30,13 +30,14 @@ struct Array{ double zMin, zMax; void DeduceAbsolutePos(){ + nDet = pos.size(); detPos.clear(); for(int id = 0; id < nDet; id++){ if( firstPos > 0 ) detPos.push_back(firstPos + pos[id]); if( firstPos < 0 ) detPos.push_back(firstPos - pos[nDet - 1 - id]); - ///printf("%d | %f, %f \n", id, pos[id], detPos[id]); + // printf("%d | %f, %f \n", id, pos[id], detPos[id]); } zMin = TMath::Min(detPos.front(), detPos.back()) - (firstPos < 0 ? detLength : 0); @@ -65,8 +66,10 @@ struct Array{ class DetGeo { public: - DetGeo(); - ~DetGeo(); + DetGeo(){}; + DetGeo(TString detGeoTxt){ LoadDetectorGeo(detGeoTxt, false);} + DetGeo(TMacro * macro){ LoadDetectorGeo(macro, false);} + ~DetGeo(){}; double Bfield; /// T int BfieldSign ; /// sign of B-field @@ -100,24 +103,19 @@ private: }; -inline DetGeo::DetGeo(){ - -} - -inline DetGeo::~DetGeo(){ - -} - inline bool DetGeo::LoadDetectorGeo(TString fileName, bool verbose){ TMacro * haha = new TMacro(); if( haha->ReadFile(fileName) > 0 ) { if( LoadDetectorGeo(haha, verbose) ){ + delete haha; return true; }else{ + delete haha; return false; } }else{ + delete haha; return false; } } @@ -142,7 +140,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){ std::vector str = AnalysisLib::SplitStr(macro->GetListOfLines()->At(i)->GetName(), " "); - //printf("%3d | %s\n", i, str[0].c_str()); + // printf("%3d | %s\n", i, str[0].c_str()); if( str[0].find("####") != std::string::npos ) break; if( str[0].find("#===") != std::string::npos ) { diff --git a/Armory/ClassReactionConfig.h b/Armory/ClassReactionConfig.h index ebab454..a90b89a 100644 --- a/Armory/ClassReactionConfig.h +++ b/Armory/ClassReactionConfig.h @@ -17,6 +17,8 @@ class ReactionConfig { public: ReactionConfig(){} + ReactionConfig(TString configTxt){ LoadReactionConfig(configTxt);} + ReactionConfig(TMacro * macro){ LoadReactionConfig(macro);} ~ReactionConfig(){} int beamA, beamZ; @@ -83,11 +85,14 @@ inline bool ReactionConfig::LoadReactionConfig(TString fileName){ TMacro * haha = new TMacro(); if( haha->ReadFile(fileName) > 0 ) { if( LoadReactionConfig(haha) ){ + delete haha; return true; }else{ + delete haha; return false; } }else{ + delete haha; return false; } } diff --git a/Cleopatra/ClassHelios.h b/Cleopatra/ClassHelios.h index 6b19fdf..212c915 100644 --- a/Cleopatra/ClassHelios.h +++ b/Cleopatra/ClassHelios.h @@ -35,26 +35,27 @@ struct trajectory{ int detID, detRowID; int loop; double effLoop; + + void PrintTrajectory(){ + printf("=====================\n"); + printf(" theta : %f deg\n", theta*TMath::RadToDeg()); + printf(" phi : %f deg\n", phi*TMath::RadToDeg()); + printf(" vt : %f mm/ns\n", vt); + printf(" vp : %f mm/ns\n", vp); + printf(" rho : %f mm\n", rho); + printf(" z0 : %f mm\n", z0); + printf(" t0 : %f ns\n", t0); + printf("(x, y, z) : (%f, %f. %f) mm\n", x, y, z); + printf(" R : %f mm\n", R); + printf(" t : %f ns\n", t); + printf(" effLoop : %f cycle\n", effLoop); + printf(" Loop : %d cycle\n", loop); + printf(" detRowID : %d \n", detRowID); + printf(" detID : %d \n", detID); + + } }; -void PrintTrajectory(trajectory a){ - printf("=====================\n"); - printf(" theta : %f deg\n", a.theta*TMath::RadToDeg()); - printf(" phi : %f deg\n", a.phi*TMath::RadToDeg()); - printf(" vt : %f mm/ns\n", a.vt); - printf(" vp : %f mm/ns\n", a.vp); - printf(" rho : %f mm\n", a.rho); - printf(" z0 : %f mm\n", a.z0); - printf(" t0 : %f ns\n", a.t0); - printf("(x, y, z) : (%f, %f. %f) mm\n", a.x, a.y, a.z); - printf(" R : %f mm\n", a.R); - printf(" t : %f ns\n", a.t); - printf(" effLoop : %f cycle\n", a.effLoop); - printf(" Loop : %d cycle\n", a.loop); - printf(" detRowID : %d \n", a.detRowID); - printf(" detID : %d \n", a.detID); - -} class HELIOS{ public: @@ -86,11 +87,11 @@ public: } int DetAcceptance(); - int CalArrayHit(TLorentzVector Pb, int Zb); + int CalArrayHit(TLorentzVector Pb, int Zb, bool debug = false); int CalRecoilHit(TLorentzVector PB, int ZB); //int CalHit(TLorentzVector Pb, int Zb, TLorentzVector PB, int ZB, double xOff = 0, double yOff = 0 ); // return 0 for no hit, 1 for hit - void CalTrajectoryPara(TLorentzVector P, int Z, int id); // id = 0 for Pb, id = 1 for PB. + void CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil); int GetNumberOfDetectorsInSamePos(){return array.mDet;} double GetEnergy(){return e;} @@ -98,25 +99,28 @@ public: /// clockwise rotation for B-field along the z-axis, sign = 1. double XPos(double Zpos, double theta, double phi, double rho, int sign){ + if( TMath::IsNaN(Zpos) ) return TMath::QuietNaN(); return rho * ( TMath::Sin( TMath::Tan(theta) * Zpos / rho - sign * phi ) + sign * TMath::Sin(phi) ) + xOff; } double YPos(double Zpos, double theta, double phi, double rho, int sign){ + if( TMath::IsNaN(Zpos) ) return TMath::QuietNaN(); return rho * sign * (TMath::Cos( TMath::Tan(theta) * Zpos / rho - sign * phi ) - TMath::Cos(phi)) + yOff; } double RPos(double Zpos, double theta, double phi, double rho, int sign){ - double x = XPos(Zpos, theta, phi, rho, sign) ; - double y = YPos(Zpos, theta, phi, rho, sign) ; - return sqrt(x*x+y*y); + if( TMath::IsNaN(Zpos) ) return TMath::QuietNaN(); + double x = XPos(Zpos, theta, phi, rho, sign) ; + double y = YPos(Zpos, theta, phi, rho, sign) ; + return sqrt(x*x+y*y); } - 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 GetXPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } + double GetYPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } + double GetR(double ZPos) { return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : 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, 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 GetRecoilXPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : XPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } + double GetRecoilYPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } + double GetRecoilR(double ZPos) { return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } double GetBField() {return detGeo.Bfield;} double GetDetRadius() {return array.detPerpDist;} @@ -125,6 +129,9 @@ public: trajectory GetTrajectory_B() {return orbitB;} DetGeo GetDetectorGeometry() {return detGeo;} + + TString GetHitMessage() const {return hitMessage;} + TString GetAcceptanceMessage() const {return accMessage;} private: @@ -156,10 +163,12 @@ private: double eB; ///energy of heavy recoil bool isDetReady; + + TString hitMessage; + TString accMessage; //acceptance check double xOff, yOff; // beam position - bool overrideDetDistance; bool overrideFirstPos; bool isCoincidentWithRecoil; @@ -182,6 +191,9 @@ HELIOS::HELIOS(){ isDetReady = false; + hitMessage = ""; + accMessage = ""; + overrideDetDistance = false; overrideFirstPos = false; isCoincidentWithRecoil = false; @@ -220,40 +232,43 @@ int HELIOS::DetAcceptance(){ if( isDetReady == false ) return 0; // -1 ========= when recoil direction is not same side of array - if( array.firstPos < 0 && orbitb.z > 0 ) return -1; - if( array.firstPos > 0 && orbitb.z < 0 ) return -1; + if( array.firstPos < 0 && orbitb.z > 0 ) {accMessage = "array at upstream, z is downstream."; return -1;} + if( array.firstPos > 0 && orbitb.z < 0 ) {accMessage = "array at downstream, z is upstream."; return -1;} // -11 ======== rho is too small - if( 2 * orbitb.rho < array.detPerpDist ) return -11; + if( 2 * orbitb.rho < array.detPerpDist ) { accMessage = "rho is too small"; return -11;} // -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist) - if( orbitb.detRowID == -1 ) return -15; + if( orbitb.detRowID == -1 ) {accMessage = "det Row ID == -1"; return -15;} // -10 =========== when rho is too big . - if( detGeo.bore < 2 * orbitb.rho) return -10; + if( detGeo.bore < 2 * orbitb.rho) { accMessage = "rho is too big"; return -10;} // -14 ========== check particle-B hit radius on recoil dectector - if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) return -14; + if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) { + accMessage = "heavy recoil does not hit recoil detector"; + return -14; + } //if( isCoincidentWithRecoil && (orbitB.R > rhoRecoilout || orbitB.R < rhoRecoilin) ) return -14; // -12 ========= check is particle-b was blocked by recoil detector 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; + if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { accMessage = "light recoil blocked by recoil detector"; return -12;} + if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { accMessage = "light recoil blocked by recoil detector"; return -12;} // -13 ========= not more than 3 loops - if( orbitb.loop > 3 ) return -13; + if( orbitb.loop > 3 ) {accMessage = "more than 3 loops."; return -13;} // -2 ========= calculate the "y"-distance from detector center - if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) return -2; + if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) { accMessage = "hit at the XY gap."; return -2;} // -3 ==== when zPos further the range of whole array, more loop would not save - 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; + if( array.firstPos < 0 && orbitb.z < array.detPos[0] - array.detLength ) { accMessage = "hit more upstream than the array length"; return -3; } + if( array.firstPos > 0 && orbitb.z > array.detPos[array.nDet-1] + array.detLength ) { accMessage = "hit more downstream than the array length"; return -3;} // -4 ======== Hit on blacker - 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; + if( array.blocker != 0 && array.firstPos > 0 && array.detPos[0] - array.blocker < orbitb.z && orbitb.z < array.detPos[0] ) { accMessage = "hit blocker"; 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 ) { accMessage = "hit blocker"; return -4;} // 2 ====== when zPos less then the nearest position, more loop may hit int increaseLoopFlag = 0; @@ -264,6 +279,7 @@ int HELIOS::DetAcceptance(){ orbitb.effLoop += 1.0; orbitb.loop += 1; orbitb.t = orbitb.t0 * orbitb.effLoop; + accMessage = " hit less than the nearest array. increase loop "; return 2; } @@ -273,6 +289,7 @@ int HELIOS::DetAcceptance(){ if( array.detPos[i] - array.detLength <= orbitb.z && orbitb.z <= array.detPos[i]) { orbitb.detID = i; detX = ( orbitb.z - (array.detPos[i] + array.detLength/2 ))/ array.detLength * 2 ;// range from -1 , 1 + accMessage = "hit array"; return 1; } } @@ -282,6 +299,7 @@ int HELIOS::DetAcceptance(){ ///printf(" %d | %f < z = %f < %f \n", i, array.detPos[i], orbitb.z, array.detPos[i]+length); orbitb.detID = i; detX = ( orbitb.z - (array.detPos[i] - array.detLength/2 ))/ array.detLength*2 ;// range from -1 , 1 + accMessage = "hit array"; return 1; } } @@ -291,11 +309,11 @@ int HELIOS::DetAcceptance(){ // -5 ======== check hit array gap 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; + if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) { accMessage = "hit array Z-gap"; return -5; }//increaseLoopFlag = 3; } }else{ 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( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) { accMessage = "hit array Z-gap"; return -5; }//increaseLoopFlag = 3; } } if (increaseLoopFlag == 3 ) { @@ -303,16 +321,17 @@ int HELIOS::DetAcceptance(){ orbitb.effLoop += 1.0; orbitb.loop += 1; orbitb.t = orbitb.t0 * orbitb.effLoop; + accMessage = " try one more loop. "; return 3; } - + accMessage = " unknown reason "; return -20; // for unknown reason } -void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){ +void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil){ - if( id == 0 ){ + if( isLightRecoil ){ orbitb.theta = P.Theta(); orbitb.phi = P.Phi(); orbitb.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm @@ -323,10 +342,8 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){ orbitb.detID = -1; orbitb.detRowID = -1; - - } - - if( id == 1 ){ + + }else{ orbitB.theta = P.Theta(); orbitB.phi = P.Phi(); orbitB.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm @@ -337,29 +354,27 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){ orbitB.detID = -1; orbitB.detRowID = -1; - } } -int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ +int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb, bool debug){ e = Pb.E() - Pb.M(); detX = TMath::QuietNaN(); rhoHit = TMath::QuietNaN(); - CalTrajectoryPara(Pb, Zb, 0); + CalTrajectoryPara(Pb, Zb, true); int targetLoop = 1; int inOut = array.detFaceOut == true ? 1: 0; //1 = from Outside, 0 = from inside - bool debug = false; - if( debug ) { printf("===================================\n"); printf("theta : %f deg, phi : %f deg \n", orbitb.theta * TMath::RadToDeg(), orbitb.phi * TMath::RadToDeg()); printf("z0: %f mm, rho : %f mm \n", orbitb.z0, orbitb.rho); printf(" inOut : %d = %s \n", inOut, inOut == 1 ? "Out" : "in"); printf(" z range : %.2f - %.2f \n", detGeo.zMin, detGeo.zMax); + printf(" B-field sign : %d\n", detGeo.BfieldSign); printf("-----------------------------------\n"); } @@ -373,7 +388,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ double phiD = TMath::TwoPi()/array.mDet * i ; double dphi = orbitb.phi - phiD; double aEff = array.detPerpDist - (xOff * TMath::Cos(phiD) + yOff * TMath::Sin(phiD)); - double hahaha = asin( aEff/ orbitb.rho - detGeo.BfieldTheta * sin(dphi)); + double hahaha = asin( aEff/ orbitb.rho - detGeo.BfieldSign * sin(dphi)); int n = 2*targetLoop + inOut; @@ -390,22 +405,12 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ } ///Selection - if( !TMath::IsNaN(zP) && 0< zP/orbitb.z0 && TMath::Max(0, targetLoop-1) < zP/orbitb.z0 && zP/orbitb.z0 < targetLoop ) { + if( !TMath::IsNaN(zP) && 0 < zP/orbitb.z0 && TMath::Max(0, targetLoop-1) < zP/orbitb.z0 && zP/orbitb.z0 < targetLoop ) { zPossible.push_back(zP); dID.push_back(i); } } - /* - if( zPossible.size() == 0 ){ // will not happen - zHit = TMath::QuietNaN(); - xPos = TMath::QuietNaN(); - yPos = TMath::QuietNaN(); - loop = -1; - detID = -1; - detRowID = -1; - return -1 ; - }*/ - + if( debug ) printf("-----------------------------------\n"); double dMin = 1; for( int i = 0; i < (int) zPossible.size(); i++){ @@ -434,6 +439,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ orbitb.z = TMath::QuietNaN(); orbitb.loop = -1; orbitb.detRowID = -1; + hitMessage = "wrong direction."; return - 2; } @@ -447,6 +453,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ orbitb.z = TMath::QuietNaN(); orbitb.loop = -1; orbitb.detRowID = -1; + hitMessage = "not on the detector plan."; return -3; } } @@ -459,12 +466,13 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ orbitb.y = GetYPos(orbitb.z) ; orbitb.R = GetR(orbitb.z); + hitMessage = "successful hit."; return 1; // return 1 when OK } int HELIOS::CalRecoilHit(TLorentzVector PB, int ZB){ - CalTrajectoryPara(PB, ZB, 1); + CalTrajectoryPara(PB, ZB, false); orbitB.z = detGeo.recoilPos; orbitB.x = GetRecoilXPos(detGeo.recoilPos) ; diff --git a/Cleopatra/ClassTransfer.h b/Cleopatra/ClassTransfer.h index 8bf4de8..b2cb4da 100644 --- a/Cleopatra/ClassTransfer.h +++ b/Cleopatra/ClassTransfer.h @@ -21,7 +21,20 @@ //======================================================= class TransferReaction { public: - TransferReaction(); + TransferReaction(){Inititization();}; + TransferReaction(int beamA, int beamZ, + int targetA, int targetZ, + int recoilA, int recoilZ, float beamEnergy_AMeV){ + Inititization(); + SetReactionSimple(beamA, beamZ, + targetA, targetZ, + recoilA, recoilZ, beamEnergy_AMeV); + } + TransferReaction(string reactionConfigFile){ + Inititization(); + SetReactionFromFile(reactionConfigFile); + } + ~TransferReaction(); void SetA(int A, int Z, double Ex); @@ -41,15 +54,13 @@ public: TString GetReactionName(); TString GetReactionName_Latex(); - ReactionConfig GetRectionConfig() { return reaction;} + ReactionConfig GetRectionConfig() { return reactionConfig;} double GetMass_A() const {return mA + ExA;} double GetMass_a() const {return ma;} double GetMass_b() const {return mb;} double GetMass_B() const {return mB + ExB;} - int GetCharge_b() const {return reaction.recoilLightZ;} - double GetCMTotalKE() {return Etot - mA - ma;} double GetQValue() {return mA + ExA + ma - mb - mB - ExB;} double GetMaxExB() {return Etot - mb - mB;} @@ -59,6 +70,8 @@ public: TLorentzVector GetPb(){return Pb;} TLorentzVector GetPB(){return PB;} + void PrintFourVectors(); + void CalReactionConstant(); void Event(double thetaCM_rad, double phiCM_rad); @@ -72,7 +85,7 @@ public: private: - ReactionConfig reaction; + ReactionConfig reactionConfig; string nameA, namea, nameb, nameB; double thetaIN, phiIN; @@ -92,10 +105,12 @@ private: TLorentzVector PA, Pa, Pb, PB; TString format(TString name); + + void Inititization(); }; -TransferReaction::TransferReaction(){ +void TransferReaction::Inititization(){ thetaIN = 0.; phiIN = 0.; @@ -104,7 +119,7 @@ TransferReaction::TransferReaction(){ Setb(1,1); SetB(13,6); TA = 6; - T = TA * reaction.beamA; + T = TA * reactionConfig.beamA; ExA = 0; ExB = 0; @@ -126,8 +141,8 @@ TransferReaction::~TransferReaction(){ void TransferReaction::SetA(int A, int Z, double Ex = 0){ Isotope temp (A, Z); mA = temp.Mass; - reaction.beamA = A; - reaction.beamZ = Z; + reactionConfig.beamA = A; + reactionConfig.beamZ = Z; ExA = Ex; nameA = temp.Name; isReady = false; @@ -137,8 +152,8 @@ void TransferReaction::SetA(int A, int Z, double Ex = 0){ void TransferReaction::Seta(int A, int Z){ Isotope temp (A, Z); ma = temp.Mass; - reaction.targetA = A; - reaction.targetZ = Z; + reactionConfig.targetA = A; + reactionConfig.targetZ = Z; namea = temp.Name; isReady = false; isBSet = false; @@ -146,8 +161,8 @@ void TransferReaction::Seta(int A, int Z){ void TransferReaction::Setb(int A, int Z){ Isotope temp (A, Z); mb = temp.Mass; - reaction.recoilLightA = A; - reaction.recoilLightZ = Z; + reactionConfig.recoilLightA = A; + reactionConfig.recoilLightZ = Z; nameb = temp.Name; isReady = false; isBSet = false; @@ -155,8 +170,8 @@ void TransferReaction::Setb(int A, int Z){ void TransferReaction::SetB(int A, int Z){ Isotope temp (A, Z); mB = temp.Mass; - reaction.recoilHeavyA = A; - reaction.recoilHeavyZ = Z; + reactionConfig.recoilHeavyA = A; + reactionConfig.recoilHeavyZ = Z; nameB = temp.Name; isReady = false; isBSet = true; @@ -164,7 +179,7 @@ void TransferReaction::SetB(int A, int Z){ void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){ this->TA = KEA; - this->T = TA * reaction.beamA; + this->T = TA * reactionConfig.beamA; this->thetaIN = theta; this->phiIN = phi; isReady = false; @@ -174,15 +189,15 @@ void TransferReaction::SetReactionSimple(int beamA, int beamZ, int targetA, int targetZ, int recoilA, int recoilZ, float beamEnergy_AMeV){ - reaction.SetReactionSimple(beamA, beamZ, + reactionConfig.SetReactionSimple(beamA, beamZ, targetA, targetZ, recoilA, recoilZ, beamEnergy_AMeV); - SetA(reaction.beamA, reaction.beamZ); - Seta(reaction.targetA, reaction.targetZ); - Setb(reaction.recoilLightA, reaction.recoilLightZ); - SetB(reaction.recoilHeavyA, reaction.recoilHeavyZ); - SetIncidentEnergyAngle(reaction.beamEnergy, 0, 0); + SetA(reactionConfig.beamA, reactionConfig.beamZ); + Seta(reactionConfig.targetA, reactionConfig.targetZ); + Setb(reactionConfig.recoilLightA, reactionConfig.recoilLightZ); + SetB(reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ); + SetIncidentEnergyAngle(reactionConfig.beamEnergy, 0, 0); CalReactionConstant(); @@ -200,13 +215,13 @@ void TransferReaction::SetExB(double Ex){ void TransferReaction::SetReactionFromFile(string reactionConfigFile){ - if( reaction.LoadReactionConfig(reactionConfigFile) ){ + if( reactionConfig.LoadReactionConfig(reactionConfigFile) ){ - SetA(reaction.beamA, reaction.beamZ); - Seta(reaction.targetA, reaction.targetZ); - Setb(reaction.recoilLightA, reaction.recoilLightZ); - SetB(reaction.recoilHeavyA, reaction.recoilHeavyZ); - SetIncidentEnergyAngle(reaction.beamEnergy, 0, 0); + SetA(reactionConfig.beamA, reactionConfig.beamZ); + Seta(reactionConfig.targetA, reactionConfig.targetZ); + Setb(reactionConfig.recoilLightA, reactionConfig.recoilLightZ); + SetB(reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ); + SetIncidentEnergyAngle(reactionConfig.beamEnergy, 0, 0); CalReactionConstant(); @@ -245,9 +260,9 @@ TString TransferReaction::GetReactionName_Latex(){ void TransferReaction::CalReactionConstant(){ if( !isBSet){ - reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA; - reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ; - Isotope temp (reaction.recoilHeavyA, reaction.recoilHeavyZ); + reactionConfig.recoilHeavyA = reactionConfig.beamA + reactionConfig.targetA - reactionConfig.recoilLightA; + reactionConfig.recoilHeavyZ = reactionConfig.beamZ + reactionConfig.targetZ - reactionConfig.recoilLightZ; + Isotope temp (reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ); mB = temp.Mass; isBSet = true; } @@ -267,6 +282,21 @@ void TransferReaction::CalReactionConstant(){ isReady = true; } +void TransferReaction::PrintFourVectors(){ + + printf("A : %10.2f %10.2f %10.2f %10.2f\n", PA.E(), PA.Px(), PA.Py(), PA.Pz()); + printf("a : %10.2f %10.2f %10.2f %10.2f\n", Pa.E(), Pa.Px(), Pa.Py(), Pa.Pz()); + printf("b : %10.2f %10.2f %10.2f %10.2f\n", Pb.E(), Pb.Px(), Pb.Py(), Pb.Pz()); + printf("B : %10.2f %10.2f %10.2f %10.2f\n", PB.E(), PB.Px(), PB.Py(), PB.Pz()); + printf("-------------------------------------------------------\n"); + printf("B : %10.2f %10.2f %10.2f %10.2f\n", + PA.E() + Pa.E() - Pb.E() - PB.E(), + PA.Px() + Pa.Px() - Pb.Px() - PB.Px(), + PA.Py() + Pa.Py() - Pb.Py() - PB.Py(), + PA.Pz() + Pa.Pz() - Pb.Pz() - PB.Pz()); + +} + void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){ if( isReady == false ){ @@ -325,7 +355,7 @@ std::pair TransferReaction::CalExThetaCM(double e, double z, dou double mass = mb; double massB = mB; double y = e + mass; - double slope = 299.792458 * reaction.recoilLightZ * abs(Bfield) / TMath::TwoPi() * beta / 1000.; // MeV/mm; + double slope = 299.792458 * reactionConfig.recoilLightZ * abs(Bfield) / TMath::TwoPi() * beta / 1000.; // MeV/mm; double alpha = slope/beta; double G = alpha * gamma * beta * perpDist ; double Z = alpha * gamma * beta * z; diff --git a/Cleopatra/Transfer.h b/Cleopatra/Transfer.h index 1a41fcd..861b15b 100644 --- a/Cleopatra/Transfer.h +++ b/Cleopatra/Transfer.h @@ -32,16 +32,16 @@ void Transfer( TString filename = "reaction.dat"){ /// when no file, no output. //############################################# Set Reaction - TransferReaction reaction; - reaction.SetReactionFromFile(basicConfig); + TransferReaction transfer; + transfer.SetReactionFromFile(basicConfig); printf("*****************************************************************\n"); - printf("*\e[1m\e[33m %27s \e[0m*\n", reaction.GetReactionName().Data()); + printf("*\e[1m\e[33m %27s \e[0m*\n", transfer.GetReactionName().Data()); printf("*****************************************************************\n"); printf("----- loading reaction setting from %s. \n", basicConfig.c_str()); printf("\e[32m#################################### Beam \e[0m\n"); - const ReactionConfig reactionConfig = reaction.GetRectionConfig(); + const ReactionConfig reactionConfig = transfer.GetRectionConfig(); reactionConfig.PrintReactionConfig(); @@ -56,12 +56,12 @@ void Transfer( const DetGeo detGeo = helios.GetDetectorGeometry(); printf("==================================== E-Z plot slope\n"); - double betaRect = reaction.GetReactionBeta() ; - double gamma = reaction.GetReactionGamma(); - double mb = reaction.GetMass_b(); - double pCM = reaction.GetMomentumbCM(); + double betaRect = transfer.GetReactionBeta() ; + double gamma = transfer.GetReactionGamma(); + double mb = transfer.GetMass_b(); + double pCM = transfer.GetMomentumbCM(); double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass - double slope = 299.792458 * reaction.GetCharge_b() * abs(helios.GetBField()) / TMath::TwoPi() * betaRect / 1000.; /// MeV/mm + double slope = 299.792458 * reactionConfig.recoilLightZ * abs(helios.GetBField()) / TMath::TwoPi() * betaRect / 1000.; /// MeV/mm printf(" e-z slope : %f MeV/mm\n", slope); double intercept = q/gamma - mb; // MeV printf(" e-z intercept (ground state) : %f MeV\n", intercept); @@ -72,11 +72,11 @@ void Transfer( keyParaOut = fopen (filename.Data(), "w+"); printf("=========== save key reaction constants to %s \n", filename.Data()); - fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetMass_b(), "mass_b"); - fprintf(keyParaOut, "%-15d //%s\n", reaction.GetCharge_b(), "charge_b"); - fprintf(keyParaOut, "%-15.8f //%s\n", reaction.GetReactionBeta(), "betaCM"); - fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetCMTotalEnergy(), "Ecm"); - fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetMass_B(), "mass_B"); + fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_b(), "mass_b"); + fprintf(keyParaOut, "%-15d //%s\n", reactionConfig.recoilLightZ, "charge_b"); + fprintf(keyParaOut, "%-15.8f //%s\n", transfer.GetReactionBeta(), "betaCM"); + fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetCMTotalEnergy(), "Ecm"); + fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_B(), "mass_B"); fprintf(keyParaOut, "%-15.4f //%s\n", slope/betaRect, "alpha=slope/betaRect"); fflush(keyParaOut); @@ -147,9 +147,9 @@ void Transfer( printf("%3s | %7s | %5s | %3s | %10s | %5s \n", "", "Ex[MeV]", "Xsec", "SF", "sigma[MeV]", "y0[MeV]"); printf("----+---------+------+-----+------------+--------\n"); for(int i = 0; i < n ; i++){ - reaction.SetExB(ExKnown[i]); - reaction.CalReactionConstant(); - kCM.push_back(reaction.GetMomentumbCM()); + transfer.SetExB(ExKnown[i]); + transfer.CalReactionConstant(); + kCM.push_back(transfer.GetMomentumbCM()); y0.push_back(TMath::Sqrt(mb*mb + kCM[i]*kCM[i])/gamma - mb); if( reactionConfig.isDecay ) { TLorentzVector temp(0,0,0,0); @@ -169,9 +169,9 @@ void Transfer( ExKnown.push_back(0.0); ExStrength.push_back(1.0); ExWidth.push_back(0.0); - reaction.SetExB(ExKnown[0]); - reaction.CalReactionConstant(); - kCM.push_back(reaction.GetMomentumbCM()); + transfer.SetExB(ExKnown[0]); + transfer.CalReactionConstant(); + kCM.push_back(transfer.GetMomentumbCM()); y0.push_back(TMath::Sqrt(mb*mb + kCM[0]*kCM[0])/gamma - mb); } @@ -215,7 +215,7 @@ void Transfer( TMacro reactionData(filename.Data()); double KEAmean = reactionConfig.beamEnergy; TString str; - str.Form("%s @ %.2f MeV/u", reaction.GetReactionName_Latex().Data(), KEAmean); + str.Form("%s @ %.2f MeV/u", transfer.GetReactionName_Latex().Data(), KEAmean); config.SetName(str.Data()); config.Write("reactionConfig"); detGeoTxt.Write("detGeo"); @@ -486,7 +486,7 @@ void Transfer( //==== Set Ex of A ExAID = gRandom->Integer(nExA); ExA = ExAList[ExAID]; - reaction.SetExA(ExA); + transfer.SetExA(ExA); //==== Set Ex of B if( ExKnown.size() == 1 ) { @@ -496,7 +496,7 @@ void Transfer( ExID = exDist->GetRandom(); Ex = ExKnown[ExID]+ (ExWidth[ExID] == 0 ? 0 : gRandom->Gaus(0, ExWidth[ExID])); } - reaction.SetExB(Ex); + transfer.SetExB(Ex); //==== Set incident beam KEA = reactionConfig.beamEnergy; @@ -514,9 +514,9 @@ void Transfer( phi = 0.0; //==== for taregt scattering - reaction.SetIncidentEnergyAngle(KEA, theta, 0.); - reaction.CalReactionConstant(); - TLorentzVector PA = reaction.GetPA(); + transfer.SetIncidentEnergyAngle(KEA, theta, 0.); + transfer.CalReactionConstant(); + TLorentzVector PA = transfer.GetPA(); //depth = 0; if( isTargetScattering ){ @@ -525,9 +525,9 @@ void Transfer( msA.SetTarget(density, depth); TLorentzVector PAnew = msA.Scattering(PA); KEAnew = msA.GetKE()/reactionConfig.beamA; - reaction.SetIncidentEnergyAngle(KEAnew, theta, phi); - reaction.CalReactionConstant(); - Ecm = reaction.GetCMTotalKE(); + transfer.SetIncidentEnergyAngle(KEAnew, theta, phi); + transfer.CalReactionConstant(); + Ecm = transfer.GetCMTotalKE(); } //==== Calculate thetaCM, phiCM @@ -541,9 +541,9 @@ void Transfer( double phiCM = TMath::TwoPi() * gRandom->Rndm(); //==== Calculate reaction - reaction.Event(thetaCM, phiCM); - TLorentzVector Pb = reaction.GetPb(); - TLorentzVector PB = reaction.GetPB(); + transfer.Event(thetaCM, phiCM); + TLorentzVector Pb = transfer.GetPb(); + TLorentzVector PB = transfer.GetPB(); //==== Calculate energy loss of scattered and recoil in target if( isTargetScattering ){ @@ -593,7 +593,7 @@ void Transfer( ///printf(" thetaCM : %f \n", thetaCM * TMath::RadToDeg()); if( Tb > 0 || TB > 0 ){ - helios.CalArrayHit(Pb, reaction.GetCharge_b()); + helios.CalArrayHit(Pb, reactionConfig.recoilLightZ); helios.CalRecoilHit(PB, new_zB); hit = 2; while( hit > 1 ){ hit = helios.DetAcceptance(); } /// while hit > 1, goto next loop; @@ -640,17 +640,17 @@ void Transfer( //other recoil detectors if ( detGeo.recoilPos1 != 0 ){ - xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1); - yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1); + xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1); + yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1); rhoRecoil1 = helios.GetRecoilR(detGeo.recoilPos1); } if ( detGeo.recoilPos2 != 0 ){ - xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2); - yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2); + xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2); + yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2); rhoRecoil2 = helios.GetRecoilR(detGeo.recoilPos2); } - std::pair ExThetaCM = reaction.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius()); + std::pair ExThetaCM = transfer.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius()); ExCal = ExThetaCM.first; thetaCMCal = ExThetaCM.second; diff --git a/README.md b/README.md index 3ad59b9..6ee3403 100644 --- a/README.md +++ b/README.md @@ -12,9 +12,9 @@ Analysis ├── SOLARIS.sh // bash script to define some env variable and functions -├── armory // analysis codes, independent from experiment. +├── Armory // analysis codes, independent from experiment. -├── Cleopatra // Swaper for DWBA code Ptolomey +├── Cleopatra // Swaper for DWBA code Ptolomey and simulation ├── data_raw // should be the symbolic link to the raw data, created by SetUpNewExp