...
 
Commits (20)
......@@ -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
......@@ -41,14 +42,14 @@ using namespace std; // std namespace: so you can do things like 'cou
ClassImp(AliAnalysisTaskMyTask) // classimp: necessary for root
AliAnalysisTaskMyTask::AliAnalysisTaskMyTask() : AliAnalysisTaskSE(),
fESD(0), fOutputList(0), fHistPt(0), fHistEvents(0)
fESD(0), fOutputList(0), fHistPt(0), fHistEvents(0), fHistMass(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)
fESD(0), fOutputList(0), fHistPt(0), fHistEvents(0), fHistMass(0)
{
// constructor
DefineInput(0, TChain::Class()); // define the input of the analysis: in this case we take a 'chain' of events
......@@ -85,10 +86,16 @@ void AliAnalysisTaskMyTask::UserCreateOutputObjects()
// if requested (dont worry about this now)
// example of a histogram
fHistPt = new TH1F("fHistPt", "fHistPt", 100, 0, 10); // create your histogra
fHistPt = new TH1F("fHistPt", "fHistPt", 100, 0, 10); // create your histogram
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!
// my particle mass histogram
Double_t fHistMassEdges[12] = {0.0,0.0005,0.0405,0.08,0.12,0.13,0.17,0.48,0.52,0.92,0.96,1.0}; // 11 bins =>> has 11+1 = 12 edges
fHistMass = new TH1F("fHistMass","Particle Histogram;M_{particle}", 11, fHistMassEdges);
fOutputList->Add(fHistMass);
// Histograms for dimuons
fHistEvents = new TH1F("fHistEvents","fHistEvents;N_{events}",100,0.,10000.);
......@@ -118,7 +125,7 @@ void AliAnalysisTaskMyTask::UserExec(Option_t *)
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
......@@ -127,20 +134,52 @@ void AliAnalysisTaskMyTask::UserExec(Option_t *)
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();
// 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;
summaryfile << " Event " << esd_event_id << ": iTracks = " << iTracks << endl;
summaryfile << " Magnetic Field: " << MagneticField << endl;
summaryfile << " Vertex: " << endl << " X = " << Vx << endl << " Y = " << Vy << endl << " Z = " << Vz << endl << endl;
eventsdetail << "\n" << "********** Event " << esd_event_id << ": iTracks = " << iTracks << " **********" << endl << endl;
eventsdetail << "Event\tTrack\tMass\t\tEnergy\t\tPx\t\tPy\t\tPt(file)\tPt(calculado)\tPz\t\tCharge" << endl << 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 EVENT, TRACK NUMBER, MASS, ENERGY, corresponding MOMENTUM (x, y, transversal, z) and CHARGE to eventsdetail txt file
eventsdetail << esd_event_id << "\t" << i << "\t";
eventsdetail << fixed << Mass << "\t" << Energy << "\t";
eventsdetail << fixed << Px << "\t" << Py << " \t " << Pt << " \t " << PtCal << " \t " << Pz << "\t" << Charge << 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
fHistMass->Fill(Mass);
}
Event++;
......
......@@ -22,6 +22,7 @@ class AliAnalysisTaskMyTask : public AliAnalysisTaskSE
AliESDEvent* fESD; //! input event
TList* fOutputList; //! output list
TH1F* fHistPt; //! dummy histogram
TH1F* fHistMass; //!my particle mass histogram !! :D
TH1F* fHistEvents; // Number of events
AliAnalysisTaskMyTask(const AliAnalysisTaskMyTask&); // not implemented
......
File added
#+TITLE: Alice LabBook
#+AUTHOR: LMS
#+STARTUP: overview indent
#+TAGS: noexport(n) deprecated(d) Marcelo(m) Lucas(l)
#+EXPORT_SELECT_TAGS: export
#+EXPORT_EXCLUDE_TAGS: noexport
#+SEQ_TODO: TODO(t!) STARTED(s!) WAITING(w!) | DONE(d!) CANCELLED(c!) DEFERRED(f!)
This document is in English.
* Introduction
Most of the code in this Laboratory Notebook is written in R.
#+begin_src R :results output :session :exports both
suppressMessages(library(tidyverse))
sessionInfo()
#+end_src
#+RESULTS:
#+begin_example
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux bullseye/sid
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.3.7.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
[5] readr_1.3.1 tidyr_0.8.3 tibble_2.1.3 ggplot2_3.2.1
[9] tidyverse_1.2.1
loaded via a namespace (and not attached):
[1] Rcpp_1.0.2 cellranger_1.1.0 pillar_1.4.2 compiler_3.6.1
[5] tools_3.6.1 zeallot_0.1.0 jsonlite_1.6 lubridate_1.7.4
[9] gtable_0.3.0 nlme_3.1-141 lattice_0.20-38 pkgconfig_2.0.2
[13] rlang_0.4.0 cli_1.1.0 rstudioapi_0.10 haven_2.1.1
[17] withr_2.1.2 xml2_1.2.2 httr_1.4.1 generics_0.0.2
[21] vctrs_0.2.0 hms_0.5.1 grid_3.6.1 tidyselect_0.2.5
[25] glue_1.3.1 R6_2.4.0 readxl_1.3.1 modelr_0.1.5
[29] magrittr_1.5 backports_1.1.4 scales_1.0.0 rvest_0.3.4
[33] assertthat_0.2.1 colorspace_1.4-1 stringi_1.4.3 lazyeval_0.2.2
[37] munsell_0.5.0 broom_0.5.2 crayon_1.3.4
#+end_example
* 2019-10-11 Plot data created by sol-analitica
** Data :ATTACH:
:PROPERTIES:
:Attachments: positions.zip
:ID: f33a9852-afa9-4cda-b9e8-34b1f2dfe368
:END:
** Unzip and Read
#+begin_src R :results output :session :exports both
FILE <- "data/f3/3a9852-afa9-4cda-b9e8-34b1f2dfe368/positions.zip"
unzip(FILE)
read_delim("positions.txt", delim=" ", col_names=FALSE) %>%
rename(Track = X1,
Time = X2,
X = X3,
Y = X4,
Z = X5) %>%
print -> df.exp0
#+end_src
#+RESULTS:
#+begin_example
Parsed with column specification:
cols(
X1 = col_double(),
X2 = col_double(),
X3 = col_double(),
X4 = col_double(),
X5 = col_double()
)
# A tibble: 109,900 x 5
Track Time X Y Z
<dbl> <dbl> <dbl> <dbl> <dbl>
 1 0 0 -0.0000773 0.00170 0.195
 2 0 0.0001 0.0000994 0.00165 0.195
 3 0 0.0002 0.000276 0.00160 0.195
 4 0 0.000300 0.000453 0.00155 0.195
 5 0 0.0004 0.000630 0.00150 0.195
 6 0 0.0005 0.000806 0.00146 0.195
 7 0 0.000600 0.000983 0.00141 0.195
 8 0 0.0007 0.00116 0.00136 0.195
 9 0 0.0008 0.00134 0.00131 0.195
