Commit PWG2 Makefile. To compile the code of the PWG2 module one has to do: make...
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnReader.cxx
CommitLineData
b35947a8 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 AliRsnReader
18//
19// Reader for conversion of ESD or Kinematics output into AliRsnEvent
20//
21// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
22//-------------------------------------------------------------------------
23
24#include <Riostream.h>
25
26#include <TH1.h>
27#include <TH3.h>
28#include <TFile.h>
29#include <TChain.h>
30#include <TArrayF.h>
31#include <TParticle.h>
32#include <TRandom.h>
33#include <TObjString.h>
34#include <TObjectTable.h>
35#include <TOrdCollection.h>
36
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 "AliGenEventHeader.h"
44#include "AliGenPythiaEventHeader.h"
45
46#include "AliRsnDaughter.h"
47#include "AliRsnEvent.h"
48#include "AliRsnReader.h"
49
50ClassImp(AliRsnReader)
51
52//--------------------------------------------------------------------------------------------------------
53AliRsnReader::AliRsnReader()
54//
55// Constructor.
56// Initializes all working parameters to default values:
57// - ESD particle identification
58// - rejection of non-ITS-refitted tracks
59// - maximum distance allowed from primary vertex = 3.0 cm (beam pipe)
60//
61{
62 fPIDMethod = kESDPID;
63 Int_t i;
64 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = 1.0;
65 fPtLimit4PID = 4.0;
66 fProbThreshold = 0.0;
67 fMaxRadius = 3.0; // the beam pipe
68
69 fEvents = 0;
70 for (i = AliPID::kElectron; i <= AliPID::kProton; i++) {
71 fEffPos[i] = 0;
72 fEffNeg[i] = 0;
73 }
74 fResPt = 0;
75 fResLambda = 0;
76 fResPhi = 0;
77}
78//--------------------------------------------------------------------------------------------------------
79AliRsnReader::AliRsnReader(const AliRsnReader &copy) : TObject(copy)
80//
81// Copy constructor.
82// Initializes all working parameters to same values of another simlar object.
83// TTree data member is not created.
84//
85{
86 fPIDMethod = copy.fPIDMethod;
87 Int_t i;
88 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = copy.fPrior[i];
89 fPtLimit4PID = copy.fPtLimit4PID;
90 fProbThreshold = copy.fProbThreshold;
91 fMaxRadius = copy.fMaxRadius;
92
93 fEvents = 0;
94 for (i = AliPID::kElectron; i <= AliPID::kProton; i++) {
95 fEffPos[i] = copy.fEffPos[i];
96 fEffNeg[i] = copy.fEffNeg[i];
97 }
98
99 fResPt = copy.fResPt;
100 fResLambda = copy.fResLambda;
101 fResPhi = copy.fResPhi;
102}
103//--------------------------------------------------------------------------------------------------------
104void AliRsnReader::Clear(Option_t *option)
105//
106// Clear collection of filenames.
107// If requested with the option "DELETELIST",
108// the collection object is also deleted.
109//
110{
111 TString opt(option);
112
113 if (!opt.CompareTo("TREE", TString::kIgnoreCase)) {
114 fEvents->Reset();
115 if (!opt.CompareTo("DELTREE", TString::kIgnoreCase)) {
116 delete fEvents;
117 fEvents = 0;
118 }
119 }
120
121 Int_t i;
122 for (i = AliPID::kElectron; i <= AliPID::kProton; i++) {
123 fEffPos[i] = 0;
124 fEffNeg[i] = 0;
125 }
126 fResPt = 0;
127 fResLambda = 0;
128 fResPhi = 0;
129}
130//--------------------------------------------------------------------------------------------------------
131Bool_t AliRsnReader::EffSim(Int_t pdg, Double_t pt, Double_t eta, Double_t z)
132//
133// If efficiency histogram are present, they are used to simulate efficiency.
134// Pt, Eta and Z are used to find reconstruction efficiency value to be used
135// and PDG is used to select efficiency histogram for a given particle.
136// An extraction is done, and it is supposed that particle must be accepted
137// only when this function retunrs kTRUE (= successful extraction).
138//
139{
140 // find particle sign from PDG code
141 Int_t sign;
142 if (TMath::Abs(pdg) >= 211) {
143 if (pdg > 0) sign = 1; else sign = -1;
144 }
145 else {
146 if (pdg > 0) sign = -1; else sign = 1;
147 }
148
149 // convert PDG code into one value in AliPID::kSPECIES
150 // (if returned value is not a charged particle, procedure is skipped)
151 Int_t index = FindType(pdg);
152 TH3D *eff = 0;
153 if (index >= AliPID::kElectron && index <= AliPID::kProton) {
154 if (sign > 0) eff = fEffPos[index]; else eff = fEffNeg[index];
155 }
156
157 // if no efficiency histogram is defined, method returns a 'fail' result
158 if (!eff) return kFALSE;
159
160 // otherwise, a random extraction is made
161 Int_t ibin = eff->FindBin(z, eta, pt);
162 Double_t ref = (Double_t)eff->GetBinContent(ibin);
163 Double_t ran = gRandom->Rndm();
164 return (ran <= ref);
165}
166//--------------------------------------------------------------------------------------------------------
167Double_t* AliRsnReader::GetPIDprobabilities(AliRsnDaughter track)
168//
169// Computes PID probabilites starting from priors and weights
170//
171{
172 Double_t *prob = new Double_t[AliPID::kSPECIES];
173
174 Int_t i;
175
176 // step 1 - compute the normalization factor
177 Double_t sum = 0.0;
178 for (i = 0; i < AliPID::kSPECIES; i++) {
179 prob[i] = fPrior[i] * track.GetPIDweight(i);
180 sum += prob[i];
181 }
182 if (sum <= 0.0) return 0;
183
184 // step 2 - normalize PID weights by the relative prior probability
185 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
186 prob[i] /= sum;
187 }
188
189 return prob;
190}
191//--------------------------------------------------------------------------------------------------------
192void AliRsnReader::Identify(AliRsnDaughter &track)
193//
194// Identifies a track according to method selected
195//
196{
197 Bool_t doESDPID = (fPIDMethod == kESDPID);
198 Bool_t keepRecSign = (fPIDMethod == kPerfectPIDwithRecSign);
199
200 if (doESDPID) {
201 // when doing ESD PID it is supposed that charge sign
202 // comes from ESD track and it is not modified
203 Double_t pt = track.GetPt();
204 if (pt <= fPtLimit4PID) {
205 Double_t *prob = GetPIDprobabilities(track);
206 if (!prob) track.SetPDG(0);
207 Int_t idx[AliPID::kSPECIES];
208 TMath::Sort(AliPID::kSPECIES, prob, idx);
209 Int_t imax = idx[0];
210 Double_t maxprob = prob[imax];
211 if (maxprob >= fProbThreshold) {
212 track.SetPDG((UShort_t)AliPID::ParticleCode(imax));
213 }
214 delete [] prob;
215 }
216 else {
217 track.SetPDG(0);
218 }
219 }
220 else {
221 Short_t truePDG = track.GetTruePDG();
222 track.SetPDG((UShort_t)TMath::Abs(truePDG));
223 if (!keepRecSign) {
224 if (TMath::Abs(truePDG) <= 13) {
225 if (truePDG > 0) track.SetSign(-1); else track.SetSign(1);
226 }
227 else {
228 if (truePDG > 0) track.SetSign(1); else track.SetSign(-1);
229 }
230 }
231 }
232}
233//--------------------------------------------------------------------------------------------------------
234TTree* AliRsnReader::ReadParticles(const char *path, Option_t *option)
235//
236// Opens the file "galice.root" and kinematics in the path specified,
237// loads Kinematics informations and fills and AliRsnEvent object
238// with data coming from simulation.
239// If required, an efficiency simulation can be done which cuts off some tracks
240// depending on a MonteCarlo extraction weighted on the efficiency histograms
241// passed to the class.
242// Moreover, a smearing can be done if requested, using the resolution histograms
243// passed to the class.
244// Allowed options:
245// - "E" --> do efficiency simulation
246// - "P" --> do momentum smearing
247//
248{
249 fPIDMethod = kPerfectPID;
250
251 TTree *events = new TTree("selection", "AliRsnEvents");
252 TTree::SetBranchStyle(1);
253 AliRsnEvent *event = new AliRsnEvent;
254 TBranch *branch = events->Branch("events", "AliRsnEvent", &event, 32000, 1);
255 branch->SetAutoDelete(kFALSE);
256
257 // get path
258 TString strPath(path);
259 if (strPath.Length() < 1) return 0;
260 strPath += '/';
261
262 // evaluate options
263 TString opt(option);
264 opt.ToUpper();
265 Bool_t doEffSim = (opt.Contains("E"));
266 Bool_t doSmearing = (opt.Contains("P"));
267
268 // initialize run loader
269 AliRunLoader *runLoader = OpenRunLoader(path);
270 if (!runLoader) return 0;
271 Int_t nEvents = runLoader->GetNumberOfEvents();
272 TArrayF vertex(3);
273
274 // loop on events
275 Int_t i, iev, index = 0, imum, nParticles, nStoredParticles = 0;
276 Double_t vtot, px, py, pz, arg, eta;
277 TParticle *particle = 0, *mum = 0;
278 for (iev = 0; iev < nEvents; iev++) {
279
280 // load given event
281 cout << "\rEvent " << iev << " of " << nEvents << flush;
282 runLoader->GetEvent(iev);
283
284 // primary vertex
285 runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
286 //cout << "Process type: " << ((AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader())->ProcessType() << endl;
287
288 // kinematics
289 AliStack *stack = runLoader->Stack();
290 if (!stack) continue;
291
292 // get number of particles to collect
293 nParticles = stack->GetNtrack();
294 if (!nParticles) return 0;
295
296 // create new AliRsnEvent object
297 event->Clear("DELETE");
298 event->SetESD(kFALSE);
299 event->SetPath(strPath);
300 event->Init();
301
302 // loop on tracks
303 index = 0;
304 for (i = 0; i < nParticles; i++) {
305
306 // get particle
307 particle = stack->Particle(i);
308 if (!particle) continue;
309
310 // check against maximum allowed distance from primary vertex
311 vtot = (particle->Vx() - vertex[0])*(particle->Vx() - vertex[0]);
312 vtot += (particle->Vy() - vertex[1])*(particle->Vy() - vertex[1]);
313 vtot += (particle->Vz() - vertex[2])*(particle->Vz() - vertex[2]);
314 vtot = TMath::Sqrt(vtot);
315 if (vtot > fMaxRadius) continue;
316
317 // efficiency selection
318 if (doEffSim) {
319 arg = 0.5 * TMath::ATan2(particle->Pz(), particle->Pt());
320 if (arg > 0.) eta = -TMath::Log(arg); else continue;
321 Bool_t result = EffSim(particle->GetPdgCode(), (Double_t)vertex[2], eta, (Double_t)particle->Pt());
322 if (result == kFALSE) continue;
323 }
324
325 // smear kinematic parameters (if requested)
326 px = particle->Px();
327 py = particle->Py();
328 pz = particle->Pz();
329 if (doSmearing) SmearMomentum(px, py, pz);
330
331 // add to collection of tracks
332 nStoredParticles++;
333 AliRsnDaughter *track = new AliRsnDaughter;
334 if (!track->Adopt(particle)) {
335 delete track;
336 continue;
337 }
338 track->SetIndex(index);
339 track->SetLabel(i);
340 track->SetPxPyPz(px, py, pz);
341 imum = particle->GetFirstMother();
342 if (imum >= 0) {
343 mum = stack->Particle(imum);
344 track->SetMotherPDG( (Short_t)mum->GetPdgCode() );
345 }
346 Identify(*track);
347
348 event->AddTrack(*track);
349 delete track;
350
351 } // end of loop over tracks
352
353 // compute total multiplicity
354 event->GetMultiplicity();
355
356 // link to events tree and fill
357 events->Fill();
358 }
359
360 runLoader->UnloadKinematics();
361 delete runLoader;
362
363 return events;
364}
365//--------------------------------------------------------------------------------------------------------
366TTree* AliRsnReader::ReadTracks(const char *path, Option_t *option)
367//
368// Reads AliESD in a given path, with and "experimental" method.
369// When using this method, no Kinematics information is assumed
370// and particle identification is always done with the bayesian method.
371// No Kinematics informations are stored.
372// Allowed options (actually):
373// - "R" : reject tracks which do not have the flag "kITSRefit" set to TRUE
374// - "F" : reject 'fake' tracks (negative label)
375//
376{
377 fPIDMethod = kESDPID;
378
379 TTree *events = new TTree("selection", "AliRsnEvents");
380 TTree::SetBranchStyle(1);
381 AliRsnEvent *event = new AliRsnEvent;
382 TBranch *branch = events->Branch("events", "AliRsnEvent", &event, 32000, 1);
383 branch->SetAutoDelete(kFALSE);
384
385 // get path
386 TString strPath(path);
387 if (strPath.Length() < 1) return 0;
388 strPath += '/';
389
390 // evaluate options
391 TString opt(option);
392 opt.ToUpper();
393 Bool_t checkITSrefit = (opt.Contains("R"));
394 Bool_t rejectFakes = (opt.Contains("F"));
395
396 // opens ESD file
397 TFile *fileESD = TFile::Open(Form("%s/AliESDs.root", strPath.Data()));
398 if (!fileESD) return 0;
399 if (fileESD->IsOpen() == kFALSE) return 0;
400 TTree* treeESD = (TTree*)fileESD->Get("esdTree");
401 AliESD *esd = 0;
402 treeESD->SetBranchAddress("ESD", &esd);
403 Int_t nev = (Int_t)treeESD->GetEntries();
404
405 // loop on events
406 Int_t i, nSelTracks = 0;
407 Double_t vertex[3];
408 for (i = 0; i < nev; i++) {
409
410 // message
411 cout << "\rEvent " << i << flush;
412 treeESD->GetEntry(i);
413
414 // primary vertex
415 vertex[0] = esd->GetVertex()->GetXv();
416 vertex[1] = esd->GetVertex()->GetYv();
417 vertex[2] = esd->GetVertex()->GetZv();
418
419 // create new AliRsnEvent object
420 event->Clear("DELETE");
421 event->Init();
422 event->SetPath(strPath);
423 event->SetESD();
424
425 // get number of tracks
426 Int_t ntracks = esd->GetNumberOfTracks();
427 if (!ntracks) continue;
428
429 // store tracks from ESD
430 Int_t index, label;
431 Double_t vtot, v[3];
432 for (index = 0; index < ntracks; index++) {
433
434 // get track
435 AliESDtrack *esdTrack = esd->GetTrack(index);
436
437 // check against vertex constraint
438 esdTrack->GetXYZ(v);
439 vtot = (v[0] - vertex[0])*(v[0] - vertex[0]);
440 vtot += (v[1] - vertex[1])*(v[1] - vertex[1]);
441 vtot += (v[2] - vertex[2])*(v[2] - vertex[2]);
442 vtot = TMath::Sqrt(vtot);
443 if (vtot > fMaxRadius) continue;
444
445 // check for ITS refit
446 if (checkITSrefit) {
447 if (!(esdTrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
448 }
449
450 // check for fakes
451 label = esdTrack->GetLabel();
452 if (rejectFakes && (label < 0)) continue;
453
454 // create AliRsnDaughter (and make Bayesian PID)
455 AliRsnDaughter track;
456 if (!track.Adopt(esdTrack, checkITSrefit)) continue;
457 track.SetIndex(index);
458 track.SetLabel(label);
459 Identify(track);
460
461 // store in TClonesArray
462 event->AddTrack(track);
463 nSelTracks++;
464 }
465
466 // compute total multiplicity
467 event->GetMultiplicity();
468
469 // link to events tree and fill
470 events->Fill();
471 }
472
473 fileESD->Close();
474 return events;
475}
476//--------------------------------------------------------------------------------------------------------
477TTree* AliRsnReader::ReadTracksAndParticles(const char *path, Option_t *option)
478//
479// Reads AliESD in a given path, getting also informations from Kinematics.
480// In this case, the PID method used is the one selected with apposite setter.
481// Allowed options (actually):
482// - "R" : reject tracks which do not have the flag "kITSRefit" set to TRUE
483// - "F" : reject 'fake' tracks (negative label)
484// - "M" : use 'true' momentum instead of reconstructed one
485//
486{
487 TTree *events = new TTree("selection", "AliRsnEvents");
488 TTree::SetBranchStyle(1);
489 AliRsnEvent *event = new AliRsnEvent;
490 TBranch *branch = events->Branch("events", "AliRsnEvent", &event, 32000, 1);
491 branch->SetAutoDelete(kFALSE);
492
493 // get path
494 TString strPath(path);
495 if (strPath.Length() < 1) return 0;
496
497 // evaluate options
498 TString opt(option);
499 opt.ToUpper();
500 Bool_t checkITSrefit = (opt.Contains("R"));
501 Bool_t rejectFakes = (opt.Contains("F"));
502 Bool_t copyMomentum = (opt.Contains("M"));
503
504 // opens ESD file
505 TFile *fileESD = TFile::Open(Form("%s/AliESDs.root", strPath.Data()));
506 if (!fileESD) return 0;
507 if (fileESD->IsOpen() == kFALSE) return 0;
508 TTree* treeESD = (TTree*)fileESD->Get("esdTree");
509 AliESD *esd = 0;
510 treeESD->SetBranchAddress("ESD", &esd);
511 Int_t nevRec = (Int_t)treeESD->GetEntries();
512
513 // initialize run loader
514 AliRunLoader *runLoader = OpenRunLoader(path);
515 if (!runLoader) return 0;
516 Int_t nevSim = runLoader->GetNumberOfEvents();
517
518 // check number of reconstructed and simulated events
519 if ( (nevSim != 0 && nevRec != 0) && (nevSim != nevRec) ) {
520 cerr << "Count mismatch: sim = " << nevSim << " -- rec = " << nevRec << endl;
521 return 0;
522 }
523 else if (nevSim == 0 && nevRec == 0) {
524 cerr << "Count error: sim = rec = 0" << endl;
525 return 0;
526 }
527
528 // loop on events
529 Int_t i, procType, ntracks, nSelTracks = 0;
530 Double_t vertex[3];
531 for (i = 0; i < nevRec; i++) {
532
533 // get event
534 cout << "\rEvent " << i << " " << flush;
535 treeESD->GetEntry(i);
536 runLoader->GetEvent(i);
537
538 // reject event if it is diffractive
539 procType = ((AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader())->ProcessType();
540 if (procType == 92 || procType == 93 || procType == 94) {
541 cout << "Skipping diffractive event" << endl;
542 continue;
543 }
544
545 // get particle stack
546 AliStack *stack = runLoader->Stack();
547
548 // primary vertex
549 vertex[0] = esd->GetVertex()->GetXv();
550 vertex[1] = esd->GetVertex()->GetYv();
551 vertex[2] = esd->GetVertex()->GetZv();
552
553 // multiplicity
554 ntracks = esd->GetNumberOfTracks();
555 if (!ntracks) {
556 Warning("ReadTracksAndParticles", "Event %d has no tracks!", i);
557 continue;
558 }
559
560 // create new AliRsnEvent object
561 event->Clear("DELETE");
562 event->Init();
563 event->SetPath(strPath);
564 event->SetESD();
565
566 // store tracks from ESD
567 Int_t index, label;
568 Double_t vtot, v[3];
569 for (index = 0; index < ntracks; index++) {
570
571 // get track
572 AliESDtrack *esdTrack = esd->GetTrack(index);
573
574 // check against vertex constraint
575 esdTrack->GetXYZ(v);
576 vtot = (v[0] - vertex[0])*(v[0] - vertex[0]);
577 vtot += (v[1] - vertex[1])*(v[1] - vertex[1]);
578 vtot += (v[2] - vertex[2])*(v[2] - vertex[2]);
579 vtot = TMath::Sqrt(vtot);
580 if (vtot > fMaxRadius) continue;
581
582 // check for fakes
583 label = esdTrack->GetLabel();
584 if (rejectFakes) {
585 if (label < 0) continue;
586 }
587
588 // create AliRsnDaughter (and make Bayesian PID)
589 AliRsnDaughter track;
590 if (!track.Adopt(esdTrack, checkITSrefit)) continue;
591 track.SetIndex(index);
592
593 // retrieve particle and get Kine info
594 TParticle *part = stack->Particle(TMath::Abs(label));
595 track.SetTruePDG(part->GetPdgCode());
596 Int_t mother = part->GetFirstMother();
597 track.SetMother(mother);
598 if (mother >= 0) {
599 TParticle *mum = stack->Particle(mother);
600 track.SetMotherPDG(mum->GetPdgCode());
601 }
602 if (copyMomentum) track.SetPxPyPz(part->Px(), part->Py(), part->Pz());
603
604 // identification
605 Identify(track);
606
607 // store in TClonesArray
608 event->AddTrack(track);
609 nSelTracks++;
610 }
611
612 // compute total multiplicity
613 event->GetMultiplicity();
614
615 // link to events tree and fill
616 events->Fill();
617 }
618
619 runLoader->UnloadKinematics();
620 delete runLoader;
621 fileESD->Close();
622
623 return events;
624}
625//--------------------------------------------------------------------------------------------------------
626void AliRsnReader::SetEfficiencyHistogram(AliPID::EParticleType type, TH3D *h, Bool_t pos)
627//
628// Sets efficiency histogram for a given particle species.
629// If third argument is "true", histo is assigned to positive particles,
630// otherwise it is assigned to negative particles.
631//
632{
633 if (type >= AliPID::kElectron && type < AliPID::kPhoton) {
634 if (pos) fEffPos[type] = h; else fEffNeg[type] = h;
635 }
636}
637//--------------------------------------------------------------------------------------------------------
638void AliRsnReader::SetPriorProbabilities(Double_t *prior)
639//
640// Set prior probabilities to be used in case of ESD PID.
641//
642{
643 if (!prior) return;
644
645 Int_t i = 0;
646 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = prior[i];
647}
648//--------------------------------------------------------------------------------------------------------
649void AliRsnReader::SetPriorProbability(AliPID::EParticleType type, Double_t value)
650//
651// Sets prior probability referred to a single particle species.
652//
653{
654 if (type >= AliPID::kElectron && type < AliPID::kPhoton) fPrior[type] = value;
655}
656//--------------------------------------------------------------------------------------------------------
657void AliRsnReader::SmearMomentum(Double_t &px, Double_t &py, Double_t &pz)
658//
659// Use resolution histograms to do smearing of momentum
660// (to introduce reconstruction effects)
661//
662{
663 Double_t pt = TMath::Sqrt(px*px + py*py);
664 Double_t lambda = TMath::ATan2(pz, pt);
665 Double_t phi = TMath::ATan2(py, px);
666
667 Int_t ibin;
668 Double_t ref;
669 if (fResPt) {
670 ibin = fResPt->FindBin(pt);
671 ref = (Double_t)fResPt->GetBinContent(ibin);
672 pt *= 1.0 + ref;
673 }
674 if (fResLambda) {
675 ibin = fResLambda->FindBin(pt);
676 ref = (Double_t)fResLambda->GetBinContent(ibin);
677 lambda *= 1.0 + ref;
678 }
679 if (fResPhi) {
680 ibin = fResPhi->FindBin(pt);
681 ref = (Double_t)fResPhi->GetBinContent(ibin);
682 phi *= 1.0 + ref;
683 }
684
685 px = pt * TMath::Cos(phi);
686 py = pt * TMath::Sin(phi);
687 pz = pt * TMath::Tan(lambda);
688}
689//--------------------------------------------------------------------------------------------------------
690AliPID::EParticleType AliRsnReader::FindType(Int_t pdg)
691//
692// Finds the enum value corresponding to a PDG code
693//
694{
695 pdg = TMath::Abs(pdg);
696 switch (pdg) {
697 case 11: return AliPID::kElectron; break;
698 case 13: return AliPID::kMuon; break;
699 case 211: return AliPID::kPion; break;
700 case 321: return AliPID::kKaon; break;
701 case 2212: return AliPID::kProton; break;
702 default : return AliPID::kPhoton;
703 }
704}
705//--------------------------------------------------------------------------------------------------------
706AliRunLoader* AliRsnReader::OpenRunLoader(const char *path)
707//
708// Open the Run loader with events in a given path
709//
710{
711 // clear gALICE
712 if (gAlice) {
713 delete gAlice;
714 gAlice = 0;
715 }
716
717 // initialize run loader
718 TString name(path);
719 name += "/galice.root";
720 AliRunLoader *runLoader = AliRunLoader::Open(name.Data());
721 if (runLoader) {
722 runLoader->LoadgAlice();
723 gAlice = runLoader->GetAliRun();
724 runLoader->LoadKinematics();
725 runLoader->LoadHeader();
726 }
727
728 return runLoader;
729}