Creat repository with both standard example and 'sol-analitica'

AliAnalysisTaskMyTask.cxx

@@ -28,8 +28,9 @@
 #include "AliESDInputHandler.h"
 #include "AliAnalysisTaskMyTask.h"
 #include "AliESDtrack.h"
-#include "Riostream.h"
 #include "AliESDVertex.h"
+#include "AliVertex.h"
+#include "Riostream.h"
 
 
 Int_t esd_event_id = 0; // global variable to store unique event id
@@ -110,44 +111,64 @@ void AliAnalysisTaskMyTask::UserExec(Option_t *)
 
     Int_t Event=0;
     Int_t TrigEvent=0;
-    ofstream summaryfile, eventsdetail;
-    summaryfile.open ("events_summary.txt",std::ofstream::app);
-    eventsdetail.open ("events_detail.txt",std::ofstream::app);
-
+    ofstream summary, detail;
+    summary.open ("esd-summary.dat",std::ofstream::app);
+    detail.open ("esd-detail.dat",std::ofstream::app);
 
 
     fESD = dynamic_cast<AliESDEvent*>(InputEvent());    // get an event (called fESD) from the input file
 
-                                                        // there's another event format (ESD) which works in a similar wya
+                                                        // there's another event format (ESD) which works in a similar way
 
-                                                  // but is more cpu/memory unfriendly. for now, we'll stick with aod's
+                                                 
     if(!fESD) return;                                   // if the pointer to the event is empty (getting it failed) skip this event
         // example part: i'll show how to loop over the tracks in an event
         // and extract some information from them which we'll store in a histogram
 
     Int_t iTracks(fESD->GetNumberOfTracks());           // see how many tracks there are in the event
 
-    Double_t vX = fESD->GetPrimaryVertex()->GetX();
-    Double_t vY = fESD->GetPrimaryVertex()->GetY();
-    Double_t vZ = fESD->GetPrimaryVertex()->GetZ();
+    Double_t Vx = fESD->GetPrimaryVertex()->GetX();	// gets vertexes from individual events
+    Double_t Vy = fESD->GetPrimaryVertex()->GetY();
+    Double_t Vz = fESD->GetPrimaryVertex()->GetZ();
+    Double_t MagneticField = fESD->GetMagneticField();
+
+    summary << iTracks << " " << Vx << " " << Vy << " " << Vz << " " << MagneticField << endl;
+
+/*
 
-    // Add Event details to summaryfile and eventsdetail
-    summaryfile << " Event " << esd_event_id << ": iTracks = "  << iTracks << "  Vertex: X = " << vX << " Y = " << vY << " Z = " << vZ << endl;
-    eventsdetail << "********** Event " << esd_event_id << ": iTracks = "  << iTracks << " **********" << endl;
+Assumed Units: Mass (GeV/c^2)[CONFIRMED] || Energy (GeV) || Momentum (GeV/c) || Charge (* 1.6*10^-19 C)
+
+*/
+
+    for(Int_t i(0); i < iTracks; i++) {                 // loop over all these tracks
 
-    for(Int_t i(0); i < iTracks; i++) {                 // loop ove rall these tracks
         AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(i));         // get a track (type AliESDtrack) from the event
+
         if(!track) continue;                            // if we failed, skip this track
-        // Add track number and corresponting pt to eventsdetail txt file
-        eventsdetail << esd_event_id << " " << i << " " << track->Pt() << endl;
+	
+	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 Pz = track->Pz();
+	Double_t Charge = track->Charge();
+
+        // Add MASS, CHARGE, MOMENTUM (x, y, z) to esd-detail.dat file
+	detail << fixed << Mass << " " << Charge << " ";
+	detail << fixed << Px << " " << Py << " " << Pz << endl;
+
+	//'fixed' fixes the number of decimal places so numbers are vertically aligned
+
         fHistPt->Fill(track->Pt());                     // plot the pt value of the track in a histogram
     }
 
     Event++;
     esd_event_id++; // Increment global esd_event_id
     fHistEvents->Fill(Event);
-    summaryfile.close();
-    eventsdetail.close();
+    summary.close();
+    detail.close();
 
                                                        // continue until all the tracks are processed
     PostData(1, fOutputList);                           // stream the results the analysis of this event to

notes.txt

+fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+
+cria um ponteiro para um objeto do tipo AliESDEvent, que é uma classe filha do tipo AliVEvent.
+Sabemos disso pois é definida assim:
+
+class AliESDEvent : public AliVEvent {
+
+(...)
+
+}
+
+
+A classe AliESDEvent tem um método chamado GetPrimaryVertex(), que retorna um ponteiro para objeto do tipo AliESDVertex.
+Sabemos disso pois a função GetPrimaryVertex() está definida assim, no corpo de AliESDEvent:
+
+const AliESDVertex *GetPrimaryVertex() const;
+
+
+A classe AliESDVertex, por sua vez, é filha do tipo AliVertex:
+
+class AliESDVertex : public AliVertex {
+
+(...)
+
+}
+
+
+A classe AliVertex tem métodos para os vértices de cada evento: GetX(), GetY() e GetZ().
+Portanto, faz sentido escrever (como em AliAnalysisTaskMyTask.cxx):
+
+Double_t Vx = fESD->GetPrimaryVertex()->GetX();
+
+PESQUISAR E ADICIONAR NA WIKI: https://github.com/alisw/AliRoot/tree/master/STEER/ESD
+*CLUSTERS* + *ZDC*

runAnalysis.C

 void runAnalysis()
 {
     // Erase output txt files
-    ofstream summaryfile, eventsdetail;
-    summaryfile.open ("events_summary.txt");
-    summaryfile << " " << endl;
-    summaryfile.close();
+    ofstream summary, detail;
+    summary.open ("esd-summary.dat");
+    summary << " " << endl;
+    summary.close();
 
-    eventsdetail.open ("events_detail.txt");
-    eventsdetail << " " << endl;
-    eventsdetail.close();
+    detail.open ("esd-detail.dat");
+    detail << " " << endl;
+    detail.close();
 
 
 
@@ -36,7 +36,8 @@ void runAnalysis()
     // if you want to run locally, we need to define some input
     TChain* chain = new TChain("esdTree");
     // add a few files to the chain (change this so that your local files are added)
-    chain->Add("AliESDs.root");
+    chain->Add("AliESDs.root"); // Breno put it on the same directory that was cloned from Pezzi's 					// repository: AliESD_Example
+
     //chain->Add("../root_files/AliAOD.Muons2.root");
     //chain->Add("../root_files/AliAOD.Muons3.root");
     //chain->Add("../root_files/AliAOD.Muons4.root");

sol-analitica/AddMyTask.C

+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* AddMyTask
+ *
+ * empty task which can serve as a starting point for building an analysis
+ * as an example, one histogram is filled
+ */
+ // Author: Redmer A. Bertens, Utrecht University, 2012
+ 
+class AliAnalysisDataContainer;
+
+AliAnalysisTaskMyTask* AddMyTask(TString name = "name")
+{
+    // get the manager via the static access member. since it's static, you don't need
+    // an instance of the class to call the function
+    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+    if (!mgr) {
+        return 0x0;
+    }
+    // get the input event handler, again via a static method.
+    // this handler is part of the managing system and feeds events
+    // to your task
+    if (!mgr->GetInputEventHandler()) {
+        return 0x0;
+    }
+    // by default, a file is open for writing. here, we get the filename
+    TString fileName = AliAnalysisManager::GetCommonFileName();
+    fileName += ":MyTask";      // create a subfolder in the file
+    // now we create an instance of your task
+    AliAnalysisTaskMyTask* task = new AliAnalysisTaskMyTask(name.Data());
+    if(!task) return 0x0;
+    // add your task to the manager
+    mgr->AddTask(task);
+    // your task needs input: here we connect the manager to your task
+    mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+    // same for the output
+    mgr->ConnectOutput(task,1,mgr->CreateContainer("MyOutputContainer", TList::Class(), AliAnalysisManager::kOutputContainer, fileName.Data()));
+    // in the end, this macro returns a pointer to your task. this will be convenient later on
+    // when you will run your analysis in an analysis train on grid
+    return task;
+}

sol-analitica/AliAnalysisTaskMyTask.cxx

+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* AliAnaysisTaskMyTask
+ *
+ * empty task which can serve as a starting point for building an analysis
+ * as an example, one histogram is filled
+ */
+
+#include "TChain.h"
+#include "TH1F.h"
+#include "TList.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliAnalysisTaskMyTask.h"
+#include "AliESDtrack.h"
+#include "AliESDVertex.h"
+#include "AliVertex.h"
+#include "Riostream.h"
+
+
+Int_t esd_event_id = 0; // global variable to store unique event id
+
+class AliAnalysisTaskMyTask;    // your analysis class
+
+using namespace std;            // std namespace: so you can do things like 'cout'
+
+ClassImp(AliAnalysisTaskMyTask) // classimp: necessary for root
+
+AliAnalysisTaskMyTask::AliAnalysisTaskMyTask() : AliAnalysisTaskSE(),
+    fESD(0), fOutputList(0), fHistPt(0),    fHistEvents(0)
+{
+    // default constructor, don't allocate memory here!
+    // this is used by root for IO purposes, it needs to remain empty
+}
+//_____________________________________________________________________________
+AliAnalysisTaskMyTask::AliAnalysisTaskMyTask(const char* name) : AliAnalysisTaskSE(name),
+    fESD(0), fOutputList(0), fHistPt(0),  fHistEvents(0)
+{
+    // constructor
+    DefineInput(0, TChain::Class());    // define the input of the analysis: in this case we take a 'chain' of events
+                                        // this chain is created by the analysis manager, so no need to worry about it,
+                                        // it does its work automatically
+    DefineOutput(1, TList::Class());    // define the ouptut of the analysis: in this case it's a list of histograms
+                                        // you can add more output objects by calling DefineOutput(2, classname::Class())
+                                        // if you add more output objects, make sure to call PostData for all of them, and to
+                                        // make changes to your AddTask macro!
+}
+//_____________________________________________________________________________
+AliAnalysisTaskMyTask::~AliAnalysisTaskMyTask()
+{
+    // destructor
+    if(fOutputList) {
+        delete fOutputList;     // at the end of your task, it is deleted from memory by calling this function
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskMyTask::UserCreateOutputObjects()
+{
+    // create output objects
+    //
+    // this function is called ONCE at the start of your analysis (RUNTIME)
+    // here you ceate the histograms that you want to use
+    //
+    // the histograms are in this case added to a tlist, this list is in the end saved
+    // to an output file
+    //
+    fOutputList = new TList();          // this is a list which will contain all of your histograms
+                                        // at the end of the analysis, the contents of this list are written
+                                        // to the output file
+    fOutputList->SetOwner(kTRUE);       // memory stuff: the list is owner of all objects it contains and will delete them
+                                        // if requested (dont worry about this now)
+
+    // example of a histogram
+    fHistPt = new TH1F("fHistPt", "fHistPt", 100, 0, 10);       // create your histogra
+    fOutputList->Add(fHistPt);          // don't forget to add it to the list! the list will be written to file, so if you want
+                                        // your histogram in the output file, add it to the list!
+
+
+    // Histograms for dimuons
+    fHistEvents = new TH1F("fHistEvents","fHistEvents;N_{events}",100,0.,10000.);
+
+    fOutputList->Add(fHistEvents);
+
+    PostData(1, fOutputList);           // postdata will notify the analysis manager of changes / updates to the
+                                        // fOutputList object. the manager will in the end take care of writing your output to file
+                                        // so it needs to know what's in the output
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskMyTask::UserExec(Option_t *)
+{
+    // user exec
+    // this function is called once for each event
+    // the manager will take care of reading the events from file, and with the static function InputEvent() you
+    // have access to the current event.
+    // once you return from the UserExec function, the manager will retrieve the next event from the chain
+
+    Int_t Event=0;
+    Int_t TrigEvent=0;
+    ofstream summary, detail;
+    summary.open ("esd-summary.dat",std::ofstream::app);
+    detail.open ("esd-detail.dat",std::ofstream::app);
+
+
+    fESD = dynamic_cast<AliESDEvent*>(InputEvent());    // get an event (called fESD) from the input file
+
+                                                        // there's another event format (ESD) which works in a similar way
+
+                                                 
+    if(!fESD) return;                                   // if the pointer to the event is empty (getting it failed) skip this event
+        // example part: i'll show how to loop over the tracks in an event
+        // and extract some information from them which we'll store in a histogram
+
+    Int_t iTracks(fESD->GetNumberOfTracks());           // see how many tracks there are in the event
+
+    Double_t Vx = fESD->GetPrimaryVertex()->GetX();	// gets vertexes from individual events
+    Double_t Vy = fESD->GetPrimaryVertex()->GetY();
+    Double_t Vz = fESD->GetPrimaryVertex()->GetZ();
+    Double_t MagneticField = fESD->GetMagneticField();
+
+    summary << iTracks << " " << Vx << " " << Vy << " " << Vz << " " << MagneticField << endl;
+
+/*
+
+Assumed Units: Mass (GeV/c^2)[CONFIRMED] || Energy (GeV) || Momentum (GeV/c) || Charge (* 1.6*10^-19 C)
+
+*/
+
+    for(Int_t i(0); i < iTracks; i++) {                 // loop over all these tracks
+
+        AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(i));         // get a track (type AliESDtrack) from the event
+
+        if(!track) continue;                            // if we failed, skip this track
+	
+	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 Pz = track->Pz();
+	Double_t Charge = track->Charge();
+
+        // Add MASS, CHARGE, MOMENTUM (x, y, z) to esd-detail.dat file
+	detail << fixed << Mass << " " << Charge << " ";
+	detail << fixed << Px << " " << Py << " " << Pz << endl;
+
+	//'fixed' fixes the number of decimal places so numbers are vertically aligned
+
+        fHistPt->Fill(track->Pt());                     // plot the pt value of the track in a histogram
+    }
+
+    Event++;
+    esd_event_id++; // Increment global esd_event_id
+    fHistEvents->Fill(Event);
+    summary.close();
+    detail.close();
+
+                                                       // continue until all the tracks are processed
+    PostData(1, fOutputList);                           // stream the results the analysis of this event to
+                                                        // the output manager which will take care of writing
+                                                        // it to a file
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskMyTask::Terminate(Option_t *)
+{
+    // terminate
+    // called at the END of the analysis (when all events are processed)
+}
+//_____________________________________________________________________________

sol-analitica/AliAnalysisTaskMyTask.h

+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. */
+/* See cxx source for full Copyright notice */
+/* $Id$ */
+
+#ifndef AliAnalysisTaskMyTask_H
+#define AliAnalysisTaskMyTask_H
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskMyTask : public AliAnalysisTaskSE
+{
+    public:
+                                AliAnalysisTaskMyTask();
+                                AliAnalysisTaskMyTask(const char *name);
+        virtual                 ~AliAnalysisTaskMyTask();
+
+        virtual void            UserCreateOutputObjects();
+        virtual void            UserExec(Option_t* option);
+        virtual void            Terminate(Option_t* option);
+
+    private:
+        AliESDEvent*            fESD;           //! input event
+        TList*                  fOutputList;    //! output list
+        TH1F*                   fHistPt;        //! dummy histogram
+        TH1F*                   fHistEvents; // Number of events
+
+        AliAnalysisTaskMyTask(const AliAnalysisTaskMyTask&); // not implemented
+        AliAnalysisTaskMyTask& operator=(const AliAnalysisTaskMyTask&); // not implemented
+
+        ClassDef(AliAnalysisTaskMyTask, 1);
+};
+
+#endif

sol-analitica/README.md

+# Sample macros for processing AliESDs from CERN Open Data Portal
+
+## Requirements
+
+Execute with aliroot5. Place `AliESDs.root` files in reach of the `chain->Add("AliESDs.root");` line on `runAnalysis.C` and execute.
+
+
+## Instructions
+
+1) Install aliBuild. Follow instructions on https://alice-doc.github.io/alice-analysis-tutorial/building/custom.html
+
+2) Initialize AliPhysics
+
+```bash
+mkdir -p ~/alice
+cd ~/alice
+aliBuild init AliPhysics@master
+```
+3) Verify dependencies (Optional)
+
+```bash
+$ aliDoctor AliPhysics
+```
+4) Build AliPhysics with aliroot5 (this may take a long time)
+```bash
+aliBuild build AliPhysics --defaults user -z aliroot5
+```
+5) Enter AliPhysics environment
+```bash
+alienv enter AliPhysics/latest-aliroot5-user
+```
+6) Run the macro
+
+```bash
+aliroot runAnalysis.C
+```
+Results will be saved on two text files: `events_summary.txt` and `events_detail.txt`.
+
+## Credits
+
+This example is based on the ALICE analysis tutorial by Redmer A. Bertens.

