sol-analitica-tpc.c 5 KB
Newer Older
1
2
3
4
5
6
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.141592
#define TPCRADIUS 2.466
7
#define c 3.0E8
8
9
10
11

// Reads file with data in the following order:
//
// number of tracks, vx, vy, vz, magnetic field => esd-summary
12
// mass, charge, px, py, pz, Vx, Vy, Vz => esd-detail
13
14
15
//
// For now, we are considering classical mechanics, so p = mv (without relativistic considerations)

16
17
18
19
20

// 
// REMEMBER: Vx is for velocity; vx is for vertex.
// 

21
22
23
void main() {

	int i, j, numberOfTracks;
24
	double vx, vy, vz, w, Rx, Ry, Rz, B, radius, prevRadius, t, dt, skip[8];
25
26
27
28
29
30
31
	double gamaX[2], gamaY[2];

	FILE *detail, *summary;

	summary = fopen("esd-summary.dat","r"); // filename must be "esd-summary.dat"
	detail = fopen("esd-detail.dat","r"); // filename must be "esd-detail.dat"

32
33
34
35
	i=0;
	j=0;

	do { // Reads summary file until number of tracks is between 15 and 30, which feels like an appropriate amount
36
37

		fscanf(summary,"%d %lf %lf %lf %lf\n",&numberOfTracks, &vx,&vy,&vz,&B);
38
39
40
41
42
43
44
45
46
47
48
49
		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;

	}

50
51
52
53
54
55
56
57
58


	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



59
60
61
62
63
64
65
	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
66
67
68

	}

69
	for (i=0;i<numberOfTracks;i++) { // Reads only data regarding whichever event was the first to have an appropriate number of tracks
70

71
		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]);
72

73
		lorentz[i] = 1 / sqrt( 1 - Vx[i]*Vx[i]+Vy[i]*Vy[i]+Vz[i]*Vz[i] ); // Lorentz factor [SOMETHING FUNNY HERE...]
74

75
76
	}
	
77

78
	printf("\nNumber of tracks: %d \n\nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n\n", numberOfTracks, vx, vy, vz, B);
79
80


81
82
83
	//test
	for(i=0;i<24;i++) {
		printf("\nlorentz[%d] = %lf\n\n", i,lorentz[i]);
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
	}



	// The following loops fill the position array for each track and for many different instants of time
	// gamaX is for x and y coordinates of a particle with only px momentum
	// gamaY is for x and y coordinates of a particle with only py momentum
	//
	// gama ( particle position vector ):
	//       X coordinate = x coordinate of gamaX + x coordinate of gamaY + primary x vertex
	//       Y coordinate = y coordinate of gamaX + y coordinate of gamaY + primary y vertex
	//       X coordinate = z coordinate of particle (straight line) + primary z vertex

        FILE *positions;

        positions = fopen("positions.txt","w");

	for(i=0;i<numberOfTracks;i++) { // Loops in each track

	
		// Here we consider the momentum to be m * v
		// If momentum is zero, R is zero!!
		Rx = px[i] / ( q[i] * B ) * 3.335641; // unit conversion to meters
		Ry = py[i] / ( q[i] * B ) * 3.335641; // unit conversion to meters

109
		w = ( q[i] * B ) / ( lorentz[i] * m[i] ); // Angular frequency
110
111

		j = 0;
112
113
114
		radius = 0;
		//dt = fabs(1/w) / 1000.0; // time step is proportional to (1/w) so it takes plenty of steps
		dt = 0.0001;
115
116
117
118
119
120
121
122
123
124
125
126
127
128

		fprintf(positions,"\n"); // skip a line, for organization

		do { // Loops for every instant of time 

		    t = (double)(j * dt);


		    // "particle" with X momentum only:

		    gamaX[0] = Rx * sin(w*t); 			// 0 = x coordinate
		    gamaX[1] = Rx * ( cos(w*t) - 1 ); 		// 1 = y coordinate


129
		    // "particle" with Y momentum only:
130
131
132
133
134
135
136
137

		    gamaY[1] = Ry * sin(w*t); 			// 0 = x coordinate
		    gamaY[0] = -Ry * ( cos(w*t) - 1 ); 		// 1 = y coordinate


		    // 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

138
		    fprintf(positions, "%d %e %e %e %e\n", i, t, gama[0], gama[1], gama[2]); // Print using scientific notation (very small numbers)
139

140
141
		    prevRadius = radius;
		    radius = sqrt( pow(gama[0]-vx,2) + pow(gama[1]-vy,2) ); // distance of particle from vertex ("radial distance")
142
143
144
145

		    j++;


146
		} 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)
147
148
149
150
151
152
153
154

	}

	fclose(positions);


} // Program works for one specific event.

155
// Remember to organize AliAnalysisTaskMyTask.cxx 11834