diff --git a/CHANGELOG.md b/CHANGELOG.md index f7e7cfaa7..8c2c35b6b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ it in future. * Remove Millepede * Remove outdated example shipEvent_ex.py +* Remove ALPACA generator ## 24.11 diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 9aa006f51..82196a57f 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -3,7 +3,6 @@ import sys import ROOT ROOT.gSystem.Load('libEGPythia8') -import makeALPACAEvents import shipunit as u import shipRoot_conf @@ -86,7 +85,6 @@ parser.add_argument("--Genie", dest="genie", help="Genie for reading and processing neutrino interactions", required=False, action="store_true") parser.add_argument("--NuRadio", dest="nuradio", help="misuse GenieGenerator for neutrino radiography and geometry timing test", required=False, action="store_true") parser.add_argument("--Ntuple", dest="ntuple", help="Use ntuple as input", required=False, action="store_true") -group.add_argument("--ALPACA", dest="ALPACA", help="Use ALPACA as input", required=False, action="store_true") parser.add_argument("--MuonBack",dest="muonback", help="Generate events from muon background file, --Cosmics=0 for cosmic generator data", required=False, action="store_true") parser.add_argument("--FollowMuon",dest="followMuon", help="Make muonshield active to follow muons", required=False, action="store_true") parser.add_argument("--FastMuon", dest="fastMuon", help="Only transport muons for a fast muon only background estimate", required=False, action="store_true") @@ -155,7 +153,6 @@ if options.genie: simEngine = "Genie" if options.nuradio: simEngine = "nuRadiography" if options.ntuple: simEngine = "Ntuple" -if options.ALPACA: simEngine = "ALPACA" if options.muonback: simEngine = "MuonBack" if options.nuage: simEngine = "Nuage" if options.mudis: simEngine = "muonDIS" @@ -327,24 +324,6 @@ # P8gen.SetId(-211) primGen.AddGenerator(P8gen) -if simEngine == "ALPACA": - print('Generating ALP events of mass {} GeV with the photon coupling coefficient {} GeV^-1'.format(options.theMass, options.theDPepsilon)) - target = ship_geo.target - startZ = target.z0 - lengthZ = target.length - endZ = startZ + lengthZ - SmearBeam = 1*u.cm # finite beam size - Lmin = ((ship_geo.Chamber1.z - ship_geo.chambers.Tub1length) - ship_geo.target.z0)/100. - Lmax = (ship_geo.TrackStation1.z - ship_geo.target.z0)/100. - print('ALPACA is initializing.') - inputFile = makeALPACAEvents.runEvents(options.theMass,options.theDPepsilon,options.nEvents,Lmin,Lmax,startZ,endZ,SmearBeam) - if inputFile: print('ALPACA is done.') - ut.checkFileExists(inputFile) - ALPACAgen = ROOT.ALPACAGenerator() - ALPACAgen.Init(inputFile) - print('ALPACAGenerator is reading the ALPACA events') - primGen.AddGenerator(ALPACAgen) - if simEngine == "FixedTarget": P8gen = ROOT.FixedTargetGenerator() P8gen.SetTarget("volTarget_1",0.,0.) diff --git a/python/makeALPACAEvents.py b/python/makeALPACAEvents.py deleted file mode 100755 index 1caf58e70..000000000 --- a/python/makeALPACAEvents.py +++ /dev/null @@ -1,144 +0,0 @@ -import os,sys -import numpy,math -import ROOT - -hGeV = 6.58211928*pow(10.,-16)* pow(10.,-9) # no units or it messes up!! -c_light = 2.99792458e+10 -mres="1.0d0" -gax="1d-7" -nev="100" - -def ALPACAFormatting(s): - s=str(s) - if s.find('d')==-1: - s = numpy.format_float_scientific(float(s)) - s= s.replace('e','d') - s= s.replace('+','') - return s - - -def Ctau(mres,gax): - return c_light*hGeV*65.*math.pi/((float(gax)*float(gax))*(float(mres)*float(mres)*float(mres))) - -def Decaylength(e,p,ctau): - beta=p/e - gamma=e/math.sqrt(e*e-p*p) - return beta*gamma*ctau - -def Decayweight(Lmin,Lmax,Decaylength,LS): - return math.exp(-LS/Decaylength)*((Lmax-Lmin)/Decaylength) - -def inputWrite(mres,gax,nev,Lmin,Lmax):#need to apply Rmin & Rmax as well.. - Mres = ALPACAFormatting(mres) - Gax = ALPACAFormatting(gax) - L = ALPACAFormatting(float(Lmax)-float(Lmin)) - Lmin = ALPACAFormatting(Lmin) - Lmax = ALPACAFormatting(Lmax) - with open("input.DAT","w") as f: - f.write("****************************************************************************************\n") - f.write("*********** RE-RUN ./init IF FIRST FIVE PARAMETERS ARE CHANGED: **********************\n") - f.write("****************************************************************************************\n") - f.write("************************* Experimental parameters *************************************\n") - f.write("****************************************************************************************\n") - f.write("4d2 ! [ebeam] : Beam kinetic energy (GeV)\n") - f.write("'prot' ! [btype] : Beam type ('prot' = proton, 'elec' = electron)\n") - f.write("95.95d0 ! [aa] : Target mass number\n") - f.write("42d0 ! [az] : Target atomic number\n") - f.write("%s ! [lsh] : Shielding length (m)\n"%(Lmin)) - f.write("%s ! [dvol] : Decay volume (m)\n"%(L)) - f.write("1.5d-1 ! [rmin] : Inner radius (m)\n") - f.write("2.5d0 ! [rmax] : Outer radius (m)\n") - f.write("0.1d0 ! [dmin] : Minimum photon separation (m)\n") - f.write("1d0 ! [emin] : Minimum photon energy (GeV)\n") - f.write("1d-2 ! [athetamax] : Maximum ALP theta (when adecay is false)\n") - f.write("****************************************************************************************\n") - f.write("********************************** ALP parameters **************************************\n") - f.write("****************************************************************************************\n") - f.write("%s ! [mres] : ALP mass (GeV)\n"%(Mres)) - f.write("%s ! [gax] : ALP photon coupling (GeV^-1)\n"%(Gax)) - f.write("****************************************************************************************\n") - f.write("****************************** Output parameters **************************************\n") - f.write("****************************************************************************************\n") - f.write("'_%s_%s' ! [outtag] : For output file\n"%(mres,gax)) - f.write("****************************************************************************************\n") - f.write("************************* Integration parameters ***************************************\n") - f.write("****************************************************************************************\n") - f.write("100000 ! [ncall] : Number of calls for preconditioning\n") - f.write("10 ! [itmx] : Number of iterations for preconditioning\n") - f.write("1d0 ! [prec] : Relative accuracy (in %) in main run\n") - f.write("100000 ! [ncall1] : Number of calls in first iteration\n") - f.write("100000 ! [inccall] : Number of increase calls per iteration\n") - f.write("50 ! [itend] : Maximum number of iterations\n") - f.write("6 ! [iseed] : Random number seed (integer > 0)\n") - f.write("****************************************************************************************\n") - f.write("******************************* Unweighted events **************************************\n") - f.write("****************************************************************************************\n") - f.write(".true. ! [genunw] : Generate unweighted events\n") - f.write("%s ! [nev] : Number of events ( < 1000000 recommended)\n"%(nev)) - f.write("'hepevt' ! [erec] : Event record format ('lhe' = Les Houches, 'hepevt' = HEPEVT)\n") - f.write("****************************************************************************************\n") - f.write("******************************* Cuts ************************************************\n") - f.write(".true. ! [gencuts] : Generate cuts - [rmin], [rmax], [dmin], [emin]\n") - f.write(".true. ! [adecay] : Include ALP decay\n") - f.write("****************************************************************************************\n") - f.write("****************************************************************************************\n") - f.write("****************************************************************************************\n") - -def ntupleWrite(ctau,fn,Lmin,Lmax,startZ,endZ,SmearBeam,mres,gax,old): - Lmin = float(Lmin) - Lmax = float(Lmax) - startZ = float(startZ) - endZ = float(endZ) - SmearBeam = float(SmearBeam) - wdir = os.environ['ALPACABIN'] - with open(wdir+"/outputs/output_"+str(mres)+"_"+str(gax)+".dat","r") as fp: - for line in fp: - if "Cross section" in line: - line=line.split() - xs=line[3]+"_"+line[5] - fn=fn.readlines() - inputFile=old+"/"+"alp_m"+str(mres)+"_g"+str(gax)+"_xs"+str(xs)+".root" - f=ROOT.TFile(inputFile,"recreate") - ntup=ROOT.TNtuple("MCTrack","Track Informations","event:track:pdg:px:py:pz:x:y:z:parent:decay:e:tof:w") - for i,j in enumerate(range(5,len(fn),9)): - LS = ROOT.gRandom.Uniform(Lmin*100.,Lmax*100.) - zinter = ROOT.gRandom.Uniform(startZ,endZ) - dx, dy = ROOT.gRandom.Uniform(-1,+1)*SmearBeam, ROOT.gRandom.Uniform(-1,+1)*SmearBeam - tr=fn[j].split() - dau1=fn[j+1].split() - dau2=fn[j+2].split() - px,py,pz=float(tr[7]),float(tr[8]),float(tr[9]) - p = math.sqrt(px**2.+py**2.+pz**2.) - lam = LS/p - daux,dauy,dauz= dx+lam*px,dy+lam*py,zinter+lam*pz - dl=Decaylength(float(tr[10]),p,ctau) - w = Decayweight(Lmin*100.,Lmax*100.,dl,LS) - ntup.Fill(int(i),int(0),int(9900015),px,py,pz,dx,dy,zinter,int(-1),float(0),float(tr[10]),float(0),w) - #ntup.Fill(int(i),int(1),int(dau1[1]),float(dau1[7]),float(dau1[8]),float(dau1[9]),float(dau1[12])/10.+dx,float(dau1[13])/10.+dy,float(dau1[14])/10.+zinter,int(0),float(1),float(dau1[10]),float(dau1[15])/10./c_light,w) - #ntup.Fill(int(i),int(2),int(dau2[1]),float(dau2[7]),float(dau2[8]),float(dau2[9]),float(dau2[12])/10.+dx,float(dau2[13])/10.+dy,float(dau2[14])/10.+zinter,int(0),float(1),float(dau2[10]),float(dau2[15])/10./c_light,w) - ntup.Fill(int(i),int(1),int(dau1[1]),float(dau1[7]),float(dau1[8]),float(dau1[9]),daux,dauy,dauz,int(0),float(1),float(dau1[10]),float(dau1[15])/10./c_light,w) - ntup.Fill(int(i),int(2),int(dau2[1]),float(dau2[7]),float(dau2[8]),float(dau2[9]),daux,dauy,dauz,int(0),float(1),float(dau2[10]),float(dau2[15])/10./c_light,w) - ntup.Write() - f.Close() - return inputFile - -def runEvents(mres,gax,nev,Lmin,Lmax,startZ,endZ,SmearBeam): - print('ALPACA is starting for mass of {} GeV with photon coupling coeffiecient {} GeV^-1.'.format(mres,gax)) - pdg = ROOT.TDatabasePDG.Instance() - pdg.AddParticle('a','ALP', mres, False, gax, 0., 'a', 9900015) - wdir = os.environ['ALPACABIN'] - old=os.getcwd() - os.chdir(wdir) - inputWrite(mres,gax,nev,Lmin,Lmax) - rn="./alpaca < input.DAT" - os.system(rn) - print('ALPACA generated the events.') - pa="./evrecs/evrec_"+str(mres)+"_"+str(gax)+".dat" - fn=open(pa,"r") - ctau=Ctau(mres,gax) - print('Events are recording into a ntuple.') - inputFile = ntupleWrite(ctau,fn,Lmin,Lmax,startZ,endZ,SmearBeam,mres,gax,old) - os.chdir(old) - print('Ntuple is ready for the reading.') - return inputFile -#runEvents() diff --git a/shipgen/ALPACAGenerator.cxx b/shipgen/ALPACAGenerator.cxx deleted file mode 100644 index a1890f2b2..000000000 --- a/shipgen/ALPACAGenerator.cxx +++ /dev/null @@ -1,95 +0,0 @@ -#include -#include "TROOT.h" -#include "TFile.h" -#include "FairPrimaryGenerator.h" -#include "ALPACAGenerator.h" -#include "TDatabasePDG.h" // for TDatabasePDG -#include "TMath.h" // for Sqrt - -const Double_t cm = 10.; // pythia units are mm -const Double_t c_light = 2.99792458e+10; // speed of light in cm/sec (c_light = 2.99792458e+8 * m/s) - -using std::cout; -using std::endl; -// read events from ntuples produced - -// ----- Default constructor ------------------------------------------- -ALPACAGenerator::ALPACAGenerator() {} -// ------------------------------------------------------------------------- -// ----- Default constructor ------------------------------------------- -Bool_t ALPACAGenerator::Init(const char* fileName) { - return Init(fileName, 0); -} -// ----- Default constructor ------------------------------------------- -Bool_t ALPACAGenerator::Init(const char* fileName, const int firstEvent) { - cout << "Info ALPACAGenerator: Opening input file " << fileName << endl; - fInputFile = new TFile(fileName); - if (fInputFile->IsZombie()) { - cout << "-E ALPACAGenerator: Error opening the Signal file" << fileName << endl; - } - fTree = (TTree *)fInputFile->Get("MCTrack"); - fNevents = fTree->GetEntries(); - fn=firstEvent; - fTree->SetBranchAddress("event",&event); //event no. - fTree->SetBranchAddress("track",&track); //track no. - fTree->SetBranchAddress("pdg",&pdg); // particle pdg - fTree->SetBranchAddress("parent",&parent); // parent track number - fTree->SetBranchAddress("tof",&tof); // time of flight - fTree->SetBranchAddress("e",&e); // incoming muon energy - fTree->SetBranchAddress("w",&w); // weight of event - fTree->SetBranchAddress("x",&x); // position in x - fTree->SetBranchAddress("y",&y); // position in y - fTree->SetBranchAddress("z",&z); // position in z - fTree->SetBranchAddress("px",&px); // momentum in x - fTree->SetBranchAddress("py",&py); // momentum in y - fTree->SetBranchAddress("pz",&pz); // momentum in z - return kTRUE; -} -// ------------------------------------------------------------------------- - - -// ----- Destructor ---------------------------------------------------- -ALPACAGenerator::~ALPACAGenerator() -{ - // cout << "destroy Ntuple" << endl; - fInputFile->Close(); - fInputFile->Delete(); - delete fInputFile; -} -// ------------------------------------------------------------------------- - -// ----- Passing the event --------------------------------------------- -Bool_t ALPACAGenerator::ReadEvent(FairPrimaryGenerator* cpg) -{ - if (fnGetEntry(fn); // ALP is getting, - if (track==0&&parent==-1) { - cpg->AddTrack(pdg,px,py,pz,x,y,z,-1.,false,e,tof,w); // ALP track is added. - } - fTree->GetEntry(fn+1); // First daughter/photon is getting, - if (track==1&&parent==0) { - cpg->AddTrack(pdg,px,py,pz,x,y,z,0,true,e,tof,w); // First daughter/photon track is added. - } - fTree->GetEntry(fn+2); // Second daughter/photon is getting, - if (track==2&&parent==0) { - cpg->AddTrack(pdg,px,py,pz,x,y,z,0,true,e,tof,w); // Second daughter/photon track is added. - } - else { - return kFALSE; - } - fn=fn+3; // Goes to next event, more specifically next ALP, - cout << "Event Number:"<