sol-analitica-tpc.c 5.27 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
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
int main(int argc, char **argv) {
25
26
27


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

	FILE *detail, *summary;

33
34
35
36
37
38
39
40
41
42
43
44
	// 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);
	}
45

46
47
48
	i=0;
	j=0;

49

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

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

62
		return 1;
63

64
65
66

	}

67
68
69
70
71
72
73
74
75


	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



76
77
78
79
80
81
82
	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
83
84
85

	}

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

88
		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]);
89

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
90
		lorentz[i] = 1 / sqrt( 1 - (Vx[i]*Vx[i]+Vy[i]*Vy[i]+Vz[i]*Vz[i]) ); // Lorentz factor
91

92
93
	}
	
94

95
	printf("\nNumber of tracks: %d \n\nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n\n", numberOfTracks, vx, vy, vz, B);
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

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
109
        FILE *positions, *histogram;
110
111

        positions = fopen("positions.txt","w");
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
112
	histogram = fopen("histogram.txt","w");
113
114
115

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

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
116
117
118
119
120
		// First, print histogram.txt information

		fprintf(histogram,"%d %lf\n",i,m[i]);	


121
122
123
124
125
		// 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

126
		w = ( q[i] * B ) / ( lorentz[i] * m[i] ); // Angular frequency
127
128

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

		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


146
		    // "particle" with Y momentum only:
147
148
149
150
151
152
153
154

		    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

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

157
158
		    prevRadius = radius;
		    radius = sqrt( pow(gama[0]-vx,2) + pow(gama[1]-vy,2) ); // distance of particle from vertex ("radial distance")
159
160
161
162

		    j++;


Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
163
		} 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 vertex distance is, then, greater than current vertex distance)
164
165
166
167
168
169

	}

	fclose(positions);


170
	return 0;
171
172
} // Program works for one specific event.

173
// Remember to organize AliAnalysisTaskMyTask.cxx 11834