AliAnalysisTaskMyTask.cxx 11.9 KB
Newer Older
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
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
37
/**************************************************************************
 * 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

38
Double_t highAvgPz, lowAvgPz, highAvgPzEvent, lowAvgPzEvent; //variables to store info about highest and lowest <Pz> values
39
40
Double_t highAvgPt, highAvgPtEvent; //variables to store info about highest and lowest <Pt> values

41

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
42
43
44
45
46
47
48
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(),
49
    fESD(0), fOutputList(0), fHistPt(0), fHistAvgPz(0), fHistAvgPt(0), fHistMass(0)
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
50
51
52
53
54
55
{
    // 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),
56
    fESD(0), fOutputList(0), fHistPt(0), fHistAvgPz(0), fHistAvgPt(0), fHistMass(0)
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
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
{
    // 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
93
    fHistPt = new TH1F("fHistPt", "fHistPt", 100, 0, 10);       // create histogram
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
94
95
96
    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!

97
    // |<Pz>| histogram: absolute value of average Pz (or Pz per track) for each event
98
    fHistAvgPz = new TH1F("fHistAvgPz", "fHistAvgPz", 100, 0, 10);       // create histogram
99
100
    fOutputList->Add(fHistAvgPz);

101
102
103
104
    // <Pt> histogram: value of average Pt (or Pt per track) for each event
    fHistAvgPt = new TH1F("fHistAvgPt", "fHistAvgPt", 100, 0, 10);       // create histogram
    fOutputList->Add(fHistAvgPt);

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
105
106
107
108
109
110
111
112
113
114
115
    // 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);

    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
}
//_____________________________________________________________________________
116
void AliAnalysisTaskMyTask::export_to_our_ESD_textual_format (Int_t selectedEventID)
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
117
{
118
      ofstream esd_detail, uniqueid;
119
120
121
      std::stringstream esd_filename;
      esd_filename << "esd_detail-event_" << selectedEventID << ".dat";
      esd_detail.open (esd_filename.str(),std::ofstream::app);
122
      uniqueid.open ("uniqueid.txt");
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
123
124


125
      fESD = dynamic_cast<AliESDEvent*>(InputEvent());    // get an event (called fESD) from the input file
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
126
127


128
129
130
      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
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
131

132
133
134
135
136
137
138
139
140
141
142
143
144
145

      // ___________________________________________
      // Write ESD unique_id information on file

      Int_t run_number = fESD->GetRunNumber();
      Int_t orbit_number = fESD->GetOrbitNumber();
      Int_t bunch_cross_number = fESD->GetBunchCrossNumber();

      if(run_number != 0 && orbit_number != 0 && bunch_cross_number != 0) { // If all numbers are different than zero, write'em to file
        uniqueid << "Run" << run_number << "_Orbit" << orbit_number << "_BunchCross" << bunch_cross_number << endl;
      }

      //___________________________________________

146
      Int_t iTracks(fESD->GetNumberOfTracks());           // see how many tracks there are in the event
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
147

148
149
150
151
      Double_t Vx = 0.01 * fESD->GetPrimaryVertex()->GetX();	// gets vertexes from individual events, in METERS
      Double_t Vy = 0.01 * fESD->GetPrimaryVertex()->GetY();
      Double_t Vz = 0.01 * fESD->GetPrimaryVertex()->GetZ();
      Double_t MagneticField = 0.1 * fESD->GetMagneticField();	// gets magnetic field, in TESLA
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
152

153
  /*
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
154

155
  Assumed Units: Mass (GeV/c^2) || Energy (GeV) || Momentum (GeV/c) || Charge (* 1.6*10^-19 C)
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
156

157
  */
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
158

159
      if(selectedEventID == esd_event_id) { // when we get to the selected event, fill histograms and write data
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
160

161
        Double_t PzSum = 0;
162
        Double_t absPzSum = 0;
163

164
165
        Double_t PtSum = 0;

166
        for(Int_t i(0); i < iTracks; i++) {                 // loop over all these tracks
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
167

168
          		AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(i));         // get a track (type AliESDtrack) from the event
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
169

170
          		if(!track) continue;                            // if we failed, skip this track
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
171

172
173
          		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.
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
174

175
176
          		Double_t Px = track->Px();
          		Double_t Py = track->Py();
177
          		Double_t Pt = track->Pt(); // transversal momentum
178
          		Double_t Pz = track->Pz();
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
179

180
              PzSum += Pz/iTracks; // Pz sum for |<Pz>| histogram
181
              absPzSum += abs(Pz)/iTracks; //Remember: in C++, abs function overloads
182

183
184
              PtSum += Pt/iTracks; // Pt sum for <Pt> histogram

185
186
187
188
189
          		Double_t Charge = track->Charge();

            	// Add VERTEX (x, y, z), MASS, CHARGE and MOMENTUM (x, y, z) to esd-detail.dat file
            	esd_detail << Vx << " " << Vy << " " << Vz << " ";
            	esd_detail << Mass << " " << Charge << " ";
190
            	esd_detail << Px << " " << Py << " " << Pz << " " << Pt << endl;
191
192
193
194

            	fHistPt->Fill(Pt);                     // plot the pt value of the track in a histogram

            	fHistMass->Fill(Mass);
Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
195

196
          }
197
198
199
200
201

          if(absPzSum != 0) { // This will only fill |<Pz>| histogram for events with tracks (non-empty)
            fHistAvgPz->Fill(abs(PzSum));
          }

202
203
204
205
          if(PtSum != 0) { // This will only fill <Pt> histogram for events with tracks (non-empty)
            fHistAvgPt->Fill(abs(PtSum));
          }

206
          if(selectedEventID == 0) {
207
208

            // Pz:
209
210
211
212
            highAvgPzEvent = 0;
            lowAvgPzEvent = 0;
            highAvgPz = PzSum;
            lowAvgPz = PzSum;
213
214
215
216
217

            // Pt:
            highAvgPtEvent = 0;
            highAvgPt = PtSum;

218
219
220
221
222
223
224
225
226
227
228
          } else {

            if(PzSum>highAvgPz) {
              highAvgPz = PzSum;
              highAvgPzEvent = selectedEventID;
            }
            if(PzSum<lowAvgPz) {
              lowAvgPz = PzSum;
              lowAvgPzEvent = selectedEventID;
            }

229
230
231
232
233
            if(PtSum>highAvgPt) {
              highAvgPt = PtSum;
              highAvgPtEvent = selectedEventID;
            }

234
235
          }

236
237
238
      }

      esd_detail.close();
239
      uniqueid.close();
240
}
241

242
243
244
245
246
247
248
249
250
//_____________________________________________________________________________
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

251
    export_to_our_ESD_textual_format(esd_event_id);
252

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
253
    esd_event_id++; // Increment global esd_event_id
254

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
255
256
257
258
259
260
261
262
263
264
                                                       // 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)
265
    cout << endl << endl << "Lowest Pz Mean (<Pz>) = " << lowAvgPz << "    at Event " << lowAvgPzEvent;
266
267
268
269
    cout << endl << "Highest Pz Mean (<Pz>) = " << highAvgPz << "    at Event " << highAvgPzEvent << endl;

    cout << endl << "Highest Pt Mean (<Pt>) = " << highAvgPt << "    at Event " << highAvgPtEvent << endl << endl;

Breno Rilho Lemos's avatar
Breno Rilho Lemos committed
270
271
}
//_____________________________________________________________________________