AliAnalysisTaskMyTask.cxx 9.62 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
/**************************************************************************
 * 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
37
Int_t EventID = 0; // is equal to esd-event-id until the selected event (one with an appropriate number of tracks) is reached
38
39
40
41
42
43
44
45

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(),
46
    fESD(0), fOutputList(0), fHistPt(0),    fHistEvents(0), fHistMass(0)
47
48
49
50
51
52
{
    // 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),
53
    fESD(0), fOutputList(0), fHistPt(0),  fHistEvents(0), fHistMass(0)
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
{
    // 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!

94
95
96
97
98
99
    // 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);

100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

    // 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();

142
    // Add VERTEZES (x, y, z) and magnetic field to esd-summary.dat file
143
144
145
146
147
148
149
150
    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)

*/

151
152
153
154
155
    if(EventID == esd_event_id) {
	if(iTracks >= 15 && iTracks <= 30) {EventID = -1;}
	else {EventID++;}
    } else {EventID = -2;}

156
157
158
159
160
161
162
163
    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.
164

165
166
	Double_t Px = track->Px();
	Double_t Py = track->Py();
167
	Double_t Pt = track->Pt(); // transversal momentum, in case we need it
168
	Double_t Pz = track->Pz();
169
170
171
172
173
174

	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;

175
176
	Double_t Charge = track->Charge();

177
        // Add MASS, CHARGE and VELOCITIES (x, y, z) to esd-detail.dat file
178
	detail << fixed << Mass << " " << Charge << " ";
179
	detail << fixed << Px << " " << Py << " " << Pz << " " << Vx << " " << Vy << " " << Vz << endl;
180
181
182

	//'fixed' fixes the number of decimal places so numbers are vertically aligned

183
184
185
186
187
188
189
190
        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);
	
	}

191
192
    }

193

194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
    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)
}
//_____________________________________________________________________________