10 0 0.0009 0.00151 0.00126 0.195
# … with 109,890 more rows
#+end_example
** Plot
See here
https://stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2
on how to draw a circle.
#+begin_src R :results output graphics :file img/positions.png :exports both :width 600 :height 600 :session
circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
r = diameter / 2
tt <- seq(0,2*pi,length.out = npoints)
xx <- center[1] + r * cos(tt)
yy <- center[2] + r * sin(tt)
return(data.frame(x = xx, y = yy))
}
df.tpc.outer <- circleFun(c(0,0), 4.95, npoints = 100)
df.tpc.inner <- circleFun(c(0,0), 1.80, npoints = 100)
df.exp0 %>%
ggplot(aes(x=X, y=Y)) +
# Draw the tracks
geom_point(aes(color=as.factor(Track)), size=.1) +
# Draw approximate TPC
geom_path(data=df.tpc.outer, aes(x=x, y=y)) +
geom_path(data=df.tpc.inner, aes(x=x, y=y)) +
# Cosmetics
theme_bw(base_size=16) +
theme(legend.position="none",
axis.title=element_blank(),
axis.text=element_blank(),
axis.ticks=element_blank())
#+end_src
#+RESULTS:
[[file:img/positions.png]]
*
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*
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");
......
/**************************************************************************
* 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;
}
/**************************************************************************
* 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
Int_t EventID = 0; // is equal to esd-event-id until the selected event (one with an appropriate number of tracks) is reached
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), fHistMass(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), fHistMass(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!
// my mass histogram
Double_t fHistMassEdges[12] = {0.0,0.0005,0.0405,0.08,0.12,0.13,0.17,0.48,0.52,0.92,0.96,1.0}; // 11 bins =>> has 11+1 = 12 edges
fHistMass = new TH1F("fHistMass","Particle Histogram;M_{particle}", 11, fHistMassEdges);
fOutputList->Add(fHistMass);
// 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();
// Add VERTEZES (x, y, z) and magnetic field to esd-summary.dat file
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)
*/
if(EventID == esd_event_id) {
if(iTracks >= 15 && iTracks <= 30) {EventID = -1;}
else {EventID++;}
} else {EventID = -2;}
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(); // 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 and VELOCITIES (x, y, z) to esd-detail.dat file
detail << fixed << Mass << " " << Charge << " ";
detail << fixed << Px << " " << Py << " " << Pz << " " << Vx << " " << Vy << " " << Vz << endl;
//'fixed' fixes the number of decimal places so numbers are vertically aligned
fHistPt->Fill(Pt); // plot the pt value of the track in a histogram
if(EventID == -1) { // when we get to the selected event, fill Mass Histogram
fHistMass->Fill(Mass);
}
}
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)
}
//_____________________________________________________________________________
/* 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* fHistMass; //! my particle histogram!! :D
TH1F* fHistEvents; // Number of events
AliAnalysisTaskMyTask(const AliAnalysisTaskMyTask&); // not implemented
AliAnalysisTaskMyTask& operator=(const AliAnalysisTaskMyTask&); // not implemented
ClassDef(AliAnalysisTaskMyTask, 1);
};
#endif
# 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.
// 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();
}
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();
}
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#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, Vx, Vy, Vz => esd-detail
//
//
// REMEMBER: Vx is for velocity; vx is for vertex.
//
#define ESD_SUMARIO "esd-summary.dat"
#define ESD_DETALHE "esd-detail.dat"
int main(int argc, char **argv) {
int i, j, numberOfTracks;
double vx, vy, vz, w, Rx, Ry, Rz, B, radius, prevRadius, t, dt, skip[8];
double gamaX[2], gamaY[2];
FILE *detail, *summary;
// 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);
}
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 1;
}
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], 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
}
for (i=0;i<numberOfTracks;i++) { // Reads only data regarding whichever event was the first to have an appropriate number of tracks
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]);
lorentz[i] = 1 / sqrt( 1 - (Vx[i]*Vx[i]+Vy[i]*Vy[i]+Vz[i]*Vz[i]) ); // Lorentz factor
}
printf("\nNumber of tracks: %d \n\nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n\n", numberOfTracks, vx, vy, vz, B);
// 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, *histogram;
positions = fopen("positions.txt","w");
histogram = fopen("histogram.txt","w");
for(i=0;i<numberOfTracks;i++) { // Loops in each track
// First, print histogram.txt information
fprintf(histogram,"%d %lf\n",i,m[i]);
// 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
w = ( q[i] * B ) / ( lorentz[i] * m[i] ); // Angular frequency
j = 0;
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
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
// "particle" with Y momentum only:
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
fprintf(positions, "%d %e %e %e %e\n", i, t, gama[0], gama[1], gama[2]); // Print using scientific notation (very small numbers)
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 && 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)
}
fclose(positions);
return 0;
} // Program works for one specific event.
// Remember to organize AliAnalysisTaskMyTask.cxx 11834
#!/usr/bin/gnuplot -persist
#
#
# G N U P L O T
# Version 5.2 patchlevel 2 last modified 2017-11-01
#
# Copyright (C) 1986-1993, 1998, 2004, 2007-2017
# Thomas Williams, Colin Kelley and many others
#
# gnuplot home: http://www.gnuplot.info
# faq, bugs, etc: type "help FAQ"
# immediate help: type "help" (plot window: hit 'h')
# set terminal wxt 0 enhanced
# set output
unset clip points
set clip one
unset clip two
set errorbars front 1.000000
set border 31 front lt black linewidth 1.000 dashtype solid
set zdata
set ydata
set xdata
set y2data
set x2data
set boxwidth
set style fill empty border
set style rectangle back fc bgnd fillstyle solid 1.00 border lt -1
set style circle radius graph 0.02
set style ellipse size graph 0.05, 0.03 angle 0 units xy
set dummy x, y
set format x "% h"
set format y "% h"
set format x2 "% h"
set format y2 "% h"
set format z "% h"
set format cb "% h"
set format r "% h"
set ttics format "% h"
set timefmt "%d/%m/%y,%H:%M"
set angles radians
set tics back
set grid nopolar
set grid xtics nomxtics ytics nomytics noztics nomztics nortics nomrtics \
nox2tics nomx2tics noy2tics nomy2tics nocbtics nomcbtics
set grid layerdefault lt 0 linecolor 0 linewidth 0.500 dashtype solid, lt 0 linecolor 0 linewidth 0.500 dashtype solid
unset raxis
set theta counterclockwise right
set style parallel front lt black linewidth 2.000 dashtype solid
set key title "" center
set key fixed right top vertical Right noreverse enhanced autotitle nobox
set key noinvert samplen 4 spacing 1 width 0 height 0
set key maxcolumns 0 maxrows 0
set key noopaque
unset label
unset arrow
set style increment default
unset style line
unset style arrow
set style histogram clustered gap 2 title textcolor lt -1
unset object
set style textbox transparent margins 1.0, 1.0 border lt -1 linewidth 1.0
set offsets 0, 0, 0, 0
set pointsize 1
set pointintervalbox 1
set encoding default
unset polar
unset parametric
unset decimalsign
unset micro
unset minussign
set view 60, 30, 1, 1
set view azimuth 0
set rgbmax 255
set samples 100, 100
set isosamples 10, 10
set surface
unset contour
set cntrlabel format '%8.3g' font '' start 5 interval 20
set mapping cartesian
set datafile separator whitespace
unset hidden3d
set cntrparam order 4
set cntrparam linear
set cntrparam levels auto 5
set cntrparam points 5
set size ratio 1 1,1
set origin 0,0
set style data points
set style function lines
unset xzeroaxis
unset yzeroaxis
unset zzeroaxis
unset x2zeroaxis
unset y2zeroaxis
set xyplane relative 0.5
set tics scale 1, 0.5, 1, 1, 1
set mxtics default
set mytics default
set mztics default
set mx2tics default
set my2tics default
set mcbtics default
set mrtics default
set nomttics
set xtics border in scale 1,0.5 mirror norotate autojustify
set xtics norangelimit autofreq
set ytics border in scale 1,0.5 mirror norotate autojustify
set ytics norangelimit autofreq
set ztics border in scale 1,0.5 nomirror norotate autojustify
set ztics norangelimit autofreq
unset x2tics
unset y2tics
set cbtics border in scale 1,0.5 mirror norotate autojustify
set cbtics norangelimit autofreq
set rtics axis in scale 1,0.5 nomirror norotate autojustify
set rtics norangelimit autofreq
unset ttics
set title ""
set title font "" norotate
set timestamp bottom
set timestamp ""
set timestamp font "" norotate
set trange [ * : * ] noreverse nowriteback
set urange [ * : * ] noreverse nowriteback
set vrange [ * : * ] noreverse nowriteback
set xlabel ""
set xlabel font "" textcolor lt -1 norotate
set x2label ""
set x2label font "" textcolor lt -1 norotate
set xrange [ -2.50000 : 2.50000 ] noreverse nowriteback
set x2range [ -2.46678 : 2.46600 ] noreverse nowriteback
set ylabel ""
set ylabel font "" textcolor lt -1 rotate
set y2label ""
set y2label font "" textcolor lt -1 rotate
set yrange [ -2.50000 : 2.50000 ] noreverse nowriteback
set y2range [ -2.46600 : 2.46731 ] noreverse nowriteback
set zlabel ""
set zlabel font "" textcolor lt -1 norotate
set zrange [ * : * ] noreverse nowriteback
set cblabel ""
set cblabel font "" textcolor lt -1 rotate
set cbrange [ * : * ] noreverse nowriteback
set rlabel ""
set rlabel font "" textcolor lt -1 norotate
set rrange [ * : * ] noreverse nowriteback
unset logscale
unset jitter
set zero 1e-08
set lmargin -1
set bmargin -1
set rmargin -1
set tmargin -1
set locale "pt_BR.UTF-8"
set pm3d explicit at s
set pm3d scansautomatic
set pm3d interpolate 1,1 flush begin noftriangles noborder corners2color mean
set pm3d nolighting
set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB
set palette rgbformulae 7, 5, 15
set colorbox default
set colorbox vertical origin screen 0.9, 0.2 size screen 0.05, 0.6 front noinvert bdefault
set style boxplot candles range 1.50 outliers pt 7 separation 1 labels auto unsorted
set loadpath
set fontpath
set psdir
set fit brief errorvariables nocovariancevariables errorscaling prescale nowrap v5
GNUTERM = "wxt"
## Last datafile plotted: "internal-tpc.dat"
plot 'positions.txt' u 2:3 title 'Tracks', 'external-TPC.dat' u 1:2 w l title "TPC external radius", 'internal-tpc.dat' u 1:2 w l title "TPC internal radius"
# EOF