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