Embedded Files
2019 Logbook
Week 1 (Jan 31-Feb 7):
Overview of lab:
- This lab uses high intensity laser pulses to study their various interactions/effects with matter.
- Experiments involve: nonlinear optics, atomic physics, condensed matter physics, plasma physics, and quantum mechanics.
- Research topics include: spatio-temporal optical vortices (STOVs), electron acceleration in plasma, plasma waveguides, and nonlinear optics in the strong field regime.
Started reading/learning how lasers work:
-“Laser” stands for Light Amplification by Stimulated Emission of Radiation
-A laser works by exciting an atom so that it’s electrons jump to a higher energy orbit. When the electron returns to it’s ground state, it emits a photon, and these photons are what makes up a laser beam.
-Because the laser is made up from the same amount of energy at the same time, laser light only has one wavelength and the waves are in sync, making the light stay focused for long distances.
-Additionally, the photons are all emitted simultaneously, so the waves are in sync, making the light stay focused for long distances.
High Intensity Lasers
-Its very difficult to get a continuous laser wave, even with a constant power source, to reach very high average intensities.
-However, laser pulses can get to much larger intensities but emitting large bursts of energy for a merely microseconds.
-These pulses can achieve power of several TW (10^12 watts)
-They do this with CPA, which stands for Chirped Pulse Amplification
Chirped Pulse Amplification
- CPA is a technique used to transform peak powers from laser pulses from several GW to several TW.
- CPA does this by the fact that laser pulses are broadband, so short pulses have to be made out of larger wave frequencies.
-These lasers are amplified using detraction so that the intensity of the laser doesn’t damage internal laser components.
-These components are then put back together with a prism so that the laser has full intensity (depicted on next slide)
My Project
-I’m starting by just learning about the lab and how to use the machine shop.
-Then I’ll be assigned to building parts for the laser
-I’m not entirely sure what I’ll do after that, but hopefully I’ll eventually to do what previous undergraduate students have done:
-designed and built experiments and coded theoretical models in MATLAB for plasma behavior under Maxwell’s equations (among other things).
Week 2 (Feb 7 - Feb 14):
Started meeting people:
- This week I didn’t start formally like I had hoped.
- Continued to learn about lasers (checked out books from library)
- Met Dr. Milchberg’s coworkers and students in lab
- Emailed him multiple times, and although didn’t get starting date until next Monday, got invitation to dinner with Patricia Falcone.
-Dr. Falcone is the deputy director for Science and Technology at Lawrence Livermore National Laboratory (LLNL).
-She joined LLNL in 2015 after working at the White House Office of Science and Technology Policy (OSTP)
Continued Research on Lasers:
- Checked out two books from STEM Library on lasers and am half way through one of them.
- Set up meeting for Feb 18th to formally start.
Week 3 (Feb 14 - Feb 21):
Continued reading up on lasers.
Feb 18th: met with Dr. Milchberg, received forms that I need to fill out before properly starting.
- Met with grad students that got me started doing drawings on paper and then on AutoCad of machine parts.
- Drew a "dog".... unfortunately its not the cute kind (see picture)
- These are what they use to attach lasers and other machinery to the tables, so I start by just learning out to draw and build stuff.
- Watched videos about how to do technical drawings
- Watched videos about how to use AutoCAD
- Downloaded & Installed AutoCAD
AutoCAD Functions:
- to make a rectangle, start writing rectangle, and it will suggest the shape, click on the auction "rectangle" and then you have options on how to make the rectangle the correct size.
- I chose to make the rectangle's but specifying the base point coordinate, and then the coordinate of the diagonal corner of the rectangle using the base point as an origin. I think this works quite well.
- to make a circle, it's also quite simple. You start typing circle, click on the circle function, and specify the coordinate of the center and the radius.
- making the slot was more complicated. You need to make two circles at the edges for the round edges, and then make lines (using the line function) to connect the top and bottom of both circles. (instead of using coordinates for the line, use the mouse to click on where you want them. It will automatically make them tangent to the circle and connect in the center. Finally, use the trim tool to select only the parts of the circles that you want, and the rest will automatically delete.
- to make measurements, go to the "dimensions" tab and then click on what you want to measure. It's fairly intuitive.
Drawing of Dog
AutoCAD Drawing of Dog
Actual Dog
Week 4: (Feb 21 - Feb 28)
- 2/22/19: (2:00-5:00pm)
Worked in lab for the afternoon.
- got an extensive tour, obtained my own safety equipment, learned about potential projects, and got further training in machine shop.
Tour:
- Dr. Milchberg has many labs.
- The main one is called Big Lab and has a very large Sapphire laser set up.
- This laser is useful because it can change to a variety of wavelengths and can be very powerful.
- It's currently used to work on plasma waveguides, which is where they use high density gas, turn it into plasma with the laser (i.e. stripping the gas's electrons).
- They use that to make high density plasma on the edges of the laser, and low density plasma on the center (using low density gas).
- It turns out that the works similarly to a fiberoptic, maintaining the lasers focus for longer amounts of time.
- It's also useful because it high intensity lasers can't use glass fiber optics because they have too much energy.
- The New Lab:
- This is where I work, and if you go into the Big lab, right across from the entrance is a door.
- The door has a bright, lit up sign above that says "LASER ON".
- Technically they all have that sign, but this one is noticeable because one of the lasers in there is always on.
- That one is sapphire similar to the one in Big Lab, and I'm not sure what they do on it.
- The laser I'll be initially working with has very short pulses close to the inferred scale. (called an IR laser).
- I learned more about this on Tuesday (2/26/19) so I'll come back to this laser later.
Safety Procedures:
- Glasses: I learned more about this on Monday (2/25/19) but lasers are very dangerous to both skin and eyes.
- additionally, different lasers with different wavelengths require different types of glasses.
- If you ever go into a laser lab, be mindful to wear glasses approved for the wavelength and intensity of the laser.
- Also, be extra careful when aligning lasers.
- Often, people have a tendency to take off glasses to see laser better (this is the most common cause of injury).
- Skin: For skin it's a little more simple. Basically don't but your hand in the way of a laser because it will hurt, a lot.
- Electricity: there is a lot of power going into these lasers, so the getting electrocuted is a big hazard. Be careful around power sources.
- Safety courses: I'm required to take safety courses (laser safety, chemical safety, x-ray radiation machine safety). So I signed up for those.
- The sign up can be confusing. Go here and select the course you have to take.
Machine Shop:
- Learned about making metal pieces by watching some grad students build something for the first time. I still have to take a more extensive training.
- 2/25/19: (10:30-12:30pm)
- Took Laser Safety training course and Chemical safety course.
Chemical Safety course:
- This course is online and I did over the weekend, but I'm putting it here for conciseness.
- It was fairly straight forward, I just had to learn about different types of dangerous chemicals and what their labels mean, etc.
Laser Safety course:
- This course was an in class course that I had to attend.
- The director of laser safety walked us through an extensive power point and explained what lasers are, how they work, and their dangers.
- Much of the safety stuff I explained above, but additionally it's important to minimize anywhere the laser can go (where it might interact with a human)
- For example, you want to have the end of a laser going towards a wall, or perhaps you want curtains going around the laser table.
- The most valuable information here (in addition to scaring me into being careful around lasers) was the information about glasses.
- You might think that such a small thing should be trivial, but really it's quite complex.
- Often you have a range for the potential wavelength of a laser, and there's a range of wavelengths that glasses protect you from.
- However, the range on the glasses does not always apply for the same intensity.
- The laser safety officer recommends calling the glasses company to make sure they are correct for your experiment.
- 2/26/19: (11:00-1:30pm)
- Got assigned a couple projects and readings
Projects:
1) Building and Autocorrelator
- This is by far the one I'm most excited about, and it's for the IR laser that I mentioned earlier
- So because IR lasers have such short pulses, and their wavelengths are so short, nothing else can really measure them.
- This is because to do that you would need something shorter, which doesn't exist.
- However, when doing an experiment you still need to know how long the pulse is, and you do that by using the laser beam to measure itself.
- An autocorrelator is a type of optics used to make a beam measure itself.
- They do this by splitting the beam, having them go on parallel paths, and using a non-linear crystal to have them converge at a focal point
- If the pulses are in sync, the detector and the focal point will get a large intensity reading, otherwise it will get small intensity readings.
- I'm inserting a diagram below, and I will further explain the diagram next week, as I don't understand it fully yet.
2) Building boxes/covers
- This is less exciting, but I think it will really improve my skills overall, and I'm excited to be comment at machining.
- These boxes will be used to cover the lasers and prevent dust from getting on the optics, and prevent the lasers from shining off the table into someone's eyes.
3) Organizing Optical Equipement:
- I'm trying to teach myself optics, and it's quite difficult, so organizing optical equipment might help.
- This was I'll be both useful to the lab, and I'll become familiar with the equipment (hopefully giving me a better understanding of optics in general)
How Laser Pulses Work:
- Additionally, I've been trying to teach myself as much as possible in my free time. I've reviewed the following (and probably more than I didn't record):
- Fourier Transforms
- Waves & Optics (with many subfields)
- Autocorrelators
- Optical equipment
- Intense laser pulse mechanics
Autocorrelator
Graph of Laser Glasses Protection for Wavelengths
Week 5: (Feb 28 - March 7)
2/28/19:
finally finished machine shop training and finished my dog
- training included learning to use mills, lathes, band saws, drills, and tapping aluminum.
3/1/19:
began designing dust covers and ordering materials
- had to search McMaster Carr's website for everything I needed, measure dimensions of boxes to make sure I had enough material, and put together order.
- also had to email whole group to see if anyone else had to order anything.
- put together order of 10 items, including fastenings, valves and piping for the vacuum/gas chambers, conducting tape, etc.
Weekend:
-finished double checking order list (sounds simple but there was a lot of money involved and I couldn't risk ordering the wrong things so took a very long time)
3/4/19:
X-ray machine radiation training and formally ordered parts
- I asked my last couple questions and sent my order request to Dr. Milchberg who approved it and passed it along. The parts arrived on Tuesday.
- I also went to x-ray machine radiation training to learn about how much radiation exposure is safe and how to minimize it.
3/5/19:
Picked up packages, put parts where they belonged, and read previous publications from group.
- in the labs, everything needs to be meticulously organized, so when the parts arrived I had to learn how everything was organized and put the parts that I ordered away.
- Read papers about 4 micron lasers (Mid IR lasers, which is what I'm working on).
- Read papers about autocorrelators (which is what I'm going to build, after I finish dust covers)
X-Ray Machine Radiation Safety
- Key guidelines: minimize the time of your exposure, increase distance from source, and use shielding to protect yourself from source.
- Wear detector on your body to keep track of how much exposure you have received
- There are limits on how much you can safely have, so using detectors keeps you within that limit
- Detectors:
- different detectors work different ways, but they're important to use in case of more radiation leaking than expected. Also, it's safer to have them to test how well different types of shielding, distance and time will effect your exposure.
- Some detectors are designed for photons, others less so.
- Hard to get very accurate detectors in terms of the amount of radiation. Normally just tell you if it's there or not.
AutoCorrelator:
Model I'm building is two photons, which means that it doesn't use the non-linear crystal the same way as described above.
Instead, we're selecting a material where you need the energy of two photons to boost an electron to a higher energy state.
This material has never been used in this way before, so this experiment will get us a paper, and it's a fairly simple experiment.
How does it work?
- like described above, if the pulses are in sync, you'll get two large amounts of energy on the material at once, and hopefully, those pulses will excite an electron up to it's higher energy stake (which smaller amounts of energy cannot do).
- We will detect this with a voltage attacked to the material.
A picture of the set up is below:
Week 6: (March 7 - March 14)
Edited research paper & built box covers (or tried to)
3/7/19
- One of the grad students who supervises me asked me to read/edit his paper on electron avalanches to make sure the paper made sense to someone who didn't do the experiment.
- The paper took about two hours to read, and I made comments whenever I got lost somewhere or found wording confusing.
- The paper itself was really cool. It turns out these mid IR lasers are really good at accelerating electrons that then accelerate each other and make what they called an electron avalanche.
- I would post some of the graphs below, but the paper hasn't been finally published yet so I'm not going to hold off on that.
3/8/19
- Today I did a final drawing for one of the box lids, and started working on it in the machine shop. I was a little intimidated, so a grad student helped me get started.
- The drawing I did for the lid was just a hand drawing, which I had thought would be good enough given that I couldn't get the AutoCAD to work that day.... I was wrong.
- Potential grad students who are applying to UMD also came through and visited.
- We gave them a tour and answered their questions, and it was interesting to see that side of the process—I've never thought about and of the specifics of grad school, so that was fun.
3/11/19
- So today I tried to finish the box lid, which I did (sorta).
- I used a band saw, a mill, and hack saw, and then sand paper to make the edges smooth.
- The final product looked really good, so I brought it over to the box... and it didn't fit
- There were a couple protruding edges that I hadn't considered when building it, so I need to go back and make some of the cut outs bigger.
- I guess I learned to make my measurements more exact, or to leave more room for error.
That's it for this week because I was sick Tuesday
Below I'm inserting images of the tools I used, and what I used them for:
Above is a Hack Saw:
- I used this to cut out a full square from my box lid. It's useful because it can turn, where as other saws cannot because they are too thick.
Above is a Band Saw:
- These are basically electric saws were the blade just goes down. I used this to cut my plastic into a general shape.
Above is a Mill:
- I used this to drill the wholes in the plastic for the handles.
Above is a detailed drawing of the handles that I'm using
Above is the manufactures picture of the plastic I used: (Delrin® Acetal Resin Sheets)
Week 7: (March 14 - March 27)
This week I didn't do as much because of spring break.
I took some time to develop a better understanding of Fourier Transforms, did some readings on high power lasers, and started designing the gas chamber for the next experiment.
I did a lot of learning about the next experiment. The idea is to analyze how different gasses ionize, so I'm making a chamber that has tapped wholes to be mounted on, tapped holes for the gas tubes to connect to, and glass windows (vacuums sealed) for the laser to come into the chamber through.
I also learned a lot about how to aline the lasers with "irises," "posts," etc. I'll put some pictures below, but basically, between each optic, you need to make sure the beam is at the right height and the center of the optic. (If it's off, then the beam can be changed because the lenses/mirrors we use are mostly spherical instead of parabolic. (Spherical lenses/mirrors are much cheaper and easier to get than parabolic ones, and their a fairly good approximation).
Despite all the information I learned above, my highlight of the week was touring a variety of labs with the potential grad students on Friday. I got to see Lathrop's Giant sodium ball that they use for simulating the Earth's magnetic field, which was really really cool. It was 5 meters in diameter and surrounded by coils for measuring inducted fields. Apparently not too long ago, a guy crawled into their lab through the vents and stole the copper off of one of the safety machines, which is scary to think about but also really funny. Scary because if he had stollen it off the sodium ball, the sodium could have exploded and the whole building would be gone.
When the run the sodium ball, they spin it as a speed proportional to that of the earth, and the whole room shakes. Soon, they're going to drain the sphere and put bevels in the interior sphere, which means they're going to have a shipping container of liquid sodium just sitting in their lab.
The measures they take in case of malfunction were interesting as well.
I also got to tour other laser labs (which were similar to mine, but slightly different in different ways). It was interesting how many of the other laser lab referenced my grad student's research when explaining stuff to us, and I really enjoyed actually understanding more of what they said this time as they explained the lasers to me and the other perspective grad students.
Optic equipment:
Below is an iris, which you close and hold a card on the other side to see if it disfigures the beam asymmetrically. If it does, the beam needs to be adjusted.
Below is a depiction of a spherical vs. parabolic mirror. It's interesting that one of the labs I tourned only seamed to use parabolic mirror, where I think mine pretty much always used spherical.
I also learned how to clean equipment (there are very specific techniques).
Week 8: (March 27 - April 3)
This week I finished some of the dust covers, and started working on a matlab project.
The idea of the matlab project is to analyze data from the experiments.
So when we use lasers to ionize gases, electrons leave the atoms.
Where the atoms return to the gas particle, the particle releases light (or phosphorescence).
Because our laser is inferred, we can use cameras to image the focal point of the laser, and see the individual phosphorescence for each ionized gas particle.
This means that if we can count these phosphorescence, we can count the exact number of electron avalanches that occur, so this is what my matlab program is supposed to do: count the number and brightness of each phosphorescence.
Here is an example of the type of image we would get from data:
If you look very closely to the center, there are small white dots in a horizontal line in the center.
Here's a zoomed in image of the dots:
So, when approaching the matlab code, I came up with a couple different strategies.
The first one is that matlab has a function that can detect circles/patches, but they have to be more than 3 pixels in radius. I used this in the code written below, and it worked really well for photos where dusk in the laser caused a giant bright spot, but it didn't pick up any of these smaller spots.
tdata = 'test.tif'; %selects data
t = Tiff(tdata, 'r'); %unpacks Tiff format
info = imfinfo(tdata);
numImages = numel(info); %counts number of images in tiff
tnum = 17; % choose number of image in tiff to analyze
image = imread(tdata,tnum); %reads desired image
%This works to find BIG bright circles (good for dust)
%Does not see any of the smaller ones, probably because they're too small
%to look like cirlces
Rmin = 2; %min radius
Rmax = 25; %max radius
[centersBright, radiiBright] = imfindcircles(image,[Rmin Rmax],'ObjectPolarity','bright') %uses matlab fuction to find bright spots
viscircles(centersBright, radiiBright,'Color','b'); %returns image with cirles
Output:
centersBright =
647.5282 99.7289
radiiBright =
8.6344
This only returns anything because I chose an image with a giant dust spot, see image below.
Because this didn't work, I tried using the fact that the number value in the image matrix is much higher if the color is lighter. Since this is the case, I could use relative peaks to locate and evaluate each bright spot.
I tried writing one myself, but it didn't work very well, so I went to matlab file exchange and found a function called FastPeakFinder that works really well. It even takes other variables besides the data to try to reduce noise. I used this function to right the following:
tdata = 'test.tif'; %selects data
t = Tiff(tdata, 'r'); %unpacks Tiff format
info = imfinfo(tdata);
numImages = numel(info); %counts number of images in tiff
tnum = 1; %choose number of image in tiff to analyze
image = imread(tdata,tnum); %reads desired image
%Finds peaks of Matrix
thres = 1000; %minimum peak value
p=FastPeakFind(image, thres); %function that finds peaks of 2D matrix
imagesc(image); hold on
plot(p(1:2:end),p(2:2:end),'r+') %plots peaks
And it worked really well. This was my output:
And when you zoom in on the bright spots:
You get both the center of the bright spot, and an increase in intensity.
I'll probably keep working on this for next week, and try to get it good for all sorts of brightness (maybe combining the two types of detection methods) and reduces noise as much as possible. That way it'll be all set for when it comes time to take new data.
Week 9: (April 3 - April 10)
This week I made some very large strides in the code I was working on last week, and made some new code that analyzed the data in a different way.
So the code from last week would detect peaks, even if they weren't on the line. My new code runs the find peaks program twice, uses boundaries and filters, and as a result only finds the desired peaks and gets rid of almost all the noise.
Here's the new code:
tdata = 'RA22_6_copy.tif'; %selects data
t = Tiff(tdata, 'r'); %unpacks Tiff format
info = imfinfo(tdata);
numImages = numel(info); %counts number of images in tiff
tnum = 521; %choose number of image in tiff to analyze
image = imread(tdata,tnum); %reads desired image
%Finds peaks of Matrix
filt = (fspecial('gaussian', 7,1)); %filter???????
avgPix = mean(image(:));
thres = 3 * avgPix; %minimum peak value
edg = 80; %skips edge pixels
p=FastPeakFind(image, thres, filt, edg); %function that finds peaks of 2D matrix
%use found peaks to reduce noise
x = p(1:2:end); % x-coordinates of peaks
y = p(2:2:end); % y-coordinates of peaks
yLine = median(y); % y coordinate of most of the data
bounds = round(yLine - std(y)); %new bounds to get rid of noise
p=FastPeakFind(image, thres, filt, bounds, 1); %function that finds peaks of 2D matrix
x = p(1:2:end); % new x-coordinates of peaks (with less noise)
y = p(2:2:end); % new y-coordinates of peaks (with less noise)
%plot peaks on graph
imagesc(image); hold on
plot(x, y, 'r+');
Output:
This one makes a graph of all the rows of pixel values with peaks, and then an average of those rows:
tdata = 'RA22_6_copy.tif'; %selects data
t = Tiff(tdata, 'r'); %unpacks Tiff format
info = imfinfo(tdata);
numImages = numel(info); %counts number of images in tiff
tnum = 259; %choose number of image in tiff to analyze
image = imread(tdata,tnum); %reads desired image
%Finds peaks of Matrix
filt = (fspecial('gaussian', 7,1)); %filter???????
avgPix = mean(image(:));
thres = 3 * avgPix; %minimum peak value
edg = 80; %skips edge pixels
%fid = fopen([filename], 'w+');
p=FastPeakFind(image, thres, filt, edg); %function that finds peaks of 2D matrix
%use found peaks to reduce noise
x = p(1:2:end); % x-coordinates of peaks
y = p(2:2:end); % y-coordinates of peaks
yLine = median(y); % y coordinate of most of the data
bounds = round(yLine - std(y)); %new bounds to get rid of noise
p=FastPeakFind(image, thres, filt, bounds, 1); %function that finds peaks of 2D matrix
x = p(1:2:end); % new x-coordinates of peaks (with less noise)
y = p(2:2:end); % new y-coordinates of peaks (with less noise)
%Use peaks to graph relevent pixel values
imgL = 1000; %length of image
pMatrix = zeros([range(y)+1, imgL]); %matrix for relevent pixel values
sVect = zeros(1,imgL); %vector for sum of x-coordinate corrisponding relevent pixel values
cVect = zeros(imgL); % 1-imgL vector just for graphing
yMin = min(y) - 1;
if range(y) < 6 %include to avoid very long imbedded for loop
for n=1:range(y)+1 %goes through all relevent y values of image
for m=1:imgL %goes through all x values of image
pMatrix(n,m) = pMatrix(n,m) + image(yMin+n,m); %collect relevent data
cVect(m) = m; %make vector for graphing
end
plot(cVect, pMatrix(n,:)) %graph each relevent row
sVect = sVect + pMatrix(n,:); %vector sum of elements of relevent rows
hold on
end
hold off
figure
plot(cVect, sVect/n) %graph average of relevent rows
end
Output:
average
each row a different color
And his one makes a graph of all the rows of pixel values with peaks, and then an average of those rows, except it zooms in on the interesting x-values, and ignores the edges so that you can see more information:
tdata = 'RA22_6_copy.tif'; %selects data
t = Tiff(tdata, 'r'); %unpacks Tiff format
info = imfinfo(tdata);
numImages = numel(info); %counts number of images in tiff
tnum = 500; %choose number of image in tiff to analyze
image = imread(tdata,tnum); %reads desired image
%Finds peaks of Matrix
filt = (fspecial('gaussian', 7,1)); %filter???????
avgPix = mean(image(:));
thres = 3 * avgPix; %minimum peak value
edg = 80; %skips edge pixels
%fid = fopen([filename], 'w+');
p=FastPeakFind(image, thres, filt, edg); %function that finds peaks of 2D matrix
%use found peaks to reduce noise
x = p(1:2:end); % x-coordinates of peaks
y = p(2:2:end); % y-coordinates of peaks
yLine = median(y); % y coordinate of most of the data
bounds = round(yLine - std(y)); %new bounds to get rid of noise
p=FastPeakFind(image, thres, filt, bounds, 1); %function that finds peaks of 2D matrix
x = p(1:2:end); % new x-coordinates of peaks (with less noise)
y = p(2:2:end); % new y-coordinates of peaks (with less noise)
%Use peaks to graph relevent pixel values
imgL = 1000; %length of image
pMatrix = zeros([range(y)+1, imgL]); %matrix for relevent pixel values
sVect = zeros(1,imgL); %vector for sum of x-coordinate corrisponding relevent pixel values
cVect = zeros(imgL); % 1-imgL vector just for graphing
yMin = min(y) - 1;
xMin = round(min(x) - std(x));
xRange = round(range(x) + 2*std(x) + 5);
if range(y) < 4 %include to avoid very long imbedded for loop
for n=1:range(y)+1 %goes through all relevent y values of image
for m=1:xRange %goes through all x values of image
pMatrix(n,m) = pMatrix(n,m) + image(yMin+n,xMin+m); %collect relevent data
cVect(m) = m; %make vector for graphing
end
plot(cVect, pMatrix(n,:)) %graph each relevent row
sVect = sVect + pMatrix(n,:); %vector sum of elements of relevent rows
hold on
end
hold off
figure
plot(cVect, sVect/n) %graph average of relevent rows
[peaks,locations,width,prominence] = findpeaks(sVect/n); %saves peaks, locations, width, and prominence
end
Output:
each row a different color
average
I'm still trying to figure out why there is that weird incline on the bottom, but besides that these have worked well for all 300 images I tested them on!!
These are helpful because you get a more precise idea of the peaks, their magnitude and their position, plus Matlab has more readily available functions to analyze graphs like this instead of graphs like the ones from last week.
I'm afraid this is pretty much where all my time went this week.
Week 10: (April 10 - April 17)
This week I worked a lot on the Matlab project, and I started making the gas chamber that I mentioned a couple weeks ago.
So for the matlab project, I needed to find a better way to show the graphs that showed the peaks well, as well as the position. Since the images are 2D, using a 2D graph does not account for the pixel values and the x and y direction, so I decided to make a 3D graph.
The advantages to the 3D graph are that it has both the colored picture with the peaks in one rotation:
Which is smilier (albeit stretched out) to the image below which I already had the code to make:
The other advantage to these graphs is they have the peaks in the other rotations:
The log filtered image is on the left, and the normal one is on the right. As you can see, the normal one has much more pronounced peaks, but with the log its easier to see all the peaks.
I did this by playing with for loops to fill up matrices with the x, y, and pixel values, making the graph. See the code below:
tdata = 'RA22_6_copy.tif'; %selects data
t = Tiff(tdata, 'r'); %unpacks Tiff format
info = imfinfo(tdata);
numImages = numel(info); %counts number of images in tiff
tnum = 500; %choose number of image in tiff to analyze
image = imread(tdata,tnum); %reads desired image
image = double(image); %reads desired image
%Finds peaks of Matrix
filt = (fspecial('gaussian', 7,1)); %filter???????
avgPix = mean(image(:));
thres = 3 * avgPix; %minimum peak value
edg = 80; %skips edge pixels
%fid = fopen([filename], 'w+');
p=FastPeakFind(image, thres, filt, edg); %function that finds peaks of 2D matrix
%use found peaks to reduce noise
x = p(1:2:end); % x-coordinates of peaks
y = p(2:2:end); % y-coordinates of peaks
yLine = median(y); % y coordinate of most of the data
bounds = round(yLine - std(y)); %new bounds to get rid of noise
p=FastPeakFind(image, thres, filt, bounds, 1); %function that finds peaks of 2D matrix
x = p(1:2:end); % new x-coordinates of peaks (with less noise)
y = p(2:2:end); % new y-coordinates of peaks (with less noise)
%Use peaks to graph relevent pixel values
imgL = 1000; %length of image
yMin = round(min(y) - 6*std(y) - 5); %gives lower bound of y
xMin = round(min(x) - 3*std(x)); %gives lower bound of x
xRange = round(range(x) + 6*std(x)); %gives range x
yRange = round(range(y) + 12*std(y) + 10); %gives range y
pMatrix1 = zeros([yRange, xRange]); %matrix for relevent log filtered pixel values
pMatrix2 = zeros([yRange, xRange]); %matrix for relevent pixel values
sVect = zeros(1,xRange); %vector for sum of x-coordinate corrisponding relevent pixel values
cVect = zeros(xRange); % 1-imgL vector just for graphing
if range(y) < 20 %include to avoid very long imbedded for loop
for n=1:yRange %goes through all relevent y values of image
for m=1:xRange %goes through all x values of image
pMatrix1(n,m) = pMatrix1(n,m) + log(image(yMin+n,xMin+m)); %collect relevent data with log filter
pMatrix2(n,m) = pMatrix2(n,m) + (image(yMin+n,xMin+m)); %collect relevent data without log filter
cVect(m) = m; %make vector for graphing
end
end
%plot log graph
h1 = subplot(1,2,1);
mesh(pMatrix1);
colormap(h1);
%plot non-log graph on same figure
h2 = subplot(1,2,2);
mesh(pMatrix2);
colormap(h2);
numPeaks = size(x) %print number of peaks
end
Next I'll talk about my work on the gas chamber:
The gas chamber is going to be a aluminum container that they'll mount on the optic table. It is basically block of aluminum with a large whole through it. I has valves for the gas pipes, tapped holes for it to be mounted, and tt'll have 3 windows, one for the laser that is ionizing the gas to come in, and another for the camera to see the laser ionize the gas.
So far, I took a block of aluminum and machined it so that it was exactly 1 1/8 square and 10 inches long. This took about 6 hours on the mill, which is a machine that I pictured above. What's useful about it is that it has a screen that shows exactly how much you've moved it in each direction. I got within .001 inches on the square end and .08 inches on the long end, which is pretty good.
Next I'm going to drill a whole through it, which will not be easy because the longer the whole you drill, the less precise it is. I might have to learn to use a lathe to do this more precisely.
Next week I'll have a picture of hopefully something a little more interesting than a block of aluminum to show you, but that's all I did this week.
Week 11: (April 17 - April 24)
This week was one of those less interesting, more time consuming weeks.
I worked more on the matlab so that it outputs matrices instead of just graphs, because we're going to need to go through a lot of data.
For example, this:
%Finds peaks of Matrix
filt = (fspecial('gaussian', 7,1)); %filter???????
avgPix = mean(image(:));
thres = 2.5 * avgPix; %minimum peak value
edg = 80; %skips edge pixels
%fid = fopen([filename], 'w+');
p=FastPeakFind(image, thres, filt, edg); %function that finds peaks of 2D matrix
%use found peaks to reduce noise
x = p(1:2:end); % x-coordinates of peaks
y = p(2:2:end); % y-coordinates of peaks
yLine = median(y); % y coordinate of most of the data
bounds = round(yLine - std(y) - 4); %new bounds to get rid of noise
p=FastPeakFind(image, thres, filt, bounds, 1); %function that finds peaks of 2D matrix
x = p(1:2:end); % new x-coordinates of peaks (with less noise)
y = p(2:2:end); % new y-coordinates of peaks (with less noise)
%Use peaks to graph relevent pixel values
imgL = 1000; %length of image
yMin = min(y)-1;
xMin = min(x)-1;
yRange = range(y) + 1;
outputMatrix = zeros([3, size(x,1)]); %matrix for relevent pixel values
if range(y) < 4 %include to avoid very long imbedded for loop
for m=1:size(x,1) %goes through all x values of image
outputMatrix(1,m) = x(m);
outputMatrix(2,m) = y(m);
outputMatrix(3,m) = image(y(m),x(m)); %collect relevant data
end
outputMatrix
end
Which outputs the following:
outputMatrix =
450 460 501 534 554 590
103 103 102 102 102 102
1496 6562 21413 31371 4489 17038
So it basically is useful for going through a lot of data really quickly and only giving you relevant information, which is important because that's what I need to do.
I'm also working on writing a function that will detect how wide a light spot from an electron avalanche is. This is important because if there are two electron avalanches right next to each other, the peaks function may not pick it up, but it would be apparent by how wide it is.
Here's my draft on that so far. It's not done (doesn't even work yet) so please don't judge!
function [widths]= FindWidth(peaksMatrix, image)
x = peaksMatrix(1,:);
y = peaksMatrix(2,:);
peaks = peaksMatrix(3,:);
%Use peaks to graph relevent pixel values
avgPix = mean(image(:));
l = size(x,2);
width = zeros(1,l);
for n=1:l
baseWidth = round(.4*peaks(n) + .6*avgPix);
m = peaks(n);
xLeft = x(n);
xRight = x(n);
while m > baseWidth
xLeft = xLeft - 1;
m = image(y(n),xLeft);
end
o = peaks(n);
while o > baseWidth
xRight = xRight + 1;
o = image(y(n),xRight);
end
width(n) = xRight - xLeft + 1;
end
widths = width;
return
Next I also started making this gas chamber, which entailed drilling a long, 3/4'' inch hole through a piece of 10" aluminum.
In the picture above, I'm using a lathe (which I had never used before) to spin the piece of aluminum while the drill piece slowly goes in it. I needed to use the lathe because the mills in our shop weren't big enough to leave space for both the piece standing on it's end, and the drill piece.
Next for this is drilling a couple smaller wholes that are going to be the windows that the laser goes through, some tapped wholes for the pipes, and some tapped holes for mounting. Luckily, this doesn't have to be finished until mid May.
Week 12: (April 24 - May 1)
This week was a weird one. I don't have much to put in my logbook, because although I spent a lot of time at work, it was spent in the lab working on setting up the laser. Much of this process I've described in previous weeks. It's a very tedious, slow process of making sure every optic reflects the light at it's center to make the laser go to the next optic—perfectly. We used irises, poles, etc. (all described in week 7).
I also worked a lot on my paper presentation, which is due tomorrow.
The time I spent in the lab, not messing with lasers, I spent on MATLAB and learning more about the overarching goal of this project so that I could actually do my presentation. I'm not technically supposed to explain that here because this experiment is not public yet, and this website is (but it will be in my presentation).
Here's the find widths code (now finished)
function [widths]= FindWidth(peaksMatrix, image)
x = peaksMatrix(1,:);
y = peaksMatrix(2,:);
peaks = peaksMatrix(3,:);
%Use peaks to graph relevent pixel values
avgPix = mean(image(:));
l = size(x,2);
width = zeros(1,l);
%for n:1:l
%avgPeaks = avgPeaks + peaks(n);
%end
for n=1:l
baseWidth = round(.3*mean(peaks) + .7*avgPix);
if peaks(n) > 2500
baseWidth = round(.4*peaks(n) + .6*avgPix);
end
m = peaks(n);
xLeft = x(n);
xRight = x(n);
while m > baseWidth
xLeft = xLeft - 1;
m = image(y(n),xLeft);
end
o = peaks(n);
while o > baseWidth
xRight = xRight + 1;
o = image(y(n),xRight);
end
width(n) = xRight - xLeft + 1;
end
widths = width;
return
Here's my current matlab code. It takes an image and returns a matrix with the coordinates, pixel values and widths (in pixels) for each peak.
%BEFORE RUNNING THIS, MUST RUN Tiff2Image
%imshow(image);
%Finds peaks of Matrix
filt = (fspecial('gaussian', 7,1)); %filter???????
avgPix = mean(image(:));
thres = 2.5 * avgPix; %minimum peak value
edg = 80; %skips edge pixels
%fid = fopen([filename], 'w+');
p=FastPeakFind(image, thres, filt, edg); %function that finds peaks of 2D matrix
%use found peaks to reduce noise
x = p(1:2:end); % x-coordinates of peaks
y = p(2:2:end); % y-coordinates of peaks
yLine = median(y); % y coordinate of most of the data
bounds = round(yLine - std(y) - 4); %new bounds to get rid of noise
p=FastPeakFind(image, thres, filt, bounds, 1); %function that finds peaks of 2D matrix
x = p(1:2:end); % new x-coordinates of peaks (with less noise)
y = p(2:2:end); % new y-coordinates of peaks (with less noise)
%plot peaks on graph
%imagesc(image); hold on
%plot(x, y, 'r+');
%Use peaks to graph relevent pixel values
imgL = 1000; %length of image
yMin = min(y)-1;
xMin = min(x)-1;
yRange = range(y) + 1;
outputMatrix = zeros([4, size(x,1)]); %matrix for relevent pixel values
if range(y) < 4 %include to avoid very long imbedded for loop
for m=1:size(x,1) %goes through all x values of image
outputMatrix(1,m) = x(m);
outputMatrix(2,m) = y(m);
outputMatrix(3,m) = image(y(m),x(m)); %collect relevent data
end
%call widths fuction to find widths
widths = FindWidth(outputMatrix, image);
%add witdths to matrix
outputMatrix(4,:) = widths(:);
%print matrix to command window
outputMatrix
end
It output the following matrix, where each column refers to a peak (from left to right), and the first row is the x coordinate, the second row is the y coordinate (both in pixels), the 3rd row is the pixel values of the peak, and the fourth row is the number of pixels wide the peak is.
This is for the following graph/image which has lots of peaks (each marked with a red cross)
I'm afraid that's it for this week. I'm looking forward to the summer where everything will go a little more quickly.
2018 Logbook
The "pwd" command tells you which directory you're in.
The "ls" command tells you what files are in the directory you're in.
The "ls -l" command gives you more information (prints the author and dates) about the files.
The "ls -l -h" gives the size of files (in a format that people can understand).
Use the "cd" command to go to a directory.
The "mkdir new_directory" command makes new directory.
The ".." following command goes back one directory for listing or opening back directory depending on command.
The "rmdir" command deletes empty directories
Use "rm" to delete files and directories, and use "rm -r" to delete directories with files contained in them.
Use "touch" to make new empty files
9/17/18
HW #2:
1) The "touch" command makes empty files if the file does not exist. However, if the file does exist it just "touches" the file, meaning it updates the time stamp to that moment.
2) Use the " . " to search through your current directory for a file.
3) Use the " * " after a string to find any file in the current directory that has that specified string in the name.
4) Use "rm" to delete files and directories, and use "rm -r" to delete directories with files contained in them. The "rmdir" command only deletes empty directories.
5) The " rm * " command would delete everything the * was referencing, which in this case would be a space, and every single file has a space, so " rm * " would be quite a big problem.
6) The ">" operator replaces the contents of a file with the new file, from the previous command. The ">>" operator adds the new file without deleting the old one.
7) The potential danger with using the " > " operator is that you could accidentally delete files when adding new ones, especially because > is so similar to >>.
8) The commands listed replaces everything in the log.txt file with a list of everything in the current directory.
Useful Linus Commands:
1) Use lsblk to list block devices by their assigned name.
2) Use uname to find information about the machine.
3) Use "history" to see the recent history of commands in the terminal.
4) Use "cal" to view a calendar
5) Use "cp" to copy files to different locations
6) Use "curl ifconfig.me" to access IP address.
7) Use tree to see content of files in a easy to understand format
8) Use arrow keys to re-run commands without retyping them
9) Use "factor" to find all possible factors of an integer
10) Use bind -p to see all the short cuts available
9/9/18
HW 3
Part A
This is a screenshot from the script “testscript.tcsh”.
PartB
Below is the output when running the script. You can see the find command has located 3 files containing test (the initial argument).
Chmod changes the mode to allow you to access directories.
PartC
In part C I created the script below to search the home directory for an argument.
The path of all files found is then stored in homesearch_$ARG.txt. Below is an example of the script functioning properly.
PartD
Piping in linux is a way to take the output from one command or script as the input for another command or script.
PartE
With the following command, I replaced every “was” with “is”.
Additional problems:
October 2nd, 2018
HW 4
Use emacs to write C++ code
Use “g++ main.cpp” to compile and tanslate code
Use “./a/out” to run and print main.cpp
Use “g++ -dumpspecs” to see all the information find information about which C++ version I’m using, under the heading “version” (see below). There’s also lots of other information in there.
Use “g++ -dumpspecs” shows version that I’m using:
Use “ls /usr/include/c++/4.4.4/” to see all the commands you can use
The use ls /usr/include/math.h and ls /usr/include/bits/*math*.h see
Variables:
Numeric:
//tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
//"cout" uses a command that is in the "iostream" group, use to print hello world to screen, "endl" prints an end line
cout <<"Hello World!" << endl;
int i=2; //assign the value 2 to the variable i
cout << "i = " <<i<<endl; //prints the value of i to screen
double a=3.3; //assign the value 3.3 to the variable a. Unlike an int, a double is aloud to have multiple decimal places
cout <<"a = " <<a<<endl; //print the value of a to screen
int j = a*i; //assign the value of a multiplied by the value of i to the variable j. Note that it will be 6 instead of 6.6 because it will round down in order to have an integer value, because j is defined as an int
cout<< "a*i = " << j << endl; //print value of j to screen
return 0; //Exit the program
}
//semicolons used to separate commands
The C++ above prints this. Note that a*i is equal to 6 even though it should technically be 6.6 because the variable j that a*i is assigned to is a integer and can’t accommodate the double value of 6.6
//tell compiler that we're using C++ commands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int n=10;//assign value 10 to variable n
cout << "n is "<<n<<endl;//print n to the screen
n--;//equivalent to n=n-1
cout<<"n is now "<<n<<endl;//print new value of n to screen
n++;//equivalent to n=n+1
cout<<"n is now "<<n<<endl;//print new value of n to screen
return 0; //Exit the program
}
Non Numeric
/tell compiler that we're using C++ commands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
bool prop; //defines a boolean variable named prop
prop = (5>1);//function is true, so assigns value of 1 to prop
cout<<"prop is "<<prop<<endl;//prints props value
prop = (1>5);//function is false, so assigns value of 0 to prop
cout<<"prop is "<<prop<<endl;//prints props value
prop = (1 != 5);//function is true, so assigns value of 1 to prop
cout << "prop is " <<prop<<endl;//prints props value
return 0; //Exit the program
}
Loops 1
/tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int n=10;//set integer n to have a value of 10
while(n>0) { //set up a loop that runs while n is greater than zero
cout<<"n is "<<n<<endl;//print the value of n to screen
n--;//equivilant to n=n-1
}
return 0; //Exit the program
}
While loop starts at n=10, and then prints value for each loop until n=0 and then it exits.
Loops 2:
Re-write following, except replacing the while loops with for loops:
#include <iostream>
using namespace std;
int main() {
int n=0, m=0;
while(n<10) {
// this is the slow (or outer) loop
cout << "n is " << n << ": ";
m=0;
while(m<=n) {
// this is the fast (or inner) loop
// in this loop, the slow loop variable (n) is a constant
// this loop must run to completion before the slow loop
// can progress (during every iteration of the slow loop!)
cout << m;
m++;
}
// now the fast loop has finished and the slow loop can
// continue with the current iteration
cout << endl;
n++;
}
return 0 ;
}
Replace while loops with for loops and move all the information necessary into for loop requirements
//tell compiler that we're using C++ commands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
for(int n=0;n<10;n++) {
// this is the slow (or outer) loop
cout << "n is " << n << ": ";
for(int m=0;m<=n;m++) {
// this is the fast (or inner) loop
// in this loop, the slow loop variable (n) is a constant
// this loop must run to completion before the slow loop
// can progress (during every iteration of the slow loop!)
cout << m;
}
// now the fast loop has finished and the slow loop can
// continue with the current iteration
cout << endl;
}
return 0; //Exit the program
}
10/15/18
Logical statements:
//tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int n = 10;//assign integer variable n the value of 10
while (n>=10)//create a loop that continues to run as loon as n is greater or equal to 10
{
if(n>5) { //creates a logical statement so that the code will only do the following command if n is greater than 5 is true
cout<<"n is "<<n<<endl;//prints the value of n
}
else { //creates a logical statement so that the code will only do the following command if n is less than 5 is false
cout<<"n = "<<n<<endl;//prints the value of n
}
n--;//changes the value of n (had to put this outside the if/else values, because if it were inside it would have created an infinite loop, meaning the while statement would always be true and never exit
}
return 0; //Exit the program
}
//semicolons used to separate commands
Pointers:
//tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int i = 10;//assign i the value of 10
cout << "The memory address of i is " << &i << "\n";//print the memory address of i
cout << "The data stored at memory address " << &i << " is " << i << "\n";//prints the value assigned to that memory address
int* p = &i;//assignes the memory address of i to the variable p
cout << "The value of p is " << p << "\n";//prints the value of p as the memory address of i
cout << "We say that p 'points at' the memory location referenced by address " << p << "\n";//explains that p is the location (or points to the location/address) of i
cout << "The data stored at memory address " << p << " is " << *p << "\n";//prints that the value of i (10) is still stored at the location/address of i and/or the address that p is assigned to (or points at).
return 0; //Exit the program
}
PROGRAM #1:
//tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int i = 10;//assign i the value of 10
int j = i;//assign j the value of 10
cout << "i= " << i << " and j= " << j << "\n";//print i and j 's values
i=5;//reassign i the value of 5
cout << "i= " << i << " and j= " << j << "\n";//print i and j 's new values
j=1;//reassign j the value of 5
cout << "i= " << i << " and j= " << j << "\n";//print i and j's new values
return 0; //Exit the program
}
PROGRAM #2:
//tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int i = 10;//assign i the value of 10
int* p = &i;//makes p point to i and the value of 10
cout << "i= " << i << " and *p= " << *p << "\n";//print i and *p 's values
i=5;//reassign i the value of 5
cout << "i= " << i << " and *p= " << *p << "\n";//print i and *p 's new values (because p points to i's value, *p and i will always have the same value, even when one or the other is reassigned)
*p=1;//reassign *p the value of 1
cout << "i= " << i << " and *p= " << *p << "\n";//print i and *p's new values (because p points to i's value, *p and i will always have the same value, even when one or the other is reassigned)
return 0; //Exit the program
}
Assigning p to point to a random memory address with a specified value
//tell compiler that we're using C++ comands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int* p = new int(5);//assings p to point to a memory address with the value 5
cout << "p points at address " << p << "\n";//prints the memory address
cout << "The data stored in address " << p << " is " << *p << "\n";//prints the data stored at the memory address
*p = 10;//changes the data stored at the memory addres to 10
cout << "Now the data stored in address " << p << " is " << *p << "\n";//prints the new data stored at the memory address
return 0; //Exit the program
}
HW#5
Problem 1:
Projected range for copper = 4.739E+02 g/cm2 assuming 1 geV Kinetic Energy
To convert to cm, multiply by reciprocal of copper’s density: 8.96 g/cc
4.739E+02 g/cm2 * 1/(8.96 g/cc) = 52.89 cm
So it takes 52.89 cm of copper to stop a proton with 1 geV of Kinetic Energy
Problem 2:
Computing:
Problem 1:
//tell compiler that we're using C++ commands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int x = 17;//assign the value 17 to x
cout << "x = " << x << "\n";//prints the value of x
while (x>=4)//create loop that runs as long as x is greater or equal to 4
{
if (x>=10)//logical statement that runs only if x is greater or equal than 10
{
x = 14 - x; //sets x to 14 minus its old value
cout << "if x is greater or equal to 10, set x = 14 - x(old) = " << x << "\n";//prints the value of new x acd explains what happened
}
else //logical statement that runs the following commands only if x is less than 10
{
x--;//subtracts one from x
cout << "if x is less than 10, subtract one from x such that x = " << x << "\n";//prints the value of new x and explains what happened
}
}
return 0; //Exit the program
}
//semicolons used to separate commands
Ran twice. Once with x = 7 and once with x = 17
Problem 2:
#!/bin/tcsh
set month = `date | awk '{print $2}'`
echo $month
ls -lat -R ../| grep "$month"
October 28, 2018
HW 6
Arrays in C++
//tell compiler that we're using C++ commands to send output to terminal/printer/file
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main()
{
int ii[3] = {1,2,3}; //define the array ii to have 3 elements equal to {1,2,3}
int j=0; //define the integer variable j equal to zero
while (j<3) { //set up while loop for as long as j is less than 3
cout <<" ii of "<<j<<" is "<<ii[j]<<endl; //print j's value to terminal
j++; //add one to j and go to the next element of array
}
int LL[2][3] = {1,2,3,4,5,6}; //define the 2D array LL to have two columns and three rows, where the first row is {1,2} the second row is {3,4} and the third row is {5,6}.
j=0; //set the variable j back equal to zero
int k; //define an integer variable k
while (j<2) { //set up a while loop that runs as long as j is less than two (because LL has two columns). This while loop will iterate through the columns.
k=0;//set k equal to zero
while (k<3) { //set up another while loop (inside the first one) to go through the rows. It runs as long as k is less than three because there are 3 rows.
cout<<" LL of "<<j<<" "<<k<<" is "<<LL[j][k]<<endl;//print the value of LL at the specified row and column location.
k++; //add one to k so that you go to the next row
}
j++; //add one to j so that you go to the next column
}
return 0; //Exit the program
}
//semicolons used to separate commands
Reading data from a file:
//tell compiler that we're using C++ commands to send output to terminal/printer/file. Used for cin/cout
#include <iostream>
//Used for reading and writing of files. Includes ofstream, which is used below. without this header it would not compile
#include <fstream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main() {
ofstream myfile;//makes it possible to use fstream to edit files by making a new ofstream class called myfile
myfile.open("example.txt");//uses the open command from the ofstream class myfile to open the example.txt file. If example.txt does not exist it will create it.
myfile<<"write some junk.";//"<<" means write the following text into the file, so now the example.txt file will include the phrase "write some junk".
myfile.close();//closes and saves the file
return 0;//Exit the program
}
Declarations and combining code from different files
/* MAIN.CPP */
//tell compiler that we're using C++ commands to send output to terminal/printer/file. Used for cin/cout
#include <iostream>
//Used for reading and writing of files. Includes ifstream, which is used below. without this header it would not compile
#include <fstream>
// include the program dotprod.cpp so that we can find and use the dot_prod function
#include "dotprod.cpp"
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main () {
// declare the vectors, each with three elements
double vector1[3];
double vector2[3];
// open the input file
ifstream infile; //creates a ifstream type class called infile
infile.open("vectors.txt"); //uses ifstream class to call open command on vectors.txt. Opens vectors.txt because those are the vectors that we are using
// store the input in the vectors and print the vectors for the user
infile>>vector1[0]>>vector1[1]>>vector1[2]; //">>" means read, so this is just reading vector1 and each of its elements from vectors.txt
cout<<" Vector 1 is ("<<vector1[0]<<","<<vector1[1]<<","<<vector1[2]<<")"<<endl;//this is printing the vector1, showing each of its elements
infile>>vector2[0]>>vector2[1]>>vector2[2];// reads vector2 and each of it's elements
cout<<" Vector 2 is ("<<vector2[0]<<","<<vector2[1]<<","<<vector2[2]<<")"<<endl;//prints vector2 and each of it's elements
// close the input file
infile.close();
// call the dot_prod function from dotprod.cpp
dot_prod(vector1,vector2);
return 0;//exit the program
}
Create dot product program with dot product function
/* DOTPROD.CPP */
//tell compiler that we're using C++ commands to send output to terminal/printer/file. Used for cin/cout
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
double dot_prod(double v1[3],double v2[3]) { //define function dot_prod to take two vectors, each with three double type elements, and output a double
double dotdot; //defines variable dotdot
dotdot = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];//set variable dotdot equal to the computed dot product, the sum of the product of each corresponding element in two vectors
cout<<" The dot product is "<<dotdot<<endl;//prints the value of dotdot, which is now the value of the dot product
return 0;//exit the program
}
Create file with vectors to use in dot product (the first row used as vector1, second row used as vector2, and third row used as scalar)
1 2 3
4 5 6
7
My modifications:
main.cpp
/* MAIN.CPP */
//tell compiler that we're using C++ commands to send output to terminal/printer/file. Used for cin/cout
#include <iostream>
//Used for reading and writing of files. Includes ifstream, which is used below. without this header it would not compile
#include <fstream>
// include the program dotprod.cpp so that we can find and use the dot_prod function
#include "dotprod.cpp"
// include the program scalarmult.cpp so that we can find and use the scalar_mult function
#include "scalarmult.cpp"
//Define functions to use in the program (like the String function variable for example)
using namespace std;
//tell compiler where first line of code is
int main () {
// declare the vectors, each with three elements
double vector1[3];
double vector2[3];
double scalar;
// open the input file
ifstream infile; //creates a ifstream type class called infile
infile.open("vectors.txt"); //uses ifstream class to call open command on vectors.txt. Opens vectors.txt because those are the vectors that we are using
// store the input in the vectors and print the vectors for the user
infile>>vector1[0]>>vector1[1]>>vector1[2]; //">>" means read, so this is just reading vector1 and each of its elements from vectors.txt
cout<<" Vector 1 is ("<<vector1[0]<<","<<vector1[1]<<","<<vector1[2]<<")"<<endl;//this is printing the vector1, showing each of its elements
infile>>vector2[0]>>vector2[1]>>vector2[2];// reads vector2 and each of it's elements
cout<<" Vector 2 is ("<<vector2[0]<<","<<vector2[1]<<","<<vector2[2]<<")"<<endl;//prints vector2 and each of it's elements
infile>>scalar;//reads the scalar value in the file
cout<<"The scalar is "<<scalar<<endl;//prints the scalar value to terminal
// close the input file
infile.close();
// call the dot_prod function from dotprod.cpp
dot_prod(vector1,vector2);
// call the scalar_mult function from scalarmult.cpp
scalar_mult(vector1,scalar);
scalar_mult(vector2,scalar);
return 0;//exit the program
}
Scalarmult.cpp (new)
/* SCALARMULT.CPP */
//tell compiler that we're using C++ commands to send output to terminal/printer/file. Used for cin/cout
#include <iostream>
//Define functions to use in the program (like the String function variable for example)
using namespace std;
double scalar_mult(double v1[3],double x) {//define function scalar_mult to take a vector with three double type elements and one scalar
double mult[3];//define variable mult to be a vector with three elements of type double
mult = {x*v1[0],x*v1[1],x*v1[2]};//set mult equal to the computed scalar multiple, a vector where each of the elements are multiplied by the scalar
cout<<"The scalar multiple is ("<<mult[0]<<","<<mult[1]<<","<<mult[2]<<")"<<endl;//print the scalar multiple with each of its components
return 0;//exit program
}
Random numbers:
#include <iostream>
#include <math.h>
using namespace std;
// this function is the actual random number generator
// this code is stolen from the book numerical recipes for fortran
// it relies on random generated from overflow of memory locations
// and is a pseudo random number generator
const int a = 7141;
const int c = 54773;
const int mmod=256200;
double getFlatRandom(int& inew) {//the & makes it so you can get a new integer because it passes inew by reference, so when it changes inside the function, it also changes outside the function
double mranflat = 0.;
//the following changes inew so that it gives pseudo-random numbers
inew = inew%mmod;
double aa = double(inew)/double(mmod);
mranflat=aa;
inew = a*inew+c;
return mranflat;
}
// in this code, we will call the pseudo-random number generator and learn some things
// about its properties by filling and then displaying a histogram
int main() {
int num;
cout << "Enter the number of loop iterations: ";
cin >> num;
int inew = 2345; // This is the "seed" for the random number generator
// we will put the results from the call into a histogram that we can look at, to learn some of its
// properties. This histogram has 10 bins.
int histo[10] = {0,0,0,0,0,0,0,0,0,0};
double atmp; //
// call the random number generator 1000 times and fill a histogram
for (int i = 0; i<num; i++){
atmp=getFlatRandom(inew); // call the random number generator
histo[int(atmp*10)]++; // increment the histogram bin the number falls within
if(i<10) //this makes the following statements only apply the the first times
cout <<int(atmp*10)<<endl; //prints number randomly selected
}
// print the histogram to the screen
for (int i = 0; i<10; i++){
cout<<i<<": ";
for (int j=0; j<int((double(100*histo[i])/double(num))+0.5); j++) {//this multiplies the display to get 100 data points, even if there are not that many
cout << "=";
}
cout << endl;
}
return 0;
}
First run with 10 loop iterations:
Then took out the & and ran with 10 loop iterations:
The & makes it so you can get a new integer because it passes inew by reference, so when it changes inside the function, it also changes outside the function
A bit has two combinations, and each byte has 8 bits in it, so a byte would have 2^8 combinations. Then, because the binary starts at zero, you subtract one so the combination of n bytes is 2^8 * 2^n -1 = 2^(8n) -1.
When I added the following code (highlighted in context above) it showed that each loop outputs the same numbers because the seed is the same and the “random” numbers are being made from that same seed with the same formula.
if(i<10) //this makes the following statements only apply the the first times
cout <<int(atmp*10)<<endl; //prints number randomly selected
The result is shown below, each with the same seed.
[cms-opendata@localhost CMSSW_5_3_32]$ emacs random.cpp
[cms-opendata@localhost CMSSW_5_3_32]$ g++ random.cpp
[cms-opendata@localhost CMSSW_5_3_32]$ ./a.out random.cpp
Enter the number of loop iterations: 13
0
5
1
5
6
6
7
0
6
4
0: ===============
1: ========
2:
3:
4: ========
5: ===============================
6: =======================
7: ===============
8:
9:
[cms-opendata@localhost CMSSW_5_3_32]$ ./a.out random.cpp
Enter the number of loop iterations: 13
0
5
1
5
6
6
7
0
6
4
0: ===============
1: ========
2:
3:
4: ========
5: ===============================
6: =======================
7: ===============
8:
9:
If I change the seed, I get something different:
[cms-opendata@localhost CMSSW_5_3_32]$ emacs random.cpp
[cms-opendata@localhost CMSSW_5_3_32]$ g++ random.cpp
[cms-opendata@localhost CMSSW_5_3_32]$ ./a.out random.cpp
Enter the number of loop iterations: 13
0
5
5
8
5
0
8
9
8
7
0: ===============
1: ========
2:
3: ========
4:
5: =======================
6: ========
7: ========
8: =======================
9: ========
[cms-opendata@localhost CMSSW_5_3_32]$
As you increase the number of iterations, the histogram levels out.
[cms-opendata@localhost CMSSW_5_3_32]$ ./a.out random.cpp
Enter the number of loop iterations: 100
0
5
5
8
5
0
8
9
8
7
0: ========
1: ============
2: ==========
3: =======
4: ==========
5: ===========
6: =========
7: ===========
8: ============
9: ==========
[cms-opendata@localhost CMSSW_5_3_32]$ ./a.out random.cpp
Enter the number of loop iterations: 1000
0
5
5
8
5
0
8
9
8
7
0: ==========
1: ===========
2: =========
3: ==========
4: ==========
5: ==========
6: ==========
7: ==========
8: ===========
9: ==========
[cms-opendata@localhost CMSSW_5_3_32]$ ./a.out random.cpp
Enter the number of loop iterations: 100000
0
5
5
8
5
0
8
9
8
7
0: ==========
1: ==========
2: ==========
3: ==========
4: ==========
5: ==========
6: ==========
7: ==========
8: ==========
9: ==========
[cms-opendata@localhost CMSSW_5_3_32]$
Other (better) way to generate random numbers:
#include <stdio.h>
#include <stdlib.h> /* srand, rand */
#include <time.h>
#include <iostream>
using namespace std;
int main(){
cout << "Random Number" << rand()%100; /* This will print a random number between [0,100). */
return 0;
}
Calorimeter
#include "Riostream.h"
#include <iostream>
#include <fstream>
/*
some code to help us understand why and how resolutions affect our ability to measure the
mass of a particle
*/
// main routine
void resolutions(){
// set debug mode
bool idebug=0;
// set up the display
gROOT->Reset();
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0);
gStyle->SetOptFit(11);
gStyle->SetErrorX(0);
c1 = new TCanvas("c1","Our data",200,10,700,500);
c1->SetFillColor(0);
c1->SetBorderMode(0);
c1->SetBorderSize(1);
c1->SetFrameBorderMode(0);
// book some histograms
TH1F *histo1 = new TH1F("histo1","mass",1000,0, 200);
// setup random number generator
gRandom->SetSeed();
// get the secret parameters
ifstream myfile;
myfile.open("secretparameters.txt");
double mass, resE, resA;
myfile>>mass>>resE>>resA;
myfile.close();
// make N bosons at rest
// each boson will decay to 2 back-to-back daughter particles
// they will be back-to-back by conservation of momentum, since the momentum of the mother was zero
// assume that the mass of the daughter particles is small compared to the mother mass so we can
// assume that their energy will be large compared to their mass and we can treat them as massless.
// and thus their energy is equal to the magnitude of their momenta
// by conservations of energy, their energy
// work in 2D. Can always choose my coordinate system so that the 2 daughters are in the same plane
int N=1000;
double etrue1,etrue2,phitrue1,phitrue2; // true energy of the 2 daughters
double e1,px1,py1,phi1; // smeared 4-momenta of daughter 1
double e2,px2,py2,phi2; // smeared 4-momenta of daughter 2
double masssmeared;
for(int i=0;i<N;i++) {
// set true energy
etrue1=mass/2.;
etrue2=mass/2.;
if(idebug) cout<<"etrue "<<etrue1<<" "<<etrue2<<endl;
// choose phi for daughter 1 and daughter 2
phitrue1=2*TMath::Pi()*gRandom->Rndm();
phitrue2 = phitrue1+TMath::Pi();
if(idebug) cout<<"phitrue "<<phitrue1<<" "<<phitrue2<<endl;
// smear true energy with resolution of detector
e1=etrue1+resE*gRandom->Gaus(0.,1.);
e2=etrue2+resE*gRandom->Gaus(0.,1.);
if(idebug) cout<<"e "<<e1<<" "<<e2<<endl;
//smear angles with resolution of the detector
phi1=phitrue1+resA*gRandom->Gaus(0.,1.);
phi2=phitrue2+resA*gRandom->Gaus(0.,1.);
if(idebug) cout<<"phi "<<phi1<<" "<<phi2<<endl;
//calculate 4 momenta after smearing
px1=e1*cos(phi1);
py1=e1*sin(phi1);
px2=e2*cos(phi2);
py2=e2*sin(phi2);
if(idebug) cout<<"pxs "<<px1<<" "<<py1<<" "<<px2<<" "<<py2<<endl;
// calculate smeared mass
masssmeared=sqrt((e1+e2)*(e1+e2) - (px1+px2)*(px1+px2) - (py1+py2)*(py1+py2));
if(idebug) cout<<"masssmeared "<<masssmeared<<endl;
histo1->Fill(masssmeared);
histo1->Draw("");
c1->Update();
}
c1->SaveAs("c1.gif");
}
With N=1, the mass should be around 104 eV. It’s hard to tell how accurate this is because it only has one data point.
With N=10, it’s a little easier to see what the center point should be, but still quite confusing.
With N=100, you can see the beginning of a gaussian distribution
With N=1000, you get a clear picture of actual mass, as well as it’s distribution. I would guess from this graph that the mass is about 91, which is what we decided the mass would be in the other file.
After changing the true mass to 1, it is clear that even with lots of data points, you can’t have a gaussian distribution because the mass cannot be negative. Below is N=1.
Pictured below is N=10, with still not enough data points to get any real information.
Below is N=100, where you can see the distribution, but it’s very off.
Finally, below you can see what should be a nice distribution, but it’s wrong because mass cannot be negative so the histogram is skewed.
Accelerator Problem:
1)
2)
Ela Rockafellow
November 23, 2010
HW 7
Use madgraph to generate collisions with the code below
Use the following code to display the collisions
Below are the particles the current model has:
Below are the displayed multi particles:
The results are then available:
Finally, you have data available:
HOMEWORK:
The beam pipe is where the particles travel and (in the detector) collide
The tracker is one meter in radius and is used to find the paths the particles take through the beam pipe
The calorimeter is 2.5 meters in radius and is used to measure the energy of the particles as they pass through it.
The solenoid is 3.5 meters in radius and is used to create a 4 tesla uniform magnetic field.
The muon tracker is on the outside of the detector, 7 meters in radius, and is used to measure the energies of the muons. The calorimeter does not stop muons, so the muon tracker has extra chambers that are what measures their energy.
Use the cursor to click on the axes and change to xy perspective to look along the beam pipe. The electrons are bright green lines
Muons are long curved red lines
Missing energy looks like dashed purple lines
Jets are yellow cones coming out
Yes, the magnetic field is pointing into the page so you can use that to tell that a particle with a clockwise curve will have a positive charge, and vise versa.
Electron.ig
Events/Run_146644/Event_718940510 [1 of 25]
Vertices - 1
Electrons - 1
Muons - 0
Photons - 1
Jets - 13
Missing Energy - 1
Events/Run_146644/Event_718969382 [4 of 25]
Vertices - 4
Electrons - 1
Muons - 0
Photons - 1
Jets - 8
Missing Energy - 1
Mu.ig
Events/Run_146436/Event_90625265 [1 of 25]
Vertices - 1
Electrons - 0
Muons - 2
Photons - 0
Jets - 0
Missing Energy - 1
Events/Run_146436/Event_90882608 [24 of 25]
Vertices - 1
Electrons - 0
Muons - 1
Photons - 0
Jets - 10
Missing Energy - 0
MinimumBias.ig
Events/Run_148031/Event_441572355 [1 of 25]
Vertices - 1
Electrons - 0
Muons - 0
Photons - 0
Jets - 14
Missing Energy - 0
Events/Run_148031/Event_441641307 [12 of 25]
Vertices - 1
Electrons - 0
Muons - 0
Photons - 0
Jets - 1
Missing Energy - 0
Diphoton.ig
Events/Run_207924/Event_315722300 [10 of 10]
Vertices - 11
Electrons - 0
Muons - 0
Photons - 2
Jets - 22
Missing Energy - 1
Events/Run_202087/Event_923352992 [4 of 10]
Vertices - 6
Electrons - 0
Muons - 0
Photons - 2
Jets - 31
Missing Energy - 1
Leption.ig
Events/Run_178424/Event_666626491 [1 of 3]
Vertices - 6
Electrons - 0
Muons - 4
Photons - 0
Jets - 52
Missing Energy - 1
Events/Run_195099/Event_137440354 [3 of 3]
Vertices - 11
Electrons - 2
Muons - 3
Photons - 2
Jets - 94
Missing Energy - 1
HW#8:
Downloaded 4 root files:
wget http://hepcms-hn.umd.edu/~jabeen/Analysis/LHC-Higgs-Graviton.tgz
Opened one of them
Use the following command to make a class called "HiggsAnalysis" that you can run the Higgs analysis on:
HZZ4LeptonsAnalysisReduced->MakeClass("HiggsAnalysis")
Replace code in HiggsAnalysis with:
//The following code combines vectors to find mass, using TLorentzVectors
Edited code highlighted
#define HiggsAnalysis_cxx
#include "HiggsAnalysis.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void HiggsAnalysis::Loop()
{
// In a ROOT session, you can do:
// Root > .L HiggsAnalysis.C
// Root > HiggsAnalysis t
// Root > t.GetEntry(12); // Fill t data members with entry number 12
// Root > t.Show(); // Show values of entry 12
// Root > t.Show(16); // Read and show values of entry 16
// Root > t.Loop(); // Loop on all entries
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus("branchname",1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
if (fChain == 0) return;//gets every single event
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
TFile* output = TFile::Open("Dielectron_MC.root", "RECREATE"); // "RECREATE" would produce a new root file with name Dielectron_MC.root every time you run the code
TH1F* Z_ee = new TH1F("Z_ee", "Di-electron candidate invariant mass", 200, 0, 200);//histogram
TH1F* H_zz = new TH1F("H_zz", "ZZ candidate invariant mass", 200, 0, 300);//higgs histogram
double el1mt = 0.0;//defining double variables and setting them to zero
double el1pt = 0.0;
double el1eta = 0.0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
// if (Cut(ientry) < 0) continue;
TLorentzVector el1, el2, el3, el4;//One of the most useful classes in root (Lorentz vector) It combines relativity and vectors
el1.SetPtEtaPhiM(f_lept1_pt, f_lept1_eta, f_lept1_phi, 0.0); //define Lorentz 4-vector for electron 1 and sets PtEta values to four vectors (with lepton 1 point, lepton 1 eta, lepton 1 phi, and 0.0)
el2.SetPtEtaPhiM(f_lept2_pt, f_lept2_eta, f_lept2_phi, 0.0);//define lorantz 4-vector for electron 2 and does the same thing as above for lepton 2
el3.SetPtEtaPhiM(f_lept3_pt, f_lept3_eta, f_lept3_phi, 0.0);//define lorentz 4-vector for electron 3 and does the same thing as above for lepton 2
el4.SetPtEtaPhiM(f_lept4_pt, f_lept4_eta, f_lept4_phi, 0.0);//define lorentz 4-vector for electron 4 and does the same thing as above for lepton 2
TLorentzVector zCandidate = el1 + el2; //sets lorentz 4-vector canidate to the sum of el1 and el2 lorentz 4-vectors.
TLorentzVector zCandidate2 = el3 + el4;//sets lorentz 4-vector canidate 2 to the sum of el3 and el4 lorentz 4-vectors
TLorentzVector Higgs = zCandidate + zCandidate2;//sets lorentz 4-vector Higgs to the sum of the candidate lorentz 4-vectors
/////Find mass///////
H_zz->Fill(Higgs.M()); //fills Higgs Historgram with Higgs data
Z_ee->Fill(zCandidate.M()); //Fills historgram with candidate data
el1mt = el1.Mt(); //calling mt method from el1 object
cout << el1mt << endl; //print (out puts mass values)
}
Z_ee->Write(); //writes historgram
H_zz->Write(); //writes histogram
output->Close();
}
Resulting Histograms:
HIGGS MASS!!!!
According to this model, higgs mass = 124.1
HW9
HW10
Upload this python file:
# -*- coding: utf-8 -*- """Untitled0.ipynb Automatically generated by Colaboratory. Original file is located at https://colab.research.google.com/drive/16SMrsuUKlXt8QYA1wNHDk-m8Vlz3ufCr """ '''Deep Dreaming in Keras. Run the script with: ``` python deep_dream.py path_to_your_base_image.jpg prefix_for_results ``` e.g.: ``` python deep_dream.py img/mypic.jpg results/dream ``` ''' from __future__ import print_function from keras.preprocessing.image import load_img, save_img, img_to_array import numpy as np import scipy import argparse from keras.applications import inception_v3 from keras import backend as K parser = argparse.ArgumentParser(description='Deep Dreams with Keras.') parser.add_argument('base_image_path', metavar='base', type=str, help='Path to the image to transform.') parser.add_argument('result_prefix', metavar='res_prefix', type=str, help='Prefix for the saved results.') args = parser.parse_args() base_image_path = args.base_image_path result_prefix = args.result_prefix # These are the names of the layers # for which we try to maximize activation, # as well as their weight in the final loss # we try to maximize. # You can tweak these setting to obtain new visual effects. settings = { 'features': { 'mixed2': 0.2, 'mixed3': 0.5, 'mixed4': 2., 'mixed5': 1.5, }, } def preprocess_image(image_path): # Util function to open, resize and format pictures # into appropriate tensors. img = load_img(image_path) img = img_to_array(img) img = np.expand_dims(img, axis=0) img = inception_v3.preprocess_input(img) return img def deprocess_image(x): # Util function to convert a tensor into a valid image. if K.image_data_format() == 'channels_first': x = x.reshape((3, x.shape[2], x.shape[3])) x = x.transpose((1, 2, 0)) else: x = x.reshape((x.shape[1], x.shape[2], 3)) x /= 2. x += 0.5 x *= 255. x = np.clip(x, 0, 255).astype('uint8') return x K.set_learning_phase(0) # Build the InceptionV3 network with our placeholder. # The model will be loaded with pre-trained ImageNet weights. model = inception_v3.InceptionV3(weights='imagenet', include_top=False) dream = model.input print('Model loaded.') # Get the symbolic outputs of each "key" layer (we gave them unique names). layer_dict = dict([(layer.name, layer) for layer in model.layers]) # Define the loss. loss = K.variable(0.) for layer_name in settings['features']: # Add the L2 norm of the features of a layer to the loss. if layer_name not in layer_dict: raise ValueError('Layer ' + layer_name + ' not found in model.') coeff = settings['features'][layer_name] x = layer_dict[layer_name].output # We avoid border artifacts by only involving non-border pixels in the loss. scaling = K.prod(K.cast(K.shape(x), 'float32')) if K.image_data_format() == 'channels_first': loss += coeff * K.sum(K.square(x[:, :, 2: -2, 2: -2])) / scaling else: loss += coeff * K.sum(K.square(x[:, 2: -2, 2: -2, :])) / scaling # Compute the gradients of the dream wrt the loss. grads = K.gradients(loss, dream)[0] # Normalize gradients. grads /= K.maximum(K.mean(K.abs(grads)), K.epsilon()) # Set up function to retrieve the value # of the loss and gradients given an input image. outputs = [loss, grads] fetch_loss_and_grads = K.function([dream], outputs) def eval_loss_and_grads(x): outs = fetch_loss_and_grads([x]) loss_value = outs[0] grad_values = outs[1] return loss_value, grad_values def resize_img(img, size): img = np.copy(img) if K.image_data_format() == 'channels_first': factors = (1, 1, float(size[0]) / img.shape[2], float(size[1]) / img.shape[3]) else: factors = (1, float(size[0]) / img.shape[1], float(size[1]) / img.shape[2], 1) return scipy.ndimage.zoom(img, factors, order=1) def gradient_ascent(x, iterations, step, max_loss=None): for i in range(iterations): loss_value, grad_values = eval_loss_and_grads(x) if max_loss is not None and loss_value > max_loss: break print('..Loss value at', i, ':', loss_value) x += step * grad_values return x """Process: - Load the original image. - Define a number of processing scales (i.e. image shapes), from smallest to largest. - Resize the original image to the smallest scale. - For every scale, starting with the smallest (i.e. current one): - Run gradient ascent - Upscale image to the next scale - Reinject the detail that was lost at upscaling time - Stop when we are back to the original size. To obtain the detail lost during upscaling, we simply take the original image, shrink it down, upscale it, and compare the result to the (resized) original image. """ # Playing with these hyperparameters will also allow you to achieve new effects step = 0.01 # Gradient ascent step size num_octave = 3 # Number of scales at which to run gradient ascent octave_scale = 1.4 # Size ratio between scales iterations = 20 # Number of ascent steps per scale max_loss = 10. img = preprocess_image(base_image_path) if K.image_data_format() == 'channels_first': original_shape = img.shape[2:] else: original_shape = img.shape[1:3] successive_shapes = [original_shape] for i in range(1, num_octave): shape = tuple([int(dim / (octave_scale ** i)) for dim in original_shape]) successive_shapes.append(shape) successive_shapes = successive_shapes[::-1] original_img = np.copy(img) shrunk_original_img = resize_img(img, successive_shapes[0]) for shape in successive_shapes: print('Processing image shape', shape) img = resize_img(img, shape) img = gradient_ascent(img, iterations=iterations, step=step, max_loss=max_loss) upscaled_shrunk_original_img = resize_img(shrunk_original_img, shape) same_size_original = resize_img(original_img, shape) lost_detail = same_size_original - upscaled_shrunk_original_img img += lost_detail shrunk_original_img = resize_img(original_img, shape) save_img(result_prefix + '.png', deprocess_image(np.copy(img)))
Then:
And produces the picture:
Page updated
Report abuse