madx2g4bl

 function [outtext]=madx2g4bl(P_0,particle,folder)
 tansform madx file into G4BL input file.

 Simone Maria Liuzzo PhD@LNF 11-jul-2011

0001 % function [outtext]=madx2g4bl(P_0,particle,folder)
0002 % tansform madx file into G4BL input file.
0003 %
0004 % Simone Maria Liuzzo PhD@LNF 11-jul-2011
0005 
0006 format longe
0007 
0008 StepLength=10; % [mm] default is 10mm
0009 
0010 Nevents=1; % e+ events and e- events.
0011 disp(['writing G4BL file to simulate: ' num2str(Nevents) ' e+ and ' num2str(Nevents) ' e-'])
0012 det=0;
0013 SR=0;
0014 
0015 TotLength=1258358.149; %[mm]
0016 
0017 %beamstart=TotLength/2;
0018 %beamstart=581.981242*1000-140; % first QDI center
0019 
0020 beamstart=0;
0021 
0022 % beams shapes
0023 
0024 % diversi da zero muore!!! sestupoli??? KL? BRHO?
0025 sX=1*1e-3; %[mm]
0026 sY=1*1e-3; %[mm]
0027 sXp=1*1e-6; %[rad]
0028 sYp=1*1e-6; %[rad]
0029 
0030 killval=1;
0031 
0032 Xang=0.06600162900000006; % rad crossing angle
0033 
0034 
0035 %%%% HER %%%%%
0036 %%%%%%%%%%%%%%
0037 
0038 folder='SBHERV12DATA';
0039 P_0=6.699999981; % GeV/c
0040 particle='e+';
0041 
0042 sxt=1;
0043 quadH=1;  % go around in the other direction.
0044 
0045 Brho=3.3356409519815205*P_0; % positron and opposite rotation
0046 
0047 
0048 bigmondo=[...
0049     ...'box mondo height=10000 width=3000000 length=3000000 material=Vacuum \n'...
0050     ];
0051 
0052 posbeam=[ ' beam gaussian particle=' particle ...
0053     ' beamZ=' num2str(beamstart) ...
0054     ' beamX=' num2str(+1*1.532845188e-07,15) ...
0055     ' beamXp=' num2str(+1*1.683548689e-07,15) ...
0056     ' rotation=Y' num2str(0*180)... rotate to head back
0057     '  nEvents=' num2str(Nevents) ' '... % nEvents=10
0058     '    sigmaX=' num2str(sX,15) ...
0059     ' sigmaY=' num2str(sY,15)...
0060     ' sigmaXp=' num2str(sXp,15)...
0061     ' sigmaYp=' num2str(sYp,15) ' '...
0062     '    meanMomentum=' num2str(P_0*1000,15) ' sigmaP=0'...
0063     ' meanT=0 sigmaT=0 \n'...
0064     'trackcolor reference=1,1,1 \n'...
0065     ];
0066 
0067 head=['physics QGSP synchrotronRadiation=' num2str(SR) '\n'...
0068     'start x=0 y=0 z=' num2str(beamstart) ' radiusCut=0.001 ring=1\n'...
0069     'reference particle=' particle ' referenceMomentum=' num2str(P_0*1000,15) ' noEloss=1'...
0070     ' beamXp=0 beamYp=0 '...
0071     ' tuneMomentum=' num2str(P_0,15) '\n'...
0072     'box BeamVis width=200.0 height=200.0 length=0.1 material=Vacuum color=1,0,0 kill=1 \n'...
0073     'trace nTrace=' num2str(2*Nevents) ' oneNTuple=1 format=ascii \n'...
0074     ' param eventTimeLimit=60*8\n'... %  8 minutes maximum time for one simulation
0075     ' trackcuts keep=e+,e-,gamma maxTime=1000000*40\n' ... %ns --> 4micros
0076     ];
0077 % x damping time= 0.26700741E-01 => 15% =~ 40ms = 40*10^6 ns
0078 
0079 % get pos K and L for Quad, sext and bend
0080 Qfile=fopen([folder '/quad']);
0081 Sfile=fopen([folder '/sext']);
0082 Bfile=fopen([folder '/bend']);
0083 % get pos K and L for Quad, sext and bend -- sbend and not rbend.
0084 %Qfile=fopen('DIAMDATA/quad');
0085 %Sfile=fopen('DIAMDATA/sext');
0086 %Bfile=fopen('DIAMDATA/bend');
0087 
0088 % name s K L
0089 Q=textscan(Qfile,'%s %f %f %f','Headerlines',47);
0090 B=textscan(Bfile,'%s %f %f %f','Headerlines',47);
0091 S=textscan(Sfile,'%s %f %f %f','Headerlines',47);
0092 
0093 defQ= ['virtualdetector Det radius=10.0 length=' num2str(2*StepLength)...
0094     ' color=0,.2,.0 referenceParticle=1 format=rootExtended\n'];
0095 defS='';
0096 defB='';
0097 ord=1;
0098 plac={};
0099 defined={};
0100 d=1;
0101 
0102 for i=1:length(Q{1})
0103     name=(Q{1}(i));
0104     name{1}(1)='';
0105     name{1}(end)='';
0106     Q{1}(i)=name;
0107 end
0108 for i=1:length(S{1})
0109     name=S{1}(i);
0110     name{1}(1)='';
0111     name{1}(end)='';
0112     S{1}(i)=name;
0113 end
0114 for i=1:length(B{1})
0115     name=B{1}(i);
0116     name{1}(1)='';
0117     name{1}(end)='';
0118     B{1}(i)=name;
0119 end
0120 
0121 % small final focus quad
0122 ffQuad{1}=char(Q{1}(1));%{'QD0AL'};
0123 ffQuad{2}=char(Q{1}(2));%{'QD0BL'};
0124 ffQuad{3}=char(Q{1}(3));%{'QD0AR'};
0125 ffQuad{4}=char(Q{1}(4));%{'QD0AR'};
0126 ffQuad{5}=char(Q{1}(end));%{'QD0BR'};
0127 ffQuad{6}=char(Q{1}(end-1));%{'QDPMR'};
0128 ffQuad{7}=char(Q{1}(end-2));%{'QDPML'};
0129 ffQuad{8}=char(Q{1}(end-3));%{'QD0BR'};
0130 
0131 
0132 plac{ord}=[' cornerarc' ' z=0+' num2str(0) '+0+' num2str(0,15) ...
0133     ' angle=' num2str(-Xang/2*180/pi,15)  ' centerRadius='  num2str(1,15) '\n'...
0134     ];
0135 ord=ord+1;
0136 
0137 for i=1:length(Q{1})
0138     %genericquad QDpi fieldLength=280 apertureRadius=101.5 ironRadius=381 \
0139     %    ironLength=250 ironColor=0.6,.0,.6 kill=' num2str(killval) ' gradient=$KQDpi
0140     %
0141     if sum(strcmp(char(Q{1}(i)),defined))==0;
0142         inrad=60;
0143         outrad=70;
0144         
0145         if sum(strcmp(char(Q{1}(i)),ffQuad))~=0
0146             inrad=12;
0147             outrad=15;
0148             tempSL=StepLength;
0149             StepLength=0.01;
0150             if strcmp(char(Q{1}(i)),ffQuad{1}) || strcmp(char(Q{1}(i)),ffQuad{5}) %% P.M.
0151                 inrad=8;
0152                 outrad=10;
0153                 
0154             end
0155             
0156         end
0157         defQ=[ defQ ...
0158             ' genericquad ' Q{1}(i) ' '...
0159             ' fieldLength=' num2str(Q{4}(i)*1000,15) ... % mm
0160             ' apertureRadius=' num2str(inrad,15) ' '...
0161             ' ironRadius=' num2str(outrad,15) ' '...
0162             ' ironLength=' num2str(Q{4}(i)*1000,15) ' '...
0163             ' ironColor=0/255,191/255,255/255 kill=' num2str(killval) ' '...
0164             ' fringe=0'...
0165             ' gradient=' num2str(quadH*Q{3}(i)*Brho/Q{4}(i),15) ' '... K
0166             ' maxStep=' num2str(StepLength/2)  ...
0167             '\n']; %#ok<*AGROW>
0168         
0169         % add to already defined
0170         defined{d}=char(Q{1}(i));
0171         d=d+1;
0172     end
0173     
0174     StepLength=tempSL;
0175     
0176     %defQ=[defQ 'virtualdetector Det radius=400.0 color=0,.1,.1'];
0177     
0178     if (ord/2-floor(ord/2))==(det-1) % set 0 to place a detector in every quadrupole.
0179         plac{ord}=['place ' Q{1}(i) ' parent=  z=' num2str(Q{2}(i)*1000,15)...
0180             ' coordinates=centerline'...
0181             '\n'...
0182             'place ' 'Det' ...
0183             ' rename=DET# '...
0184             ' parent=  z=' num2str((Q{2}(i))*1000)...
0185             ' coordinates=centerline '...
0186             '\n']; %#ok<*SAGROW>
0187     else
0188         plac{ord}=['place ' Q{1}(i) ' parent=  z=' num2str(Q{2}(i)*1000,15)...
0189             ' coordinates=centerline'...
0190             '\n']; %#ok<*SAGROW>
0191     end
0192     %      if sum(strcmp(Q{1}(i),{'"QD0AL"';'"QD0BL"';'"QD0AR"';'"QD0BR"'}))>0
0193     %
0194     %      plac{ord}=['place ' Q{1}(i) ' parent=  z=' num2str(Q{2}(i)*1000,15)...
0195     %      ' coordinates=centerline'...
0196     %      ' rotation=Y' num2str(-0.066*180/pi,15) ' '...
0197     %      '\n'];
0198     %     end
0199     ord=ord+1;
0200 end
0201 
0202 for i=1:length(S{1})
0203     %multipole SF1 fieldLength=150 apertureRadius=101.5 ironRadius=381 \
0204     %    ironLength=140 ironColor=0,.6,.6 kill=' num2str(killval) ' sextupole=$KSF1
0205     
0206     if sum(strcmp(char(S{1}(i)),defined))==0;
0207         inrad=80;
0208         outrad=120;
0209         
0210         defS=[ defS ...
0211             ' multipole ' S{1}(i) ' '...
0212             ' fieldLength=' num2str(S{4}(i)*1000,15) ... % mm
0213             ' apertureRadius=' num2str(inrad,15) ' '...
0214             ' ironRadius=' num2str(outrad,15) ' '...
0215             ' ironLength=' num2str(S{4}(i)*1000,15) ' '...
0216             ' ironColor=173/255,255/255,47/255 kill=' num2str(killval) ' '...
0217             ' fringe=0'...
0218             ' maxStep=' num2str(StepLength/2)  ...
0219             ' sextupole=' num2str(sxt*S{3}(i)*Brho/S{4}(i),15) ' '... K2
0220             '\n'];
0221         % add to already defined
0222         defined{d}=char(S{1}(i));
0223         d=d+1;
0224     end
0225     
0226     
0227     plac{ord}=['place ' S{1}(i) ' parent=  z=' num2str(S{2}(i)*1000,15)...
0228         ' coordinates=centerline '...
0229         '\n']; %#ok<*SAGROW>
0230     ord=ord+1;
0231     
0232 end
0233 
0234 for i=1:length(B{1})
0235     
0236     %genericbend BHer fieldWidth=1660 fieldHeight=200 fieldLength=5400 \
0237     %    ironColor=1,0,0 ironWidth=1828 ironHeight=1320 ironLength=5000 \
0238     % kill=' num2str(killval) ' By=0.2778157677856496
0239     %tune B1 z0=100 z1=3000 initial=-0.6500 step=0.01 expr=Px1/Pz1 \
0240     %tolerance=0.000001
0241     %genericbend B fieldWidth=500 fieldHeight=500 fieldLength=500 \
0242     %ironWidth=800 ironHeight=800 ironLength=500 \
0243     %fieldMaterial=Vacuum place B z=2000 rename=B1 rotation=Y15 By=B1
0244     
0245     if sum(strcmp(char(B{1}(i)),defined))==0 % if not defined, define.
0246         Fw=240*1.0;
0247         Fh=70;
0248         Iw=10+Fw;
0249         Ih=30+Fh;
0250         
0251         split=1; %make 'split' parts of a single magnet
0252         defB=[ defB ...
0253             '\n param B' num2str(i) B{1}(i) '=' num2str(-Brho*B{3}(i)/B{4}(i),15) '\n'...
0254             ];
0255         if B{2}(i)>140 && B{2}(i)<1140
0256             
0257             defB=[ defB ...
0258                 %  ' tune B' num2str(i) B{1}(i) ' z0=' num2str(B{2}(i)*1000-B{4}(i)*500) ' z1=' num2str(B{2}(i)*1000+B{4}(i)*500) ' '...
0259                 %  ' initial=' num2str(Brho*B{3}(i)/B{4}(i),15) ' step=1e-13'...
0260                 %  ' expr=(Px1/Pz1) tolerance=0.0000001 \n'...
0261                 ];
0262         end
0263         defB=[ defB ...
0264             'genericbend ' B{1}(i) ' '...
0265             ' fieldLength=' num2str(B{4}(i)*1000/split,15) ... % mm
0266             ' fieldWidth=' num2str(Fw,15) ' '...
0267             ' fieldHeight=' num2str(Fh,15) ' '...
0268             ' ironLength=' num2str(B{4}(i)*1000/split,15) ... % mm
0269             ' ironWidth=' num2str(Iw,15) ' '...
0270             ' ironHeight=' num2str(Ih,15) ' '...
0271             ' fringe=0'...
0272             ' maxStep=' num2str(StepLength)  ...
0273             ... ' By=' num2str(Brho*B{3}(i)/B{4}(i),15)...' fieldColor=0.5,.1,.0 '...
0274             ' By=$B' num2str(i) B{1}(i) ...
0275             ' ironColor=1,.0,.0 kill=' num2str(killval) ' '...
0276             '\n'...
0277             ];
0278         % add to already defined
0279         defined{d}=char(B{1}(i));
0280         d=d+1;
0281         
0282     end
0283     
0284     plac{ord}=[
0285         'place ' B{1}(i) ' parent=  z=0+' num2str(0) '+0+' num2str(B{2}(i)*1000,15) ...
0286         ' x=' num2str(B{4}(i)/2*tan(B{3}(i)/2)*1000) ''... % move in by a sagitta amount
0287         ' coordinates=centerline'...' \n'...
0288         ' rotation=Y' num2str(B{3}(i)*180/pi/2,15)...' rename=' B{1}(i) ' \n'...
0289         ...' By=B' num2str(i) B{1}(i) ...
0290         '\n' ...
0291         ' cornerarc' ' z=0+' num2str(0) '+0+' num2str((B{2}(i)-B{4}(i)/2)*1000,15) ...
0292         ' angle=' num2str(B{3}(i)*180/pi,15)  ' centerRadius='  num2str((B{4}(i)/B{3}(i))*1000,15) '\n'...'#corner' ' parent=  z=' num2str((B{2}(i)+B{4}(i)/2)*1000,15) ... ' rotation=Y' num2str(B{3}(i)*180/pi/2,15) ' \n'...
0293         '\n']; %#ok<*SAGROW>
0294     ord=ord+1;
0295     
0296 end
0297 
0298 RF_on=0;
0299 % plac RFC
0300 %length[m]        voltage[MV]                lag          freq[MHz]             harmon
0301 %0.501777                0.6           0.424389        476.0054439               1998
0302 RF=['pillbox RFC maxGradient=' num2str(RF_on*0.6/0.501777,15) ...
0303     ' color=0.5,0.5,0 frequency=0.4760054439 innerLength=501.777 '...
0304     ' phaseAcc=0.424389 wallThick=100 innerRadius=1000\n'];
0305 RF1pos=577.1465514*1000;
0306 RF2pos=579.1935514*1000;
0307 plac{ord}=['  # place RFC z=' num2str(RF1pos) '\n'];
0308 plac{ord+1}=['  # place RFC z=' num2str(RF2pos) '\n'];
0309 
0310 ippos=0;
0311 Dumppos=beamstart-(+11+280);
0312 
0313 positions=[0;Q{2}; S{2}; B{2}; (RF1pos-TotLength)/1000; (RF2pos-TotLength)/1000 ;ippos;(beamstart-TotLength)/1000;(TotLength/2)/1000;TotLength*2];
0314 posString=[];
0315 for i=1:length(positions)-1
0316     posString=[posString num2str(positions(i)*1000) ','];
0317 end
0318 posString=[posString num2str(positions(i)*1000)];
0319 
0320 twissH=['profile ' ...
0321     'zloop=' num2str(beamstart*0) ':' num2str(TotLength*3) ':1000 '...
0322     ...' z=' posString ...
0323     ' file=twissLER particle=' particle ...
0324     ' coordinates=centerline\n']; %   ' coordinates=reference\n'];
0325 
0326 %twissH=['profile zloop=' num2str(beamstart*0) ':' num2str(beamstart+TotLength*2.5) ':1000 file=twissHER particle=' particle ...
0327 %      ' coordinates=centerline\n']; %   ' coordinates=reference\n'];
0328 
0329 
0330 dump=['tubs DUMP innerRadius=0.01 outerRadius=300 length=10 kill=1 material=Cu color=0.2,0.2,0\n '...
0331     ];
0332 
0333 % some hidrogen for interactions
0334 ip=['material ebunch a=0.005 z=1 density=0.01\n'...
0335     'tubs ip innerRadius=0 outerRadius=300 length=1 kill=' num2str(killval) ' material=Vacuum color=1,.0,.1 \n '...
0336     ];
0337 plac{ord+2}=['# place ip parent=  z=' num2str(ippos,15) ' \n' ];
0338 plac{ord+3}=[...
0339     posbeam 'place BeamVis z=' num2str(beamstart,15) ' y=101 \n'
0340     ];
0341 plac{ord+4}=['# place DUMP parent=  z=' num2str(Dumppos,15) ' \n'];
0342 plac{ord+5}=[' cornerarc' ' z=0+' num2str(TotLength) '+0+' num2str(0,15) ...
0343     ' angle=' num2str(Xang/2*180/pi,15)  ' centerRadius='  num2str(1,15) '\n'...
0344     ];
0345 
0346 %%%%%%                             RF1         RF2                   beam
0347 [s,ind]=sort([0;Q{2}; S{2}; B{2}; RF1pos/1000; RF2pos/1000 ;637.2728766;beamstart/1000;Dumppos/1000;TotLength/1000]); % in m
0348 
0349 
0350 HER=[head defQ defS defB RF ip dump];
0351 HERPlace=[plac{ind}];
0352 
0353 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       LER          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0356 
0357 
0358 %%% LER (solenoids to implement)
0359 folder='SBLERV12DATA';
0360 P_0=4.179999969; % GeV/c
0361 particle='e-';
0362 
0363 Brho=3.3356409519815205*P_0; % electron
0364 
0365 quadL=-1; % swich off LER quad.
0366 
0367 
0368 killval=1;
0369 
0370 beamstart=575.4563172*1000-140; % first QDI center
0371 
0372 beamstart=0;
0373 
0374 %beamstart=575.4563172*1000-140*3; % first QDI center
0375 beamstart=(beamstart+TotLength);
0376 %%%%%  %%%% %%%% both beams heading in the same direction!
0377 elebeam=[ 'beam gaussian particle=' particle ' beamZ=' num2str(beamstart) ...
0378     ...' rotation=Y' num2str(1*Xang*180/pi) ... % if starting at IP.
0379     ...' rotation=Y' num2str(180) ... % if starting in straight.
0380     ' beamX=' num2str(-1*5.634489178e-05,15) ...
0381     ' beamXp=' num2str(+1*5.584726107e-08,15) ...
0382     '  nEvents=' num2str(Nevents) ' '... % nEvents=10
0383     ' sigmaX=' num2str(sX,15) ...
0384     ' sigmaY=' num2str(sY,15)...
0385     ' sigmaXp=' num2str(sXp,15)...
0386     ' sigmaYp=' num2str(sYp,15) ' '...
0387     '    meanMomentum=' num2str(P_0*1000,15) ' sigmaP=0'...
0388     ' meanT=0 sigmaT=0 \n'...
0389     'trackcolor reference=1,1,1 \n'...
0390     ];
0391 
0392 headL=[...
0393     'reference particle=' particle ' referenceMomentum=' num2str(P_0*1000,15) ' noEloss=1'...
0394     ' beamXp=0 beamYp=0 '...
0395     ' beamZ=' num2str(TotLength) '+0 '...
0396     ' rotation=Y' num2str(0) ' '... % elements in the opposite direction
0397     ' tuneMomentum=' num2str(P_0,15) '\n'...
0398     ];
0399 
0400 
0401 %headL=[...
0402 %];
0403 % x damping time= 0.26700741E-01 => 15% =~ 40ms = 40*10^6 ns
0404 
0405 % get pos K and L for Quad, sext and bend
0406 Qfile=fopen([folder '/quad']);
0407 Sfile=fopen([folder '/sext']);
0408 Bfile=fopen([folder '/bend']);
0409 SOLfile=fopen([folder '/sol']);
0410 % get pos K and L for Quad, sext and bend -- sbend and not rbend.
0411 %Qfile=fopen('DIAMDATA/quad');
0412 %Sfile=fopen('DIAMDATA/sext');
0413 %Bfile=fopen('DIAMDATA/bend');
0414 
0415 % name s K L
0416 Q=textscan(Qfile,'%s %f %f %f','Headerlines',47);
0417 B=textscan(Bfile,'%s %f %f %f','Headerlines',47);
0418 S=textscan(Sfile,'%s %f %f %f','Headerlines',47);
0419 SOL=textscan(SOLfile,'%s %f %f %f','Headerlines',47);
0420 
0421 defQ= ['virtualdetector Det radius=10.0 length=' num2str(2*StepLength)...
0422     ' color=0,.2,.0 referenceParticle=1 format=rootExtended\n'];
0423 defS='';
0424 defSOL='';
0425 defB='';
0426 ord=1;
0427 placL={};
0428 defined={};
0429 d=1;
0430 % tilt reference by crossing angle.
0431 placL{ord}=[' cornerarc' ' z=0+' num2str(TotLength) '+0+' num2str(0,15) ...
0432     ' angle=' num2str(Xang/2*180/pi,15)  ' centerRadius='  num2str(1,15) '\n'...
0433     ];
0434 ord=ord+1;
0435 
0436 for i=1:length(Q{1})
0437     name=(Q{1}(i));
0438     name{1}(1)=' ';
0439     name{1}(end)='p';
0440     Q{1}(i)=name;
0441 end
0442 for i=1:length(S{1})
0443     name=S{1}(i);
0444     name{1}(1)=' ';
0445     name{1}(end)='p';
0446     S{1}(i)=name;
0447 end
0448 for i=1:length(B{1})
0449     name=B{1}(i);
0450     name{1}(1)='';
0451     name{1}(end)='p';
0452     B{1}(i)=name;
0453 end
0454 for i=1:length(SOL{1})
0455     name=SOL{1}(i);
0456     name{1}(1)='';
0457     name{1}(end)='p';
0458     SOL{1}(i)=name;
0459 end
0460 
0461 % small FF quad
0462 ffQuad{1}=char(Q{1}(1));%{'QD0AL'};
0463 ffQuad{2}=char(Q{1}(2));%{'QD0BL'};
0464 ffQuad{3}=char(Q{1}(3));%{'QD0AR'};
0465 ffQuad{4}=char(Q{1}(4));
0466 ffQuad{5}=char(Q{1}(end));%{'QD0BR'};
0467 ffQuad{6}=char(Q{1}(end-1));%{'QDPMR'};
0468 ffQuad{7}=char(Q{1}(end-2));%{'QDPML'};
0469 ffQuad{8}=char(Q{1}(end-3));
0470 
0471 
0472 for i=1:length(Q{1})
0473     %genericquad QDpi fieldLength=280 apertureRadius=101.5 ironRadius=381 \
0474     %    ironLength=250 ironColor=0.6,.0,.6 kill=' num2str(killval) ' gradient=$KQDpi
0475     %
0476     if sum(strcmp(char(Q{1}(i)),defined))==0;
0477         inrad=60;
0478         outrad=70;
0479         
0480         if sum(strcmp(char(Q{1}(i)),ffQuad))~=0
0481             inrad=12;
0482             outrad=15;
0483             tempSL=StepLength;
0484             StepLength=0.01;
0485             if strcmp(char(Q{1}(i)),ffQuad{1}) || strcmp(char(Q{1}(i)),ffQuad{5}) %% P.M.
0486                 inrad=8;
0487                 outrad=10;
0488                 
0489             end
0490         end
0491         
0492         defQ=[ defQ ...
0493             ' genericquad ' Q{1}(i) ' '...
0494             ' fieldLength=' num2str(Q{4}(i)*1000,15) ... % mm
0495             ' apertureRadius=' num2str(inrad,15) ' '...
0496             ' ironRadius=' num2str(outrad,15) ' '...
0497             ' ironLength=' num2str(Q{4}(i)*1000,15) ' '...
0498             ' ironColor=0/255,255/255,191/255 kill=' num2str(killval) ' '...
0499             ' fringe=0'...
0500             ' gradient=' num2str(quadL*Q{3}(i)*Brho/Q{4}(i),15) ' '... K
0501             ' maxStep=' num2str(StepLength/2)  ...
0502             '\n']; %#ok<*AGROW>
0503         
0504         % add to already defined
0505         defined{d}=char(Q{1}(i));
0506         d=d+1;
0507     end
0508     
0509     StepLength=tempSL;
0510      
0511     %defQ=[defQ 'virtualdetector Det radius=400.0 color=0,.1,.1'];
0512     
0513     if (ord/2-floor(ord/2))==(det-1) % set 0 to place a detector in every quadrupole.
0514         placL{ord}=['place ' Q{1}(i) ' parent=  z=0+' num2str(TotLength) '+0+' num2str(Q{2}(i)*1000,15)...
0515             ' y=0'...
0516             ' coordinates=centerline'...
0517             '\n'...
0518             'place ' 'Det' ' parent=  z=0+' num2str(TotLength) '+0+' num2str((Q{2}(i))*1000)...
0519             ' rename=DETp# '...
0520             ' y=0'...
0521             ' coordinates=centerline '...
0522             '\n']; %#ok<*SAGROW>
0523     else
0524         placL{ord}=['place ' Q{1}(i) ' parent=  z=0+' num2str(TotLength) '+0+' num2str(Q{2}(i)*1000,15)...
0525             ' y=0'...
0526             ' coordinates=centerline'...
0527             '\n']; %#ok<*SAGROW>
0528     end
0529     %      if sum(strcmp(Q{1}(i),{'"QD0AL"';'"QD0BL"';'"QD0AR"';'"QD0BR"'}))>0
0530     %
0531     %      plac{ord}=['place ' Q{1}(i) ' parent=  z=' num2str(Q{2}(i)*1000,15)...
0532     %      ' coordinates=centerline'...
0533     %      ' rotation=Y' num2str(-0.066*180/pi,15) ' '...
0534     %      '\n'];
0535     %     end
0536     ord=ord+1;
0537 end
0538 
0539 for i=1:length(S{1})
0540     %multipole SF1 fieldLength=150 apertureRadius=101.5 ironRadius=381 \
0541     %    ironLength=140 ironColor=0,.6,.6 kill=' num2str(killval) ' sextupole=$KSF1
0542     
0543     if sum(strcmp(char(S{1}(i)),defined))==0;
0544         inrad=80;
0545         outrad=120;
0546         
0547         defS=[ defS ...
0548             ' multipole ' S{1}(i) ' '...
0549             ' fieldLength=' num2str(S{4}(i)*1000,15) ... % mm
0550             ' apertureRadius=' num2str(inrad,15) ' '...
0551             ' ironRadius=' num2str(outrad,15) ' '...
0552             ' ironLength=' num2str(S{4}(i)*1000,15) ' '...
0553             ' ironColor=173/255,255/255,47/255 kill=' num2str(killval) ' '...
0554             ' fringe=0'...
0555             ' maxStep=' num2str(StepLength/2)  ...
0556             ' sextupole=' num2str(sxt*S{3}(i)*Brho/S{4}(i),15) ' '... K2
0557             '\n'];
0558         % add to already defined
0559         defined{d}=char(S{1}(i));
0560         d=d+1;
0561     end
0562     
0563     
0564     placL{ord}=['place ' S{1}(i) ' parent=  z=0+' num2str(TotLength) '+0+' num2str(S{2}(i)*1000,15)...
0565         ' y=0'...
0566         ' coordinates=centerline '...
0567         '\n']; %#ok<*SAGROW>
0568     ord=ord+1;
0569     
0570 end
0571 
0572 for i=1:length(SOL{1})
0573     %multipole SF1 fieldLength=150 apertureRadius=101.5 ironRadius=381 \
0574     %    ironLength=140 ironColor=0,.6,.6 kill=' num2str(killval) ' sextupole=$KSF1
0575     
0576     if sum(strcmp(char(SOL{1}(i)),defined))==0;
0577         inrad=80;
0578         outrad=120;
0579         
0580         defSOL=[ defSOL ...
0581             ' \ncoil C' SOL{1}(i) ...
0582             ' innerRadius=' num2str(inrad) ...
0583             ' outerRadius=' num2str(outrad) ...
0584             ' length=' num2str(SOL{4}(i)) ...
0585             ' material=Cu' ...
0586             ' maxR=' num2str(outrad) '\n\n'...
0587             ' solenoid  ' SOL{1}(i) ' coil=C' SOL{1}(i) ' '...
0588             ' current=' num2str(0) ... %%% fake!!!!! not ok!!!
0589             ' color=0,1,1 kill=' num2str(killval) ...
0590             '\n\n'];
0591         % add to already defined
0592         defined{d}=char(SOL{1}(i));
0593         d=d+1;
0594     end
0595     
0596     
0597     placL{ord}=['place ' SOL{1}(i) ' parent=  z=0+' num2str(TotLength) '+0+' num2str(SOL{2}(i)*1000,15)...
0598         ' y=0'...
0599         ' coordinates=centerline '...
0600         '\n']; %#ok<*SAGROW>
0601     ord=ord+1;
0602     
0603 end
0604 
0605 for i=1:length(B{1})
0606     
0607     %genericbend BHer fieldWidth=1660 fieldHeight=200 fieldLength=5400 \
0608     %    ironColor=1,0,0 ironWidth=1828 ironHeight=1320 ironLength=5000 \
0609     % kill=' num2str(killval) ' By=0.2778157677856496
0610     if sum(strcmp(char(B{1}(i)),defined))==0 % if not defined, define.
0611         
0612         Fw=240;
0613         Fh=70;
0614         Iw=250;
0615         Ih=100;
0616         
0617         defB=[ defB ...
0618             '\n param Bp' num2str(i) B{1}(i) '=' num2str(Brho*B{3}(i)/B{4}(i),15) '\n'...
0619             ];
0620         if B{2}(i)<400 || B{2}(i)>800
0621             defB=[ defB ...
0622                 %  ' tune Bp' num2str(i) B{1}(i) ' z0=' num2str(B{2}(i)*1000-B{4}(i)*500) ' z1=' num2str(B{2}(i)*1000+B{4}(i)*500) ' '...
0623                 %  ' initial=' num2str(Brho*B{3}(i)/B{4}(i),15) ' step=1e-13'...
0624                 %  ' expr=(Px1/Pz1) tolerance=0.0000001 \n'...
0625                 ];
0626         end
0627         defB=[ defB ...
0628             ' genericbend ' B{1}(i) ' '...
0629             ' fieldLength=' num2str(B{4}(i)*1000,15) ... % mm
0630             ' fieldWidth=' num2str(Fw,15) ' '...
0631             ' fieldHeight=' num2str(Fh,15) ' '...
0632             ' ironLength=' num2str(B{4}(i)*1000,15) ... % mm
0633             ' ironWidth=' num2str(Iw,15) ' '...
0634             ' ironHeight=' num2str(Ih,15) ' '...
0635             ' ironLength=' num2str(B{4}(i)*1000,15) ' '...
0636             ' fringe=0'...
0637             ' maxStep=' num2str(StepLength)  ...
0638             ...' By=' num2str(Brho*B{3}(i)/B{4}(i),15)...' fieldColor=0.5,.1,.0 '...
0639             ' By=$Bp' num2str(i) B{1}(i) ...
0640             ' ironColor=1,.0,.0 kill=' num2str(killval) ' '...
0641             '\n'...
0642             ];
0643         % add to already defined
0644         defined{d}=char(B{1}(i));
0645         d=d+1;
0646         
0647     end
0648     
0649     
0650     %'# cornerarc' ' parent=  z=' num2str((B{2}(i))*1000,15) ...    ' angle=' num2str(B{3}(i)*180/pi) ' centerRadius='  num2str((B{4}(i)/B{3}(i))*1000,15) '\n'...'# cornerarc' ' parent=  z=' num2str((B{2}(i)-B{4}(i)/2)*1000,15) ...    ' angle=' num2str(B{3}(i)*180/pi/2,15) ' centerRadius='  num2str((B{4}(i)/B{3}(i))*1000,15) '\n'...'#corner' ' parent=  z=' num2str((B{2}(i)-B{4}(i)/2)*1000,15) ...' rotation=Y' num2str(B{3}(i)*180/pi/2,15) ' \n'...
0651     placL{ord}=[
0652         'place ' B{1}(i) ' parent=  z=0+' num2str(TotLength) '+0+' num2str(B{2}(i)*1000,15) ...
0653         ' x=' num2str(B{4}(i)/2*tan(B{3}(i)/2)*1000) ''...
0654         ' coordinates=centerline'...' \n'...
0655         ' rotation=Y' num2str(B{3}(i)*180/pi/2,15) ...' rename=' B{1}(i) ' \n'...
0656         '\n' ...
0657         ' cornerarc' ' z=0+' num2str(TotLength) '+0+' num2str((B{2}(i)-B{4}(i)/2)*1000,15) ...
0658         ' angle=' num2str(B{3}(i)*180/pi,15)  ' centerRadius='  num2str((B{4}(i)/B{3}(i))*1000,15) '\n'...'#corner' ' parent=  z=' num2str((B{2}(i)+B{4}(i)/2)*1000,15) ... ' rotation=Y' num2str(B{3}(i)*180/pi/2,15) ' \n'...
0659         '\n']; %#ok<*SAGROW>
0660     ord=ord+1;
0661     
0662 end
0663 
0664 RF_on=0;
0665 
0666 % plac RFC
0667 %length[m]        voltage[MV]                lag          freq[MHz]             harmon
0668 %0.501777                0.6           0.424389        476.0054439               1998
0669 RFL=['pillbox RFC maxGradient=' num2str(RF_on*0.6/0.501777,15) ...
0670     ' color=0.5,0.5,0 frequency=0.4760054439 innerLength=501.777 '...
0671     ' phaseAcc=0.424389 wallThick=100 innerRadius=1000\n'];
0672 RF1pos=TotLength+577.1465514*1000;
0673 RF2pos=TotLength+579.1935514*1000;
0674 
0675 placL{ord}=['  # place RFC z=' num2str(RF1pos) '\n'];
0676 placL{ord+1}=['  # place RFC z=' num2str(RF2pos) '\n'];
0677 
0678 %rfc, at = 11.2402595;
0679 %rfc, at = 13.2872595;
0680 positions=[0;Q{2}; S{2}; SOL{2}; B{2}; (RF1pos-TotLength)/1000; (RF2pos-TotLength)/1000 ;ippos;(beamstart-TotLength)/1000;(TotLength/2)/1000;TotLength*2];
0681 posString=[];
0682 for i=1:length(positions)-1
0683     posString=[posString num2str(positions(i)*1000+TotLength) ','];
0684 end
0685 posString=[posString num2str(positions(i)*1000+TotLength)];
0686 
0687 twissL=['profile ' ...
0688     'zloop=' num2str(beamstart*0) ':' num2str(TotLength*3) ':1000 '...
0689     ...' z=' posString ...
0690     ' file=twissLER particle=' particle ...
0691     ' coordinates=centerline\n']; %   ' coordinates=reference\n'];
0692 
0693 % keep only part of sequence.
0694 % ind=ind(end/4:end/4*2);
0695 
0696 dumpL=['tubs DUMP innerRadius=0 outerRadius=300 length=10 kill=1 material=Cu color=0.2,0.2,0\n '...
0697     ];
0698 
0699 % some hidrogen for interactions
0700 ipL=['material ebunch a=0.005 z=1 density=0.01\n'...
0701     'tubs ip innerRadius=0 outerRadius=300 length=1 kill=' num2str(killval) ' material=Vacuum color=1,.0,.1 \n '...
0702     ];
0703 ippos=0;
0704 Dumppos=beamstart+(-11-280);
0705 placL{ord+2}=['# place ip parent=  z=' num2str(ippos,15) ' \n' ];
0706 placL{ord+3}=[elebeam 'place BeamVis z=' num2str(beamstart,15) ' y=-101 \n'];
0707 placL{ord+4}=['# place DUMP parent=  z=' num2str(Dumppos,15) ' \n'];
0708 
0709 placL{ord+5}=[' cornerarc' ' z=0+' num2str(TotLength*2) '+0+' num2str(0,15) ...
0710     ' angle=' num2str(-Xang/2*180/pi,15)  ' centerRadius='  num2str(1,15) '\n'...
0711     ];
0712 
0713 
0714 % first 0 is for cornerarc crossing angle rotation.
0715 [s,ind]=sort([0;Q{2}; S{2}; SOL{2}; B{2}; (RF1pos-TotLength)/1000; (RF2pos-TotLength)/1000 ;ippos;(beamstart-TotLength)/1000;(TotLength/2)/1000;TotLength*2]); % in m
0716 
0717 
0718 
0719 outtext=char(cell2mat([...
0720     bigmondo ...
0721     ... definitions
0722     HER ...  HER
0723     defQ defS defB defSOL ...  LER
0724     RF ip...
0725     HERPlace ...  Place HER magnets
0726     twissH ...
0727     headL  dumpL...
0728     placL{ind} twissL...  Place LER magnets
0729     ' particlecolor e-=0,0,1 e+=1,0,0 \n'...  electron is blue and positron is red.
0730     ...'g4ui when=4 ''/vis/open OGL''\n'...
0731     'g4ui when=4 ''/vis/viewer/set/viewpointThetaPhi 45 45''\n'...
0732     'g4ui when=4 ''/vis/viewer/set/style wireframe''\n'...
0733     'g4ui when=4 ''/run/beamOn ' num2str(Nevents*2) ' ''\n'...
0734     ... 'g4ui when=4 ''vis/ogl/printEPS''\n'...
0735     ]));
0736 % lattice 'g4ui when=4 ''/vis/viewer/set/style wireframe''\n'...    'g4ui ''/run/beamOn 20''\n'...
0737 out=fopen('G4blseqconv','w+');
0738 fprintf(out,outtext);
0739 
0740 fclose('all');
0741 clear defQ
0742 clear defS
0743 clear defB
0744 clear plac ind ord
0745 clear defined
0746 disp('--- DONE ---')