removed eff-c++ warnings
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnSelectorRL.cxx
CommitLineData
62c607bd 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-------------------------------------------------------------------------
17// Class AliRsnSelectorRL
18// ------------------------
19// Reader for conversion of ESD output into the internal format
20// used for resonance study.
21// ---
22// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
23// ---
24// adapted for TSelector compliance
25// by : R. Vernet (email: renaud.vernet@cern.ch)
26//-------------------------------------------------------------------------
27
28#include <Riostream.h>
29
30#include <TFile.h>
31#include <TChain.h>
32#include <TParticle.h>
33#include <TRandom.h>
34#include <TObjString.h>
35#include <TObjectTable.h>
36#include <TOrdCollection.h>
37#include "AliRun.h"
38#include "AliESD.h"
39#include "AliStack.h"
40#include "AliHeader.h"
41#include "AliESDtrack.h"
42#include "AliRunLoader.h"
43#include "AliRsnDaughter.h"
44#include "AliRsnEvent.h"
45#include "AliRsnSelectorRL.h"
46#include "TAlienFile.h"
47
48ClassImp(AliRsnSelectorRL)
49
50//--------------------------------------------------------------------------------------------------------
51AliRsnSelectorRL::AliRsnSelectorRL(TTree*) :
52 AliSelectorRL(),
53 fOutputPath(0),
54 fDebugFlag(0),
55 fStoreKineInfo(0),
56 fCheckITSRefit(0),
57 fRejectFakes(0),
58 fCopyMomentum(0),
59 fIsRunLoaderOpen(0),
60 fRsnEventTree(0),
61 fRsnEvent(0),
c37c6481 62 fRsnEventBranch(0),
63 fPIDMethod(kESDPID),
64 fPtLimit4PID(4.0),
65 fProbThreshold(0.0),
66 fMaxRadius(3.0)
62c607bd 67{
68//
69// Constructor.
70// Initializes all working parameters to default values:
71// - ESD particle identification
72// - rejection of non-ITS-refitted tracks
73// - maximum distance allowed from primary vertex = 3.0 cm (beam pipe)
74//
c37c6481 75 Int_t i;
76 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = 1.0;
62c607bd 77}
78
79AliRsnSelectorRL::~AliRsnSelectorRL() {
80 Clear();
81}
82//--------------------------------------------------------------------------------------------------------
83AliRsnSelectorRL::AliRsnSelectorRL(const AliRsnSelectorRL& obj) :
84 AliSelectorRL(), // not implemented a copy constructor for AliRsnSelectorRL
85 fOutputPath(obj.fOutputPath),
86 fDebugFlag(obj.fDebugFlag),
87 fStoreKineInfo(obj.fStoreKineInfo),
88 fCheckITSRefit(obj.fCheckITSRefit),
89 fRejectFakes(obj.fRejectFakes),
90 fCopyMomentum(obj.fCopyMomentum),
91 fIsRunLoaderOpen(obj.fIsRunLoaderOpen),
92 fRsnEventTree(obj.fRsnEventTree),
93 fRsnEvent(obj.fRsnEvent),
c37c6481 94 fRsnEventBranch(obj.fRsnEventBranch),
95 fPIDMethod(obj.fPIDMethod),
96 fPtLimit4PID(obj.fPtLimit4PID),
97 fProbThreshold(obj.fProbThreshold),
98 fMaxRadius(obj.fMaxRadius)
62c607bd 99{
100//
101// Copy constructor
102//
c37c6481 103 Int_t i;
104 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = obj.fPrior[i];
62c607bd 105}
106//--------------------------------------------------------------------------------------------------------
107AliRsnSelectorRL& AliRsnSelectorRL::operator=(const AliRsnSelectorRL& obj)
108{
109//
110// Assignment operator
111// works in the same way as the copy constructor
112//
113 if (this!=&obj) {
114 fDebugFlag = obj.fDebugFlag;
115 fOutputPath = obj.fOutputPath;
116 fStoreKineInfo = obj.fStoreKineInfo;
117 fCheckITSRefit = obj.fCheckITSRefit;
118 fRejectFakes = obj.fRejectFakes;
119 fCopyMomentum = obj.fCopyMomentum;
120 fIsRunLoaderOpen = obj.fIsRunLoaderOpen;
121 fRsnEventTree = obj.fRsnEventTree;
122 fRsnEvent = obj.fRsnEvent;
123 fRsnEventBranch = obj.fRsnEventBranch;
124 fPIDMethod=obj.fPIDMethod;
125 for (Int_t i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = obj.fPrior[i];
126 fPtLimit4PID = obj.fPtLimit4PID;
127 fProbThreshold = obj.fProbThreshold;
128 fMaxRadius = obj.fMaxRadius;
129 }
130
131 return *this;
132}
133//--------------------------------------------------------------------------------------------------------
134void AliRsnSelectorRL::Clear(Option_t * /*option*/)
135{
136 //
137 // Does nothing.
138 //
139}
140//--------------------------------------------------------------------------------------------------------
141Double_t* AliRsnSelectorRL::GetPIDprobabilities(AliRsnDaughter track)
142{
143//
144// Computes PID probabilites starting from priors and weights
145//
146 Int_t i;
147 Double_t *prob = new Double_t[AliPID::kSPECIES];
148
149 // step 1 - compute the normalization factor
150 Double_t sum = 0.0;
151 for (i = 0; i < AliPID::kSPECIES; i++) {
152 prob[i] = fPrior[i] * track.GetPIDweight(i);
153 sum += prob[i];
154 }
155 if (sum <= 0.0) return 0;
156
157 // step 2 - normalize PID weights by the relative prior probability
158 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
159 prob[i] /= sum;
160 }
161
162 return prob;
163}
164//--------------------------------------------------------------------------------------------------------
165void AliRsnSelectorRL::Identify(AliRsnDaughter &track)
166{
167//
168// Identifies a track according to method selected
169//
170 Bool_t doESDPID = (fPIDMethod == kESDPID);
171 Bool_t keepRecSign = (fPIDMethod == kPerfectPIDwithRecSign);
172
173 if (doESDPID) {
174 // when doing ESD PID it is supposed that charge sign
175 // comes from ESD track and it is not modified
176 Double_t pt = track.GetPt();
177 if (pt <= fPtLimit4PID) {
178 Double_t *prob = GetPIDprobabilities(track);
179 if (!prob) track.SetPDG(0);
180 Int_t idx[AliPID::kSPECIES];
181 TMath::Sort(AliPID::kSPECIES, prob, idx);
182 Int_t imax = idx[0];
183 Double_t maxprob = prob[imax];
184 if (maxprob >= fProbThreshold) {
185 track.SetPDG((UShort_t)AliPID::ParticleCode(imax));
186 }
187 delete [] prob;
188 }
189 else {
190 track.SetPDG(0);
191 }
192 }
193 else {
194 Short_t truePDG = track.GetTruePDG();
195 track.SetPDG((UShort_t)TMath::Abs(truePDG));
196 if (!keepRecSign) {
197 if (TMath::Abs(truePDG) <= 13) {
198 if (truePDG > 0) track.SetSign(-1); else track.SetSign(1);
199 }
200 else {
201 if (truePDG > 0) track.SetSign(1); else track.SetSign(-1);
202 }
203 }
204 }
205}
206//--------------------------------------------------------------------------------------------------------
207void AliRsnSelectorRL::SetPriorProbabilities(Double_t *prior)
208{
209//
210// Set prior probabilities to be used in case of ESD PID.
211//
212 if (!prior) return;
213 Int_t i = 0;
214 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = prior[i];
215}
216//--------------------------------------------------------------------------------------------------------
217void AliRsnSelectorRL::SetPriorProbability(AliPID::EParticleType type, Double_t value)
218{
219//
220// Sets prior probability referred to a single particle species.
221//
222 if (type >= AliPID::kElectron && type < AliPID::kPhoton) fPrior[type] = value;
223}
224//--------------------------------------------------------------------------------------------------------
225AliPID::EParticleType AliRsnSelectorRL::FindType(Int_t pdg)
226{
227//
228// Finds the enum value corresponding to a PDG code
229//
230 pdg = TMath::Abs(pdg);
231 switch (pdg) {
232 case 11: return AliPID::kElectron; break;
233 case 13: return AliPID::kMuon; break;
234 case 211: return AliPID::kPion; break;
235 case 321: return AliPID::kKaon; break;
236 case 2212: return AliPID::kProton; break;
237 default : return AliPID::kPhoton;
238 }
239}
240//--------------------------------------------------------------------------------------------------------
241
242//--------------------------------------------------------
243// The following is Selector stuff
244//--------------------------------------------------------
245
246//--------------------------------------------------------------------------------------------------------
247void AliRsnSelectorRL::Begin(TTree *)
248{
249//
250// Implementation of BEGIN method
251//
252 Info("Begin", "");
253 TString option = GetOption();
254}
255//--------------------------------------------------------------------------------------------------------
256void AliRsnSelectorRL::SlaveBegin(TTree * tree)
257{
258//
259// Implementation of secondary BEGIN which
260// is called by separate process managers
261//
262 Info("SlaveBegin", "");
263 Init(tree);
264 TString option = GetOption();
265}
266//--------------------------------------------------------------------------------------------------------
267void AliRsnSelectorRL::Init(TTree *tree)
268{
269//
270// Initializer
271// Connects the selector to a TTree and links the branch
272// which is used to make translation ESD --> Rsn
273//
274 Info("Init","");
275
276 if (!tree) return;
277 fTree = tree;
278 if (fDebugFlag) {
279 Info("Init", "fTree=%p fTree->GetCurrentFile()=%p", fTree, fTree->GetCurrentFile());
280 }
281 fTree->SetBranchAddress("ESD", &fESD);
282 fRsnEventTree = new TTree("selection", "AliRsnEvents");
283 TTree::SetBranchStyle(1);
284 fRsnEvent = new AliRsnEvent;
285 fRsnEventBranch = fRsnEventTree->Branch("events", "AliRsnEvent", &fRsnEvent, 32000, 1);
286 fRsnEventBranch->SetAutoDelete(kFALSE);
287}
288//--------------------------------------------------------------------------------------------------------
289Bool_t AliRsnSelectorRL::Process(Long64_t entry)
290{
291//
292// Main core of the Selector processing
293// Reads the ESD input and creates a TTree or AliRsnEvent objects
294//
295 if (fDebugFlag) Info("Process", "Processing event %d", entry);
296 if (!AliSelectorRL::Process(entry)) return kFALSE;
297
298 AliStack* stack = GetStack();
299 if (!stack) {
300 Warning("Process", "NULL stack: cannot get kinematics info");
301 fStoreKineInfo = kFALSE;
302 }
303
304 Int_t ntracks, nSelTracks = 0;
305
306 // primary vertex
307 Double_t vertex[3];
308 vertex[0] = fESD->GetVertex()->GetXv();
309 vertex[1] = fESD->GetVertex()->GetYv();
310 vertex[2] = fESD->GetVertex()->GetZv();
311
312 // multiplicity
313 ntracks = fESD->GetNumberOfTracks();
314 if (!ntracks) {
315 Warning("Process", "Event %d has no tracks!", entry);
316 return kFALSE;
317 }
318 if (fDebugFlag) Info("Process", "Number of ESD tracks : %d", ntracks);
319
320 // create new AliRsnEvent object
321 fRsnEvent->Clear("DELETE");
322 fRsnEvent->Init();
323
324 // store tracks from ESD
325 Int_t index, label;
326 Double_t vtot, v[3];
327 AliRsnDaughter track;
328 for (index = 0; index < ntracks; index++) {
329 if (fDebugFlag) Info("Process","Track # %d", index);
330 AliESDtrack *esdTrack = fESD->GetTrack(index);
331
332 // check against vertex constraint
333 esdTrack->GetXYZ(v);
334 vtot = (v[0] - vertex[0])*(v[0] - vertex[0]);
335 vtot += (v[1] - vertex[1])*(v[1] - vertex[1]);
336 vtot += (v[2] - vertex[2])*(v[2] - vertex[2]);
337 vtot = TMath::Sqrt(vtot);
338 if (vtot > fMaxRadius) continue;
339
340 // check for fakes
341 label = esdTrack->GetLabel();
342 if (fRejectFakes && (label < 0)) continue;
343
344 // copy ESDtrack data into RsnDaughter (and make Bayesian PID)
345 if (!track.Adopt(esdTrack, fCheckITSRefit)) continue;
346 track.SetIndex(index);
347
348 // if requested, get Kine info
349 if (fStoreKineInfo) {
350 if (fDebugFlag) Info("Process", "Getting part label=%d stack=%p", label, stack);
351 TParticle *part = stack->Particle(TMath::Abs(label));
352 if (fDebugFlag) Info("Process", "part=%p", part);
353 track.SetTruePDG(part->GetPdgCode());
354 Int_t mother = part->GetFirstMother();
355 track.SetMother(mother);
356 if (mother >= 0) {
357 TParticle *mum = stack->Particle(mother);
358 track.SetMotherPDG(mum->GetPdgCode());
359 }
360 // if requested, reconstructed momentum is replaced with true momentum
361 if (fCopyMomentum) track.SetPxPyPz(part->Px(), part->Py(), part->Pz());
362 }
363
364 // identification
365 Identify(track);
366
367 // store in TClonesArray
368 if (fDebugFlag) track.Print("SLITVPMNW");
369 fRsnEvent->AddTrack(track);
370 nSelTracks++;
371 }
372
373 // compute total multiplicity
374 fRsnEvent->ComputeMultiplicity();
375
376 // link to events tree and fill
377 fRsnEventTree->Fill();
378
379 return kTRUE;
380}
381//--------------------------------------------------------------------------------------------------------
382void AliRsnSelectorRL::SlaveTerminate()
383{
384//
385// SlaveTerminate
386// Partial termination method
387//
388 Info("SlaveTerminate", "");
389
390 // Add the histograms to the output on each slave server
391 // fOutput->Add(fRsnEventTree);
392}
393//--------------------------------------------------------------------------------------------------------
394void AliRsnSelectorRL::Terminate()
395{
396//
397// Global termination method
398//
399 Info("Terminate","");
400 // fRsnEventTree = dynamic_cast<TTree*>(fOutput->FindObject("fRsnEventTree"));
401
402 //AliSelector::Terminate();
403 cout << fOutputPath << endl;
404 Info("Terminate", Form("Saving in: %s", fOutputPath->Data()));
405 TFile* file = TFile::Open(fOutputPath->Data(), "RECREATE");
406 fRsnEventTree->Write();
407 file->Close();
408
409 delete fRsnEventTree;
410 delete file;
411 delete fOutputPath;
412}