Commit 8eb5abfd authored by Breno Rilho Lemos's avatar Breno Rilho Lemos 💬
Browse files

Commit before pulling lucas

parent d0c379b0
......@@ -132,6 +132,7 @@ void AliAnalysisTaskMyTask::UserExec(Option_t *)
Double_t Vz = fESD->GetPrimaryVertex()->GetZ();
Double_t MagneticField = fESD->GetMagneticField();
// Add VERTEZES (x, y, z) and magnetic field to esd-summary.dat file
summary << iTracks << " " << Vx << " " << Vy << " " << Vz << " " << MagneticField << endl;
/*
......@@ -148,16 +149,22 @@ Assumed Units: Mass (GeV/c^2)[CONFIRMED] || Energy (GeV) || Momentum (GeV/c) ||
Double_t Mass = track->M(); // returns the pion mass, if the particle can't be identified properly
Double_t Energy = track->E(); // Returns the energy of the particle given its assumed mass, but assumes the pion mass if the particle can't be identified properly.
Double_t Px = track->Px();
Double_t Py = track->Py();
Double_t Pt = track->Pt();
Double_t PtCal = sqrt(Px*Px+Py*Py); // here, we calculate transversal momentum to make sure of its meaning (then compare it to Pt)
Double_t Pt = track->Pt(); // transversal momentum, in case we need it
Double_t Pz = track->Pz();
Double_t Vel = 1 / sqrt( Mass*Mass / (Pt*Pt + Pz*Pz) + 1 );
Double_t Vx = Px * sqrt( 1 - Vel*Vel ) / Mass;
Double_t Vy = Py * sqrt( 1 - Vel*Vel ) / Mass;
Double_t Vz = Pz * sqrt( 1 - Vel*Vel ) / Mass;
Double_t Charge = track->Charge();
// Add MASS, CHARGE, MOMENTUM (x, y, z) to esd-detail.dat file
// Add MASS, CHARGE and VELOCITIES (x, y, z) to esd-detail.dat file
detail << fixed << Mass << " " << Charge << " ";
detail << fixed << Px << " " << Py << " " << Pz << endl;
detail << fixed << Px << " " << Py << " " << Pz << " " << Vx << " " << Vy << " " << Vz << endl;
//'fixed' fixes the number of decimal places so numbers are vertically aligned
......
......@@ -4,18 +4,24 @@
#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