diff --git a/Armory/ClassDetGeo.h b/Armory/ClassDetGeo.h index e79b3db..0a239b4 100644 --- a/Armory/ClassDetGeo.h +++ b/Armory/ClassDetGeo.h @@ -15,6 +15,8 @@ struct Array{ + bool enable; + double detPerpDist; /// distance from axis double detWidth; /// width double detLength; /// length @@ -85,11 +87,7 @@ public: double elumPos1, elumPos2; /// imaginary elum, only sensitive to light recoil //===================1st array - Array array1; - - //==================2nd array - bool use2ndArray; - Array array2; + Array array[2]; double zMin, zMax; /// range of detectors bool isCoincidentWithRecoil; @@ -129,9 +127,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){ TList * haha = macro->GetListOfLines(); int numLine = (haha)->GetSize(); - array1.pos.clear(); - array2.pos.clear(); - use2ndArray = false; + for( int i = 0; i < 2 ; i++) array[i].pos.clear(); int detFlag = 0; int detLine = 0; @@ -170,45 +166,35 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){ if ( detLine == 10 ) elumPos2 = atof(str[0].c_str()); } - if( detFlag == 1){ - if ( detLine == 0 ) array1.detPerpDist = atof(str[0].c_str()); - if ( detLine == 1 ) array1.detWidth = atof(str[0].c_str()); - if ( detLine == 2 ) array1.detLength = atof(str[0].c_str()); - if ( detLine == 3 ) array1.blocker = atof(str[0].c_str()); - if ( detLine == 4 ) array1.firstPos = atof(str[0].c_str()); - if ( detLine == 5 ) array1.eSigma = atof(str[0].c_str()); - if ( detLine == 6 ) array1.zSigma = atof(str[0].c_str()); - if ( detLine == 7 ) array1.detFaceOut = str[0] == "Out" ? true : false; - if ( detLine == 8 ) array1.mDet = atoi(str[0].c_str()); - if ( detLine >= 9 ) array1.pos.push_back(atof(str[0].c_str())); - } - - if( detFlag == 2){ - if ( detLine == 0 ) use2ndArray = str[0] == "true" ? true : false; - if ( detLine == 1 ) array2.detPerpDist = atof(str[0].c_str()); - if ( detLine == 2 ) array2.detWidth = atof(str[0].c_str()); - if ( detLine == 3 ) array2.detLength = atof(str[0].c_str()); - if ( detLine == 4 ) array2.blocker = atof(str[0].c_str()); - if ( detLine == 5 ) array2.firstPos = atof(str[0].c_str()); - if ( detLine == 6 ) array2.eSigma = atof(str[0].c_str()); - if ( detLine == 7 ) array2.zSigma = atof(str[0].c_str()); - if ( detLine == 8 ) array2.detFaceOut = str[0] == "Out" ? true : false; - if ( detLine == 9 ) array2.mDet = atoi(str[0].c_str()); - if ( detLine >= 10 ) array2.pos.push_back(atof(str[0].c_str())); + if( detFlag > 0){ + unsigned short ID = detFlag - 1; + if ( detLine == 0 ) array[ID].enable = str[0] == "true" ? true : false; + if ( detLine == 1 ) array[ID].detPerpDist = atof(str[0].c_str()); + if ( detLine == 2 ) array[ID].detWidth = atof(str[0].c_str()); + if ( detLine == 3 ) array[ID].detLength = atof(str[0].c_str()); + if ( detLine == 4 ) array[ID].blocker = atof(str[0].c_str()); + if ( detLine == 5 ) array[ID].firstPos = atof(str[0].c_str()); + if ( detLine == 6 ) array[ID].eSigma = atof(str[0].c_str()); + if ( detLine == 7 ) array[ID].zSigma = atof(str[0].c_str()); + if ( detLine == 8 ) array[ID].detFaceOut = str[0] == "Out" ? true : false; + if ( detLine == 9 ) array[ID].mDet = atoi(str[0].c_str()); + if ( detLine >= 10 ) array[ID].pos.push_back(atof(str[0].c_str())); } detLine ++; } - array1.DeduceAbsolutePos(); - array2.DeduceAbsolutePos(); + zMin = 99999; + zMax = -99999; - zMin = array1.zMin; - zMax = array1.zMax; - - if( use2ndArray) { - zMax = TMath::Min(array1.zMax, array2.zMax); - zMin = TMath::Min(array1.zMin, array2.zMin); + for( int i = 0; i < 2; i ++ ){ + array[i].DeduceAbsolutePos(); + if (array[i].enable ) { + double zmax = TMath::Max(array[i].zMin, array[i].zMax); + double zmin = TMath::Min(array[i].zMin, array[i].zMax); + if( zmax > zMax ) zMax = zmax; + if( zmin < zMin ) zMin = zmin; + } } if( verbose ) Print(false); @@ -226,12 +212,13 @@ inline void DetGeo::Print(bool printAll) const{ } printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", recoilPos, recoilInnerRadius, recoilOuterRadius); - printf("------------------------------------- Detector Position \n"); - array1.PrintArray(); + for( int i = 0; i < 2 ; i++){ + + if( printAll || array[i].enable ) { + printf("-----------------------------------%d-th Detector Position \n", i); + array[i].PrintArray(); + } - if( printAll || use2ndArray){ - printf("--------------------------------- 2nd Detector Position \n"); - array2.PrintArray(); } if( elumPos1 != 0 || elumPos2 != 0 || recoilPos1 != 0 || recoilPos2 != 0){ diff --git a/Armory/ClassReactionConfig.h b/Armory/ClassReactionConfig.h index a6a552a..145b204 100644 --- a/Armory/ClassReactionConfig.h +++ b/Armory/ClassReactionConfig.h @@ -63,14 +63,14 @@ public: float targetThickness; ///targetThickness_in_cm std::string beamStoppingPowerFile; ///stopping_power_for_beam - Recoil recoil1, recoil2; + Recoil recoil[2]; int numEvents; ///number_of_Event_being_generated bool isRedo; ///isReDo void SetReactionSimple(int beamA, int beamZ, int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV); + int recoilA, int recoilZ, float beamEnergy_AMeV, unsigned short ID); bool LoadReactionConfig(TString fileName); bool LoadReactionConfig(TMacro * macro); @@ -83,17 +83,17 @@ private: inline void ReactionConfig::SetReactionSimple(int beamA, int beamZ, int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV){ + int recoilA, int recoilZ, float beamEnergy_AMeV, unsigned short ID){ this->beamA = beamA; this->beamZ = beamZ; this->targetA = targetA; this->targetZ = targetZ; - this->recoil1.lightA = recoilA; - this->recoil1.lightZ = recoilZ; - recoil1.heavyA = this->beamA + this->targetA - recoil1.lightA; - recoil1.heavyZ = this->beamZ + this->targetZ - recoil1.lightZ; + this->recoil[ID].lightA = recoilA; + this->recoil[ID].lightZ = recoilZ; + recoil[ID].heavyA = this->beamA + this->targetA - recoil[ID].lightA; + recoil[ID].heavyZ = this->beamZ + this->targetZ - recoil[ID].lightZ; beamEnergy = beamEnergy_AMeV; beamEnergySigma = 0; @@ -168,27 +168,16 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){ if( recoilLine == 16 ) isRedo = str[0].compare("true" ) == 0 ? true : false; } - if( recoilFlag == 1 ){ + if( recoilFlag > 0 ){ - if( recoilLine == 0 ) recoil1.lightA = atoi(str[0].c_str()); - if( recoilLine == 1 ) recoil1.lightZ = atoi(str[0].c_str()); - if( recoilLine == 2 ) recoil1.lightStoppingPowerFile = str[0]; - if( recoilLine == 3 ) recoil1.heavyStoppingPowerFile = str[0]; - if( recoilLine == 4 ) recoil1.isDecay = str[0].compare("true") == 0 ? true : false; - if( recoilLine == 5 ) recoil1.decayA = atoi(str[0].c_str()); - if( recoilLine == 6 ) recoil1.decayZ = atoi(str[0].c_str()); - - } - - if( recoilFlag == 2 ){ - - if( recoilLine == 0 ) recoil2.lightA = atoi(str[0].c_str()); - if( recoilLine == 1 ) recoil2.lightZ = atoi(str[0].c_str()); - if( recoilLine == 2 ) recoil2.lightStoppingPowerFile = str[0]; - if( recoilLine == 3 ) recoil2.heavyStoppingPowerFile = str[0]; - if( recoilLine == 4 ) recoil2.isDecay = str[0].compare("true") == 0 ? true : false; - if( recoilLine == 5 ) recoil2.decayA = atoi(str[0].c_str()); - if( recoilLine == 6 ) recoil2.decayZ = atoi(str[0].c_str()); + unsigned ID = recoilFlag - 1; + if( recoilLine == 0 ) recoil[ID].lightA = atoi(str[0].c_str()); + if( recoilLine == 1 ) recoil[ID].lightZ = atoi(str[0].c_str()); + if( recoilLine == 2 ) recoil[ID].lightStoppingPowerFile = str[0]; + if( recoilLine == 3 ) recoil[ID].heavyStoppingPowerFile = str[0]; + if( recoilLine == 4 ) recoil[ID].isDecay = str[0].compare("true") == 0 ? true : false; + if( recoilLine == 5 ) recoil[ID].decayA = atoi(str[0].c_str()); + if( recoilLine == 6 ) recoil[ID].decayZ = atoi(str[0].c_str()); } @@ -196,12 +185,10 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){ } - recoil1.heavyA = beamA + targetA - recoil1.lightA; - recoil1.heavyZ = beamZ + targetZ - recoil1.lightZ; - - recoil2.heavyA = beamA + targetA - recoil2.lightA; - recoil2.heavyZ = beamZ + targetZ - recoil2.lightZ; - + for( int i = 0; i < 2; i++){ + recoil[i].heavyA = beamA + targetA - recoil[i].lightA; + recoil[i].heavyZ = beamZ + targetZ - recoil[i].lightZ; + } return true; } @@ -227,9 +214,9 @@ inline void ReactionConfig::Print() const{ printf(" beam stopping file : %s \n", beamStoppingPowerFile.c_str()); } - printf("------------------------------ Recoil-1\n"); recoil1.Print(); - - printf("------------------------------ Recoil-2\n"); recoil2.Print(); + for( int i = 0; i < 2; i ++ ){ + printf("------------------------------ Recoil-%d\n", i); recoil[i].Print(); + } printf("=====================================================\n"); diff --git a/Cleopatra/ClassHelios.h b/Cleopatra/ClassHelios.h index 212c915..30ea62b 100644 --- a/Cleopatra/ClassHelios.h +++ b/Cleopatra/ClassHelios.h @@ -25,133 +25,120 @@ //======================================================= struct trajectory{ - double theta, phi; - double vt, vp; // tranvser and perpendicular velocity - double rho; // orbit radius - double z0, t0; // position cycle - double x, y, z; // hit position - double t; //actual orbit time; - double R; //hit radius = sqrt(x^2+y^2); - int detID, detRowID; - int loop; - double effLoop; + double theta, phi; + double vt, vp; // tranvser and perpendicular velocity + double rho; // orbit radius + double z0, t0; // position cycle + double x, y, z; // hit position + double t; //actual orbit time; + double R; //hit radius = sqrt(x^2+y^2); + 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(){ + 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 Clear(){ + theta = TMath::QuietNaN(); + phi = TMath::QuietNaN(); + vt = TMath::QuietNaN(); + vp = TMath::QuietNaN(); + rho = TMath::QuietNaN(); + z0 = TMath::QuietNaN(); + t0 = TMath::QuietNaN(); + x = TMath::QuietNaN(); + y = TMath::QuietNaN(); + z = TMath::QuietNaN(); + effLoop = TMath::QuietNaN(); + detID = -1; + detRowID = -1; + loop = -1; + } }; - class HELIOS{ public: - HELIOS(); - ~HELIOS(); - - void SetCoincidentWithRecoil(bool TorF){ this->isCoincidentWithRecoil = TorF;} - bool GetCoincidentWithRecoil(){return this->isCoincidentWithRecoil;} - bool SetDetectorGeometry(std::string filename); - void SetBeamPosition(double x, double y) { xOff = x; yOff = y;} - - 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->array.firstPos = firstPos; - } - void OverrideDetectorDistance(double perpDist){ - overrideDetDistance = true; - printf("------ Overriding Detector Distance to : %8.2f mm \n", perpDist); - this->array.detPerpDist = perpDist; - } - - void SetDetectorOutside(bool isOutside){ - this->array.detFaceOut = isOutside; - printf(" Detectors are facing %s\n", array.detFaceOut ? "outside": "inside" ); - } - - int DetAcceptance(); - 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, bool isLightRecoil); - - int GetNumberOfDetectorsInSamePos(){return array.mDet;} - double GetEnergy(){return e;} - double GetDetX(){return detX;} // position in each detector, range from -1, 1 - - /// 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){ + HELIOS(); + ~HELIOS(); + + void SetCoincidentWithRecoil(bool TorF){ this->isCoincidentWithRecoil = TorF;} + bool GetCoincidentWithRecoil(){return this->isCoincidentWithRecoil;} + bool SetDetectorGeometry(std::string filename, unsigned short ID); + void SetBeamPosition(double x, double y) { xOff = x; yOff = y;} + + void OverrideMagneticField(double BField); + void OverrideMagneticFieldDirection(double BfieldThetaInDeg); + void OverrideFirstPos(double firstPos); + void OverrideDetectorDistance(double perpDist); + void OverrideDetectorFacing(bool isOutside); + + int CheckDetAcceptance(); + int CalArrayHit(TLorentzVector Pb, bool debug = false); + int CalRecoilHit(TLorentzVector PB); + void CalTrajectoryPara(TLorentzVector P, bool isLightRecoil); + + int GetNumberOfDetectorsInSamePos(){return array.mDet;} + double GetEnergy()const {return e;} + double GetDetX() const {return detX;} // position in each detector, range from -1, 1 + + /// 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){ 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 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 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;} + } + + 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 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); } + + void PrintGeometry() const; - trajectory GetTrajectory_b() {return orbitb;} - trajectory GetTrajectory_B() {return orbitB;} + double GetBField() const {return detGeo.Bfield;} + double GetDetRadius() const {return array.detPerpDist;} - DetGeo GetDetectorGeometry() {return detGeo;} - - TString GetHitMessage() const {return hitMessage;} - TString GetAcceptanceMessage() const {return accMessage;} - + trajectory GetTrajectory_b() const {return orbitb;} + trajectory GetTrajectory_B() const {return orbitB;} + + DetGeo GetDetectorGeometry() const {return detGeo;} + + TString GetHitMessage() const {return hitMessage;} + TString GetAcceptanceMessage() const {return accMessage;} + private: - void ClearTrajectory(trajectory t){ - t.theta = TMath::QuietNaN(); - t.phi = TMath::QuietNaN(); - t.vt = TMath::QuietNaN(); - t.vp = TMath::QuietNaN(); - t.rho = TMath::QuietNaN(); - t.z0 = TMath::QuietNaN(); - t.t0 = TMath::QuietNaN(); - t.x = TMath::QuietNaN(); - t.y = TMath::QuietNaN(); - t.z = TMath::QuietNaN(); - t.effLoop = TMath::QuietNaN(); - t.detID = -1; - t.detRowID = -1; - t.loop = -1; - } - DetGeo detGeo; Array array; @@ -178,25 +165,25 @@ private: HELIOS::HELIOS(){ - ClearTrajectory(orbitb); - ClearTrajectory(orbitB); + orbitb.Clear(); + orbitB.Clear(); - e = TMath::QuietNaN(); - eB = TMath::QuietNaN(); - detX = TMath::QuietNaN(); - rhoHit = TMath::QuietNaN(); - - xOff = 0.0; - yOff = 0.0; - - isDetReady = false; + e = TMath::QuietNaN(); + eB = TMath::QuietNaN(); + detX = TMath::QuietNaN(); + rhoHit = TMath::QuietNaN(); + + xOff = 0.0; + yOff = 0.0; + + isDetReady = false; - hitMessage = ""; - accMessage = ""; + hitMessage = ""; + accMessage = ""; - overrideDetDistance = false; - overrideFirstPos = false; - isCoincidentWithRecoil = false; + overrideDetDistance = false; + overrideFirstPos = false; + isCoincidentWithRecoil = false; } @@ -204,28 +191,73 @@ HELIOS::~HELIOS(){ } -bool HELIOS::SetDetectorGeometry(std::string filename){ - - if( detGeo.LoadDetectorGeo(filename)) { - - if( detGeo.use2ndArray ){ - array = detGeo.array2; - }else{ - array = detGeo.array1; - } - - isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil; - - isDetReady = true; - }else{ - printf("cannot read file %s.\n", filename.c_str()); - isDetReady = false; - } - - return isDetReady; +void HELIOS::OverrideMagneticField(double BField){ + this->detGeo.Bfield = BField; + this->detGeo.BfieldSign = BField > 0 ? 1: -1; } -int HELIOS::DetAcceptance(){ +void HELIOS::OverrideMagneticFieldDirection(double BfieldThetaInDeg){ + this->detGeo.BfieldTheta = BfieldThetaInDeg; +} + +void HELIOS::OverrideFirstPos(double firstPos){ + overrideFirstPos = true; + printf("------ Overriding FirstPosition to : %8.2f mm \n", firstPos); + this->array.firstPos = firstPos; +} + +void HELIOS::OverrideDetectorDistance(double perpDist){ + overrideDetDistance = true; + printf("------ Overriding Detector Distance to : %8.2f mm \n", perpDist); + this->array.detPerpDist = perpDist; +} + +void HELIOS::OverrideDetectorFacing(bool isOutside){ + this->array.detFaceOut = isOutside; + printf(" Detectors are facing %s\n", array.detFaceOut ? "outside": "inside" ); +} + +bool HELIOS::SetDetectorGeometry(std::string filename, unsigned short ID){ + + if( detGeo.LoadDetectorGeo(filename, false)) { + + array = detGeo.array[ID]; + isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil; + isDetReady = true; + + }else{ + printf("cannot read file %s.\n", filename.c_str()); + isDetReady = false; + } + + return isDetReady; +} + +void HELIOS::PrintGeometry() const{ + + printf("=====================================================\n"); + printf(" B-field: %8.2f T, Theta : %6.2f deg \n", detGeo.Bfield, detGeo.BfieldTheta); + if( detGeo.BfieldTheta != 0.0 ) { + 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"); + array.PrintArray(); + + if( detGeo.elumPos1 != 0 || detGeo.elumPos2 != 0 || detGeo.recoilPos1 != 0 || detGeo.recoilPos2 != 0){ + printf("=================================== Auxillary/Imaginary Detectors\n"); + } + if( detGeo.elumPos1 != 0 ) printf(" Elum 1 pos.: %f mm \n", detGeo.elumPos1); + if( detGeo.elumPos2 != 0 ) printf(" Elum 2 pos.: %f mm \n", detGeo.elumPos2); + if( detGeo.recoilPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", detGeo.recoilPos1); + if( detGeo.recoilPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", detGeo.recoilPos2); + printf("=====================================================\n"); + + +} + +int HELIOS::CheckDetAcceptance(){ //CalArrayHit and CalRecoilHit must be done before. @@ -329,12 +361,12 @@ int HELIOS::DetAcceptance(){ return -20; // for unknown reason } -void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil){ +void HELIOS::CalTrajectoryPara(TLorentzVector P, bool isLightRecoil){ if( isLightRecoil ){ orbitb.theta = P.Theta(); orbitb.phi = P.Phi(); - orbitb.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm + orbitb.rho = P.Pt() / abs(detGeo.Bfield) / P.GetUniqueID() / 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 @@ -346,7 +378,7 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil){ }else{ orbitB.theta = P.Theta(); orbitB.phi = P.Phi(); - orbitB.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm + orbitB.rho = P.Pt() / abs(detGeo.Bfield) / P.GetUniqueID() / 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 @@ -357,13 +389,13 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil){ } } -int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb, bool debug){ +int HELIOS::CalArrayHit(TLorentzVector Pb, bool debug){ e = Pb.E() - Pb.M(); detX = TMath::QuietNaN(); rhoHit = TMath::QuietNaN(); - CalTrajectoryPara(Pb, Zb, true); + CalTrajectoryPara(Pb, true); int targetLoop = 1; int inOut = array.detFaceOut == true ? 1: 0; //1 = from Outside, 0 = from inside @@ -470,9 +502,9 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb, bool debug){ return 1; // return 1 when OK } -int HELIOS::CalRecoilHit(TLorentzVector PB, int ZB){ +int HELIOS::CalRecoilHit(TLorentzVector PB){ - CalTrajectoryPara(PB, ZB, false); + CalTrajectoryPara(PB, false); orbitB.z = detGeo.recoilPos; orbitB.x = GetRecoilXPos(detGeo.recoilPos) ; diff --git a/Cleopatra/ClassTransfer.h b/Cleopatra/ClassTransfer.h index e433d2c..c4c38d5 100644 --- a/Cleopatra/ClassTransfer.h +++ b/Cleopatra/ClassTransfer.h @@ -30,9 +30,9 @@ public: targetA, targetZ, recoilA, recoilZ, beamEnergy_AMeV); } - TransferReaction(string configFile){ + TransferReaction(string configFile, unsigned short ID = 0){ Inititization(); - SetReactionFromFile(configFile); + SetReactionFromFile(configFile, ID); } ~TransferReaction(); @@ -49,7 +49,7 @@ public: void SetExA(double Ex); void SetExB(double Ex); - void SetReactionFromFile(string configFile); + void SetReactionFromFile(string configFile, unsigned short ID = 0); TString GetReactionName(); TString GetReactionName_Latex(); @@ -70,7 +70,8 @@ public: TLorentzVector GetPb(){return Pb;} TLorentzVector GetPB(){return PB;} - void PrintFourVectors(); + void PrintFourVectors() const; + void PrintReaction() const; void CalReactionConstant(); @@ -85,6 +86,7 @@ public: private: + Recoil recoil; ReactionConfig config; string nameA, namea, nameb, nameB; @@ -161,8 +163,8 @@ void TransferReaction::Seta(int A, int Z){ void TransferReaction::Setb(int A, int Z){ Isotope temp (A, Z); mb = temp.Mass; - config.recoil1.lightA = A; - config.recoil1.lightZ = Z; + recoil.lightA = A; + recoil.lightZ = Z; nameb = temp.Name; isReady = false; isBSet = false; @@ -170,8 +172,8 @@ void TransferReaction::Setb(int A, int Z){ void TransferReaction::SetB(int A, int Z){ Isotope temp (A, Z); mB = temp.Mass; - config.recoil1.heavyA = A; - config.recoil1.heavyZ = Z; + recoil.heavyA = A; + recoil.heavyZ = Z; nameB = temp.Name; isReady = false; isBSet = true; @@ -191,12 +193,14 @@ void TransferReaction::SetReactionSimple(int beamA, int beamZ, config.SetReactionSimple(beamA, beamZ, targetA, targetZ, - recoilA, recoilZ, beamEnergy_AMeV); + recoilA, recoilZ, beamEnergy_AMeV, 0); + + recoil = config.recoil[0]; SetA(config.beamA, config.beamZ); Seta(config.targetA, config.targetZ); - Setb(config.recoil1.lightA, config.recoil1.lightZ); - SetB(config.recoil1.heavyA, config.recoil1.heavyZ); + Setb(recoil.lightA, recoil.lightZ); + SetB(recoil.heavyA, recoil.heavyZ); SetIncidentEnergyAngle(config.beamEnergy, 0, 0); CalReactionConstant(); @@ -213,14 +217,17 @@ void TransferReaction::SetExB(double Ex){ isReady = false; } -void TransferReaction::SetReactionFromFile(string configFile){ +void TransferReaction::SetReactionFromFile(string configFile, unsigned short ID){ if( config.LoadReactionConfig(configFile) ){ SetA(config.beamA, config.beamZ); Seta(config.targetA, config.targetZ); - Setb(config.recoil1.lightA, config.recoil1.lightZ); - SetB(config.recoil1.heavyA, config.recoil1.heavyZ); + + recoil = config.recoil[ID]; + + Setb(recoil.lightA, recoil.lightZ); + SetB(recoil.heavyA, recoil.heavyZ); SetIncidentEnergyAngle(config.beamEnergy, 0, 0); CalReactionConstant(); @@ -260,9 +267,9 @@ TString TransferReaction::GetReactionName_Latex(){ void TransferReaction::CalReactionConstant(){ if( !isBSet){ - config.recoil1.heavyA = config.beamA + config.targetA - config.recoil1.lightA; - config.recoil1.heavyZ = config.beamZ + config.targetZ - config.recoil1.lightZ; - Isotope temp (config.recoil1.heavyA, config.recoil1.heavyZ); + recoil.heavyA = config.beamA + config.targetA - recoil.lightA; + recoil.heavyZ = config.beamZ + config.targetZ - recoil.lightZ; + Isotope temp (recoil.heavyA, recoil.heavyZ); mB = temp.Mass; isBSet = true; } @@ -278,11 +285,16 @@ void TransferReaction::CalReactionConstant(){ PA.RotateZ(phiIN); Pa.SetXYZM(0,0,0,ma); + + PA.SetUniqueID(config.beamZ); + Pa.SetUniqueID(config.targetZ); + Pb.SetUniqueID(recoil.lightZ); + PB.SetUniqueID(recoil.heavyZ); isReady = true; } -void TransferReaction::PrintFourVectors(){ +void TransferReaction::PrintFourVectors() const { 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()); @@ -297,6 +309,25 @@ void TransferReaction::PrintFourVectors(){ } +void TransferReaction::PrintReaction() const { + + printf("=====================================================\n"); + printf("------------------------------ Beam\n"); + printf(" beam : A = %3d, Z = %2d, Ex = %.2f MeV\n", config.beamA, config.beamZ, config.beamEx); + printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", config.beamEnergy, config.beamEnergySigma, config.beamEnergySigma/config.beamEnergy); + printf(" Angle : %.2f +- %.2f mrad\n", config.beamAngle, config.beamAngleSigma); + printf(" offset : (x,y) = (%.2f, %.2f) mmm \n", config.beamX, config.beamY); + + printf("------------------------------ Target\n"); + printf(" target : A = %3d, Z = %2d \n", config.targetA, config.targetZ); + + printf("------------------------------ Recoil\n"); + printf(" light : A = %3d, Z = %2d \n", recoil.lightA, recoil.lightZ); + printf(" heavy : A = %3d, Z = %2d \n", recoil.heavyA, recoil.heavyZ); + printf("=====================================================\n"); + +} + void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){ if( isReady == false ){ @@ -355,7 +386,7 @@ std::pair TransferReaction::CalExThetaCM(double e, double z, dou double mass = mb; double massB = mB; double y = e + mass; - double slope = 299.792458 * config.recoil1.lightZ * abs(Bfield) / TMath::TwoPi() * beta / 1000.; // MeV/mm; + double slope = 299.792458 * recoil.lightZ * 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/working/detectorGeo.txt b/working/detectorGeo.txt index e879c64..35730a3 100644 --- a/working/detectorGeo.txt +++ b/working/detectorGeo.txt @@ -13,6 +13,7 @@ false //is_coincident_with_recoil 0.00 //Elum_2_position_[mm]_when_Elum=0_disable_tree_branch #===============1st_Array +true ////is_this_array_exist_or_use_for_Simulation 11.5 //distance_from_axis_[mm] 10.0 //width_of_detector_[mm] 50 //length_of_detector_[mm] @@ -30,12 +31,12 @@ Out //detector_facing_Out_or_In 294.0 #===============2nd_Array -false //is_2nd_array_exist_is_use_for_Simulation +true //is_this_array_exist_or_use_for_Simulation 11.5 //distance_from_axis_[mm] 10.0 //width_of_detector_[mm] 50 //length_of_detector_[mm] 0 //blocker_length_[mm] -121 //first_position_-_for_upstream_[mm] +500 //first_position_-_for_upstream_[mm] 0.03 //energy_resolution_of_PSD_array_[MeV] 1.00 //position_resolution_of_PSD_array_[mm] Out //detector_facing_Out_or_In