sol-analitica/dumpTreeTotxt.C

+// Exemplo obtido de https://stackoverflow.com/questions/28970124/cern-root-exporting-data-to-plain-text
+#include <iostream>
+#include "TFile.h"
+#include "TTree.h"
+#include <fstream>
+using namespace std;
+
+void dumpTreeTotxt(){
+  TFile *f=new TFile("AliESDs.root"); // opens the root file
+  TTree *tr=(TTree*)f->Get("esdTree"); // creates the TTree object
+  //tr->Scan(); // prints the content on the screen
+
+  float a,b,c; // create variables of the same type as the branches you want to access
+
+  tr->SetBranchAddress("AliESDRun",&a); // for all the TTree branches you need this
+//  tr->SetBranchAddress("AliESDHeader",&b);
+  //tr->SetBranchAddress("nserr",&c);
+
+  ofstream myfile;
+  myfile.open ("example.txt");
+  //myfile << "TS ns nserr\n";
+
+  for (int i=0;i<tr->GetEntries();i++){
+    // loop over the tree
+
+    cout   << "Event: " << i <<  endl; //print to the screen
+    myfile << "Event: " << i <<  endl; //write to file
+    tr->GetEntry(i);
+  }
+  myfile.close();
+}

sol-analitica/runAnalysis.C

