Commit 1d1376cf authored by Breno Rilho Lemos's avatar Breno Rilho Lemos 💬
Browse files

Replace sol-analitica-tpc.c

parent eaab12fe
......@@ -2,20 +2,26 @@
#include <stdlib.h>
#include <math.h>
#define PI M_PI
#define PI 3.141592
#define TPCRADIUS 2.466
#define c 3.0E8
// Reads file with data in the following order:
//
// number of tracks, vx, vy, vz, magnetic field => esd-summary
// mass, charge, px, py, pz => esd-detail
// mass, charge, px, py, pz, Vx, Vy, Vz => esd-detail
//
// For now, we are considering classical mechanics, so p = mv (without relativistic considerations)
//
// REMEMBER: Vx is for velocity; vx is for vertex.
//
void main() {
int i, j, numberOfTracks;
double vx, vy, vz, w, Rx, Ry, Rz, B, radius, t, dt;
double vx, vy, vz, w, Rx, Ry, Rz, B, radius, prevRadius, t, dt, skip[8];
double gamaX[2], gamaY[2];
FILE *detail, *summary;
......@@ -23,41 +29,60 @@ void main() {
summary = fopen("esd-summary.dat","r"); // filename must be "esd-summary.dat"
detail = fopen("esd-detail.dat","r"); // filename must be "esd-detail.dat"
do { // Reads summary file until number of tracks is different than zero (there are tracks in the event)
i=0;
j=0;
do { // Reads summary file until number of tracks is between 15 and 30, which feels like an appropriate amount
fscanf(summary,"%d %lf %lf %lf %lf\n",&numberOfTracks, &vx,&vy,&vz,&B);
i++;
j+=numberOfTracks; // saves total number of tracks until it gets to the event we want - including that event
} while ( (numberOfTracks < 15 || numberOfTracks > 30) && i<100 ); //if it loops 100 times and does not find an appropriate number of tracks, stop looping
if(numberOfTracks < 15 || numberOfTracks > 30) { // if loop has ended and number of tracks is NOT between 15 and 30, finish business
printf("There are no events so that 15 < number of tracks < 30 \n\n");
return;
}
} while (numberOfTracks == 0);
vx = 0.01 * vx; // Vertexes are now in meters
vy = 0.01 * vy;
vz = 0.01 * vz;
B = 0.1 * B; // Magnetic field is now in Tesla
double m[numberOfTracks], q[numberOfTracks], px[numberOfTracks], py[numberOfTracks], pz[numberOfTracks];
for (i=0;i<numberOfTracks;i++) { // Reads only data regarding whichever event was the first to have a certain number of tracks
fscanf(detail,"%lf %lf %lf %lf %lf\n",&m[i],&q[i],&px[i],&py[i],&pz[i]);
double m[numberOfTracks], q[numberOfTracks], px[numberOfTracks], py[numberOfTracks], pz[numberOfTracks], Vx[numberOfTracks], Vy[numberOfTracks], Vz[numberOfTracks], lorentz[numberOfTracks];
j-=numberOfTracks; // now j equals the number of tracks before the selected event
for(i=0;i<j;i++) {
fscanf(detail,"%lf %lf %lf %lf %lf %lf %lf %lf \n", &skip[0], &skip[1], &skip[2], &skip[3], &skip[4], &skip[5], &skip[6], &skip[7]); // skips lines in detail file that contain information regarding other events
}
printf("Number of tracks: %d \nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n", numberOfTracks, vx, vy, vz, B);
for (i=0;i<numberOfTracks;i++) { // Reads only data regarding whichever event was the first to have an appropriate number of tracks
//Verify conservation of momentum:
fscanf(detail,"%lf %lf %lf %lf %lf %lf %lf %lf \n",&m[i],&q[i],&px[i],&py[i],&pz[i],&Vx[i],&Vy[i],&Vz[i]);
double pxSum=0.0, pySum=0.0, pzSum=0.0;
lorentz[i] = 1 / sqrt( 1 - Vx[i]*Vx[i]+Vy[i]*Vy[i]+Vz[i]*Vz[i] ); // Lorentz factor [SOMETHING FUNNY HERE...]
}
for(i=0;i<numberOfTracks;i++) {
printf("\nNumber of tracks: %d \n\nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n\n", numberOfTracks, vx, vy, vz, B);
pxSum += px[i];
pySum += py[i];
pzSum += pz[i];
//test
for(i=0;i<24;i++) {
printf("\nlorentz[%d] = %lf\n\n", i,lorentz[i]);
}
printf("Sum of X momentum: %lf \nSum of Y momentum: %lf \nSum of Z momentum: %lf \n\n",pxSum,pySum,pzSum);
// The following loops fill the position array for each track and for many different instants of time
......@@ -81,10 +106,12 @@ void main() {
Rx = px[i] / ( q[i] * B ) * 3.335641; // unit conversion to meters
Ry = py[i] / ( q[i] * B ) * 3.335641; // unit conversion to meters
w = ( q[i] * B ) / m[i]; // Angular frequency
w = ( q[i] * B ) / ( lorentz[i] * m[i] ); // Angular frequency
j = 0;
dt = fabs(1/w) / 1000.0; // time step is proportional to (1/w) so it takes plenty of steps
radius = 0;
//dt = fabs(1/w) / 1000.0; // time step is proportional to (1/w) so it takes plenty of steps
dt = 0.0001;
fprintf(positions,"\n"); // skip a line, for organization
......@@ -99,7 +126,7 @@ void main() {
gamaX[1] = Rx * ( cos(w*t) - 1 ); // 1 = y coordinate
// "particle with Y momentum only:
// "particle" with Y momentum only:
gamaY[1] = Ry * sin(w*t); // 0 = x coordinate
gamaY[0] = -Ry * ( cos(w*t) - 1 ); // 1 = y coordinate
......@@ -108,14 +135,15 @@ void main() {
// Main position vector:
double gama[3] = { gamaX[0]+gamaY[0] + vx, gamaX[1]+gamaY[1] + vy, (( pz[i] / m[i] ) * t) + vz }; // vertexes have already been converted from cm to m
fprintf(positions, "%e %e %e %e\n", t, gama[0], gama[1], gama[2]); // Print using scientific notation (very small numbers)
fprintf(positions, "%d %e %e %e %e\n", i, t, gama[0], gama[1], gama[2]); // Print using scientific notation (very small numbers)
radius = sqrt(gama[0]*gama[0] + gama[1]*gama[1]); // distance of particle from Z axis ("radial distance")
prevRadius = radius;
radius = sqrt( pow(gama[0]-vx,2) + pow(gama[1]-vy,2) ); // distance of particle from vertex ("radial distance")
j++;
} while( radius < TPCRADIUS ); // if the instant of time equals half the period, stop writing! (w = 2*PI / period)
} while( radius < TPCRADIUS && radius >= prevRadius ); // if the instant of time equals half the period, stop writing! (w = 2*PI / period) OR stop writing if it has taken half a lap (previous radius is, then, greater than current radius)
}
......@@ -124,4 +152,4 @@ void main() {
} // Program works for one specific event.
// Remember to organize AliAnalysisTaskMyTask.cxx
// Remember to organize AliAnalysisTaskMyTask.cxx 11834
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment