SOLARIS_Analysis/Cleopatra/knockoutSim.C

#include "ClassHelio.h"
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TF1.h"
#include "TTree.h"
#include "TRandom.h"
#include "TClonesArray.h"
#include <vector>
#include <fstream>

//----------- usage 
// $root sim.C+ | tee output.txt
// this will same the massage to output.txt

void knockout(){

   //================================================= User Setting
   //---- reaction
   int AA = 23, ZA = 9;
   int Aa = 1,  Za = 1;
   int A2 = 1,  Z2 = 1;
   
   bool isNormalKinematics = false;
   bool isOverRideExNegative = true;
   double maxkb = 200.;
   
   //---- beam
   double KEAmean = 100; // MeV/u
   const int nKEA = 1; 
   double KEAList[nKEA] = {300}; 
   
   double KEAsigma = 0; //KEAmean*0.001; // MeV/u , assume Guassian
   double thetaMean = 0.; // mrad 
   double thetaSigma = 0.; // mrad , assume Guassian due to small angle
   
   int numEvent = 100000;
   
   //---- HELIOS detector geometry
   string heliosDetGeoFile = "";//"detectorGeo_upstream.txt";
   double BField = 4.0; // if not detector, must set B-field, else, this value is not used.
   double BFieldTheta = 0.; // direction of B-field
   double eSigma = 0.0001 ; // detector energy sigma MeV
   double zSigma = 0.1 ; // detector position sigma mm
   
   //---- excitation of Beam 
   int nExA = 1;
   double ExAList[nExA];
   ExAList[0] = 0.000; // MeV
   //ExAList[1] = 1.567;
    
   //---- Separation energy
   string separationFile = "separation_energies.txt";
   
   //---- save root file name
   TString saveFileName = "knockout.root";
   
   //---- Auxiliary setting
   bool isTargetScattering = false;
   bool isDecay = false;
   bool isReDo = false; // redo calculation until detected. 
   
   //---- target
   double density = 0.913; // 0.913 g/cm3
   double targetThickness = 2./2. * 2.2e-4; // 2.2 um = 201 ug/cm2
   string stoppingPowerForA = "208Pb_in_CD2.txt"; // generated by SRIM
   string stoppingPowerForb = "1H_in_CD2.txt";
   string stoppingPowerForB = "209Pb_in_CD2.txt";
   
   //=============================================================
   //=============================================================
   //=============================================================
   //===== Set Reaction
   Knockout reaction;
   int AB = AA-A2, ZB = ZA-Z2;
   int A1 = Aa   , Z1 = Za;
   reaction.SetA(AA,ZA);
   reaction.Seta(Aa,Za);
   reaction.Set2(A2,Z2);
   reaction.SetIncidentEnergyAngle(KEAmean, 0, 0);
   reaction.OverRideExNegative(isOverRideExNegative);
   
   printf("===================================================\n");
   printf("=========== %s ===========\n", reaction.GetReactionName().Data());
   printf("=========== KE: %9.4f +- %5.4f MeV/u, dp/p = %5.2f %% \n", KEAmean, KEAsigma, KEAsigma/KEAmean * 50.);
   printf("======== theta: %9.4f +- %5.4f MeV/u \n", thetaMean, thetaSigma);
   printf("===================================================\n");

   
   //======== Set HELIOS
   printf("############################################## HELIOS configuration\n");   
   HELIOS helios1; // for particle-1
   HELIOS helios2; // for particle-2
   helios1.SetMagneticFieldDirection(BFieldTheta);
   helios2.SetMagneticFieldDirection(BFieldTheta);
   bool sethelios1 = helios1.SetDetectorGeometry(heliosDetGeoFile);
   bool sethelios2 = helios2.SetDetectorGeometry(heliosDetGeoFile);
   if( sethelios1 && sethelios2 ) {
		int mDet = helios1.GetNumberOfDetectorsInSamePos();
		printf("========== energy resol.: %f MeV\n", eSigma);
		printf("=========== pos-Z resol.: %f mm \n", zSigma); 
	}else{
		helios1.SetMagneticField(BField);
		helios2.SetMagneticField(BField);
		printf("======== B-field : %5.2f T, Theta: %5.2f deg\n", BField, BFieldTheta);
	}
	
	//==== Target scattering, only energy loss
   if(isTargetScattering) printf("############################################## Target Scattering\n");
   TargetScattering msA;
   TargetScattering msB;
   TargetScattering msb;

   if(isTargetScattering) printf("======== Target : (thickness : %6.2f um) x (density : %6.2f g/cm3) = %6.2f ug/cm2\n", 
                        targetThickness * 1e+4, 
                        density, 
                        targetThickness * density * 1e+6);  
            
   if( isTargetScattering ){
      msA.LoadStoppingPower(stoppingPowerForA);
      msb.LoadStoppingPower(stoppingPowerForb);
      msB.LoadStoppingPower(stoppingPowerForB);
   }
   
   //======= Decay of particle-B
   Decay decay;
   decay.SetMotherDaugther(AB, ZB, AB-1,ZB); //neutron decay
	
	
   //======= loading excitation energy
   printf("############################################## excitation energies\n");
   vector<double> SpList;
   printf("----- loading separation energies.");
   ifstream file;
   file.open(separationFile.c_str());
   string isotopeName;
   if( file.is_open() ){
      string line;
      int i = 0;
      while( file >> line){
         //printf("%d, %s \n", i,  line.c_str());
         if( line.substr(0,2) == "//" ) continue;
         if( i == 0 ) isotopeName = line; 
         if ( i >= 1 ){
            SpList.push_back(atof(line.c_str()));
         }
         i = i + 1;
      }
      file.close();
      printf("... done.\n");
      printf("========== %s\n", isotopeName.c_str());
      int n = SpList.size();
      for(int i = 0; i < n ; i++){
         if( isDecay ) {
            TLorentzVector temp(0,0,0,0);
            int decayID = decay.CalDecay(temp, SpList[i], 0);
            if( decayID == 1) {
               printf("%d, Sp: %6.2f MeV --> Decay. \n", i, SpList[i]);
            }else{
               printf("%d, Sp: %6.2f MeV\n", i, SpList[i]);
            }
         }else{
            printf("%d, Sp: %6.2f MeV \n", i, SpList[i]);
         }
      }
   }else{
       printf("... fail\n");
       return;
   }

   //====================== build tree
   TFile * saveFile = new TFile(saveFileName, "recreate");
   TTree * tree = new TTree("tree", "tree");
   
   double thetaNN, phiNN;
   double theta1, phi1, T1;
   double theta2, phi2, T2;
   double thetaB, TB;
   double Sp, kb, thetab, phib;
   int SpID;

   double ExA;
   int ExAID;
   double KEA, KEAscattered, theta, phi;
   
   double mB,mb;

   tree->Branch("theta1", &theta1, "theta1/D");
   tree->Branch("phi1", &phi1, "phi1/D");
   tree->Branch("T1", &T1, "T1/D");
   tree->Branch("theta2", &theta2, "theta2/D");
   tree->Branch("phi2", &phi2, "phi2/D");
   tree->Branch("T2", &T2, "T2/D");
   tree->Branch("thetaB", &thetaB, "thetaB/D");
   tree->Branch("TB", &TB, "TB/D");
   tree->Branch("thetaNN", &thetaNN, "thetaNN/D");
   tree->Branch("phiNN", &phiNN, "phiNN/D");
   tree->Branch("Sp", &Sp, "Sp/D");
   tree->Branch("kb", &kb, "kb/D");
   tree->Branch("thetab", &thetab, "thetab/D");
   tree->Branch("phib", &phib, "phib/D");
   
   tree->Branch("SpID", &SpID, "SpID/I");
   
   tree->Branch("ExAID", &ExAID, "ExAID/I");
   tree->Branch("KEA", &KEA, "KEA/D");
   if(isTargetScattering) tree->Branch("KEAscattered", &KEAscattered, "KEAscattered/D");
   tree->Branch("theta", &theta, "theta/D");
   tree->Branch("phi", &phi, "phi/D");

   tree->Branch("mB", &mB, "mB/D");
   tree->Branch("mb", &mb, "mb/D");
   
   tree->Branch("Bfield", &BField, "Bfield/D");
   
   TClonesArray * arr = new TClonesArray("TLorentzVector");
   tree->Branch("fV", arr, 256000);
   arr->BypassStreamer();
   
   TClonesArray * arrN = new TClonesArray("TLorentzVector");
   tree->Branch("fVN", arrN, 256000);
   arrN->BypassStreamer();
   
   TLorentzVector* fourVector = NULL;

   // the output of Helios.CalHit
   double e1, z1, t1, rho1, x1h, y1h, r1h;
   double e2, z2, t2, rho2, x2h, y2h, r2h;
	tree->Branch("e1", &e1, "e1/D");
   tree->Branch("z1", &z1, "z1/D");
   tree->Branch("t1", &t1, "t1/D");
   tree->Branch("x1h", &x1h, "x1h/D"); //at 200 mm downstream
   tree->Branch("y1h", &y1h, "y1h/D");
   tree->Branch("r1h", &r1h, "r1h/D");
   tree->Branch("rho1", &rho1, "rho1/D");
   
   tree->Branch("e2", &e2, "e2/D");
   tree->Branch("z2", &z2, "z2/D");
   tree->Branch("t2", &t2, "t2/D");
   tree->Branch("x2h", &x2h, "x2h/D"); //at 200 mm downstream
   tree->Branch("y2h", &y2h, "y2h/D");
   tree->Branch("r2h", &r2h, "r2h/D");
   tree->Branch("rho2", &rho2, "rho2/D");
   
   double Sp2;
   tree->Branch("Sp2", &Sp2, "Sp2/D");
   
   //different coordinate
   double a0, a1, a2; // in k1-k2 coordinate
   tree->Branch("a0", &a0, "a0/D");
   tree->Branch("a1", &a1, "a1/D");
   tree->Branch("a2", &a2, "a2/D");
   
   double b; // in THREEDEE coordinate
   tree->Branch("b", &b, "b/D");
   
   //========timer
   TBenchmark clock;
   bool shown ;   
   clock.Reset();
   clock.Start("timer");
   shown = false;
   printf("############################################## generating %d events \n", numEvent);
   
   //====================================================== calculate 
   int count = 0;
   for( int i = 0; i < numEvent; i++){
      bool redoFlag = true;
      if( !isReDo ) redoFlag = false;
      do{
      
         //==== Set Ex of A
         ExAID = gRandom->Integer(nExA);
         ExA = ExAList[ExAID];
         reaction.SetExA(ExA);
         
         //==== Set Sp of B
         SpID = gRandom->Integer(SpList.size());
         Sp = SpList[SpID]; 
         
         //==== Set incident beam
         if( KEAsigma == 0 ){
            KEA = KEAmean;
         }else{
            KEA = gRandom->Gaus(KEAmean, KEAsigma);
         }
         if( thetaSigma == 0 ){
            theta = thetaMean;
         }else{
            theta = gRandom->Gaus(thetaMean, thetaSigma);
         }
         
         int rKEA = gRandom->Integer(nKEA);
         KEA = KEAList[rKEA];
         
         /*
         KEA = 300*gRandom->Rndm();
         BField = 4*gRandom->Rndm();
			helios1.SetMagneticField(BField);
			helios2.SetMagneticField(BField);
			*/
         
         reaction.SetIncidentEnergyAngle(KEA, theta, 0.);
         
         
         /*
         //For target scattering
         reaction.CalIncidentChannel(isNormalKinematics); // but only need is PA
         TLorentzVector PA = reaction.GetPA();            
         
         double depth = 0;
         if( isTargetScattering ){
            //==== Target scattering, only energy loss
            depth = targetThickness * gRandom->Rndm();
            msA.SetTarget(density, depth); 
            TLorentzVector PAnew = msA.Scattering(PA);
            KEAscattered = msA.GetKE()/AA;
            reaction.SetIncidentEnergyAngle(KEAscattered, theta, phi);
         }*/
         
         //==== Calculate reaction
         thetaNN = TMath::ACos(2 * gRandom->Rndm() - 1) ; 
         phiNN = TMath::TwoPi() * gRandom->Rndm(); 
         
         kb = maxkb * gRandom->Rndm();
         thetab = TMath::ACos(2 * gRandom->Rndm() - 1);
         phib = TMath::TwoPi() * gRandom->Rndm();

         reaction.SetBSpk(Sp, kb, thetab, phib);
         reaction.CalReactionConstant(isNormalKinematics);
         reaction.Event(thetaNN, phiNN);

         TLorentzVector PA = reaction.GetPA();
         TLorentzVector Pa = reaction.GetPa();
               
         TLorentzVector P1 = reaction.GetP1();
         TLorentzVector P2 = reaction.GetP2();
         TLorentzVector Pb = reaction.GetPb();
         TLorentzVector PB = reaction.GetPB();
         
         /*
         //==== Calculate energy loss of scattered and recoil in target
         if( isTargetScattering ){
            if( Pb.Theta() < TMath::PiOver2() ){
               msb.SetTarget(density, targetThickness - depth);
            }else{
               msb.SetTarget(density, depth);
            }
            Pb = msb.Scattering(Pb);
            TbLoss = msb.GetKELoss();
            msB.SetTarget(density, targetThickness - depth);
            PB = msB.Scattering(PB);
         }else{
            TbLoss = 0;
         }
         
         //======= Decay of particle-B
         if( isDecay){
            int decayID = decay.CalDecay(PB, Sp, 0); // decay to ground state
            if( decayID == 1 ){
               PB = decay.GetDaugther_D();
               decayTheta = decay.GetAngleChange();
            }else{
               decayTheta = TMath::QuietNaN();
            }
         }
         */
                  
         //################################### tree branches
         //===== reaction
         theta1 = P1.Theta() * TMath::RadToDeg();
         theta2 = P2.Theta() * TMath::RadToDeg();
         thetaB = PB.Theta() * TMath::RadToDeg();
      
         T1 = P1.E() - P1.M();
         T2 = P2.E() - P2.M();
         TB = PB.E() - PB.M();
         
         phi1 = P1.Phi() * TMath::RadToDeg();
         phi2 = P2.Phi() * TMath::RadToDeg();
         
         //--------- diff coordinate
         b = TMath::ASin( TMath::Sin(P2.Theta()) * TMath::Sin(P1.Phi() - P2.Phi() )) * TMath::RadToDeg();
         
         TVector3 kA = (PA.Vect()).Unit();
         TVector3 k1 = (P1.Vect()).Unit();
         TVector3 k2 = (P2.Vect()).Unit();
         TVector3 n = (k1.Cross(k2)).Unit();
         TVector3 j = (kA - (kA.Dot(n))*n).Unit();
         
         a0 = TMath::ASin(n.Dot(kA)) * TMath::RadToDeg();
         a1 = TMath::ACos(k1.Dot(j)) * TMath::RadToDeg();
         a2 = TMath::ACos(k2.Dot(j)) * TMath::RadToDeg();
         
         //----------- mass
         mB = PB.M();
         mb = Pb.M();
         
         TVector3 bA = PA.BoostVector();
            
         for(int i = 0; i < 6 ; i++){
            TLorentzVector temp;
            double xyzt[4];
            switch(i){
               case 0: temp = PA; break;
               case 1: temp = Pa; break;
               case 2: temp = P1; break;
               case 3: temp = P2; break;
               case 4: temp = PB; break;
               case 5: temp = Pb; break;
            }
            
            temp.GetXYZT(xyzt);
            
            fourVector = (TLorentzVector*) arr->ConstructedAt(i);
            fourVector->SetXYZT(xyzt[0], xyzt[1], xyzt[2], xyzt[3]);
            
            //into normal kinematic
            temp.Boost(-bA);  
            temp.GetXYZT(xyzt);
            
            fourVector = (TLorentzVector*) arrN->ConstructedAt(i);
            fourVector->SetXYZT(xyzt[0], xyzt[1], xyzt[2], xyzt[3]);
            
         }
         
         
         //==== Helios
         int hit1 = helios1.CalHit(P1, Z1, PB, ZB);
         int hit2 = helios2.CalHit(P2, Z2, PB, ZB);
         
         double recoilZ = 1000;
         
         e1 = helios1.GetEnergy() + gRandom->Gaus(0, eSigma);
         z1 = helios1.GetZ() ; 
         t1 = helios1.GetTime();
         rho1 = helios1.GetRho();
         x1h = helios1.GetXPos(recoilZ);
         y1h = helios1.GetYPos(recoilZ);
         r1h = helios1.GetR(recoilZ);
         
         e2 = helios2.GetEnergy() + gRandom->Gaus(0, eSigma);
         z2 = helios2.GetZ() ; 
         t2 = helios2.GetTime();
         rho2 = helios2.GetRho();
         x2h = helios2.GetXPos(recoilZ);
         y2h = helios2.GetYPos(recoilZ);
         r2h = helios2.GetR(recoilZ);
         
         double gammaA = PA.Gamma();
         double betaA = PA.Beta();
         
         Sp2 = (1-gammaA)* 938.272 - gammaA*(e1+e2) + betaA*gammaA*(P1.Pz() + P2.Pz()) ;
         
         //printf("%f, %f | %f, %f \n", e1, z1, e2, z2);
         
         //change thetaNN into deg
         thetaNN = thetaNN * TMath::RadToDeg();
         
         if( hit1 == 1) {
            count ++;
         }
         
         if( isReDo ){
            if( hit1 == 1) {
               redoFlag = false;
            }else{
               redoFlag = true;
               //printf("%d, %2d, thetaNN : %f, theta : %f, z0: %f \n", i, hit, thetaNN * TMath::RadToDeg(), thetab, helios.GetZ0());
            }
         }else{
            redoFlag = false;
         }
         
         
      }while( redoFlag );
      tree->Fill();
      
      //#################################################################### Timer  
      clock.Stop("timer");
      Double_t time = clock.GetRealTime("timer");
      clock.Start("timer");
      
      if ( !shown ) {
         if (fmod(time, 10) < 1 ){
            printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
            shown = 1;
         }
      }else{
         if (fmod(time, 10) > 9 ){
            shown = 0;
         }
      }
   }
   
   saveFile->Write();
   saveFile->Close();
   
   printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
   
   gROOT->ProcessLine(".q");
   
   /**/

}