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

5
#define PI 3.141592
6
#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

// 
// REMEMBER: Vx is for velocity; vx is for vertex.
19
20
21
22
//

#define ESD_SUMARIO "esd-summary.dat"
#define ESD_DETALHE "esd-detail.dat"
23

24
25
26
void main() {

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

	FILE *detail, *summary;

32
33
34
35
36
37
38
39
40
41
42
43
	// filename must be "esd-summary.dat"
	summary = fopen(ESD_SUMARIO,"r");
	if (summary == NULL){
	  fprintf(stderr, "O arquivo %s não pode ser aberto.\n", ESD_SUMARIO);
	  exit(1);
	}
	// filename must be "esd-detail.dat"
	detail = fopen(ESD_DETALHE,"r");
	if (detail == NULL){
	  fprintf(stderr, "O arquivo %s não pode ser aberto.\n", ESD_DETALHE);
	  exit(1);
	}
44

45
46
47
48
	i=0;
	j=0;

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

		fscanf(summary,"%d %lf %lf %lf %lf\n",&numberOfTracks, &vx,&vy,&vz,&B);
51
52
53
54
55
56
57
58
59
60
61
62
		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;

	}

63
64
65
66
67
68
69
70
71


	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



72
73
74
75
76
77
78
	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
79
80
81

	}

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

84
		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]);
85

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

88
89
	}
	
90

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


94
95
96
	//test
	for(i=0;i<24;i++) {
		printf("\nlorentz[%d] = %lf\n\n", i,lorentz[i]);
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
	}



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

122
		w = ( q[i] * B ) / ( lorentz[i] * m[i] ); // Angular frequency
123
124

		j = 0;
125
126
127
		radius = 0;
		//dt = fabs(1/w) / 1000.0; // time step is proportional to (1/w) so it takes plenty of steps
		dt = 0.0001;
128
129
130
131
132
133
134
135
136
137
138
139
140
141

		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


142
		    // "particle" with Y momentum only:
143
144
145
146
147
148
149
150

		    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

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

153
154
		    prevRadius = radius;
		    radius = sqrt( pow(gama[0]-vx,2) + pow(gama[1]-vy,2) ); // distance of particle from vertex ("radial distance")
155
156
157
158

		    j++;


159
		} 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)
160
161
162
163
164
165
166
167

	}

	fclose(positions);


} // Program works for one specific event.

168
// Remember to organize AliAnalysisTaskMyTask.cxx 11834