+void runAnalysis()
+{
+    // Erase output txt files
+    ofstream summary, detail;
+    summary.open ("esd-summary.dat");
+    summary << " " << endl;
+    summary.close();
+
+    detail.open ("esd-detail.dat");
+    detail << " " << endl;
+    detail.close();
+
+
+
+    // since we will compile a class, tell root where to look for headers
+    gROOT->ProcessLine(".include $ROOTSYS/include");
+    gROOT->ProcessLine(".include $ALICE_ROOT/include");
+
+    // create the analysis manager
+    AliAnalysisManager *mgr = new AliAnalysisManager("testAnalysis");
+    AliESDInputHandler *esdH = new AliESDInputHandler();
+    mgr->SetInputEventHandler(esdH);
+
+    // compile the class (locally)
+    gROOT->LoadMacro("AliAnalysisTaskMyTask.cxx++g");
+    // load the addtask macro
+    gROOT->LoadMacro("AddMyTask.C");
+    // create an instance of your analysis task
+    AliAnalysisTaskMyTask *task = AddMyTask();
+
+    if(!mgr->InitAnalysis()) return;
+    mgr->SetDebugLevel(2);
+    mgr->PrintStatus();
+    mgr->SetUseProgressBar(1, 25);
+
+    // if you want to run locally, we need to define some input
+    TChain* chain = new TChain("esdTree");
+    // add a few files to the chain (change this so that your local files are added)
+    chain->Add("AliESDs.root"); // Breno put it on the same directory that was cloned from Pezzi's 					// repository: AliESD_Example
+
+    //chain->Add("../root_files/AliAOD.Muons2.root");
+    //chain->Add("../root_files/AliAOD.Muons3.root");
+    //chain->Add("../root_files/AliAOD.Muons4.root");
+
+    // start the analysis locally, reading the events from the tchain
+    mgr->StartAnalysis("local", chain);
+
+
+    exit();
+}

sol-analitica/sol-analitica-tpc.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define PI 3.141592
+#define TPCRADIUS 2.466
+
+// 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
+//
+// For now, we are considering classical mechanics, so p = mv (without relativistic considerations)
+
+void main() {
+
+	int i, j, numberOfTracks;
+	double vx, vy, vz, w, Rx, Ry, Rz, B, radius, t, dt;
+	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"
+
+	do { // Reads summary file until number of tracks is different than zero (there are tracks in the event)
+
+		fscanf(summary,"%d %lf %lf %lf %lf\n",&numberOfTracks, &vx,&vy,&vz,&B);
+
+	} 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]);
+
+	}
+
+	printf("Number of tracks: %d \nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n", numberOfTracks, vx, vy, vz, B);
+
+	//Verify conservation of momentum:
+
+	double pxSum=0.0, pySum=0.0, pzSum=0.0;
+
+