Adding AliRsnReader in the RESONANCES directory
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnReader.cxx
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
50 ClassImp(AliRsnReader)
51
52 //--------------------------------------------------------------------------------------------------------
53 AliRsnReader::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 //--------------------------------------------------------------------------------------------------------
79 AliRsnReader::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 //--------------------------------------------------------------------------------------------------------
104 void 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 //--------------------------------------------------------------------------------------------------------
131 Bool_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 //--------------------------------------------------------------------------------------------------------
167 Double_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 //--------------------------------------------------------------------------------------------------------
192 void 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 //--------------------------------------------------------------------------------------------------------
234 TTree* 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 //--------------------------------------------------------------------------------------------------------
366 TTree* 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 //--------------------------------------------------------------------------------------------------------
477 TTree* 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 //--------------------------------------------------------------------------------------------------------
626 void 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 //--------------------------------------------------------------------------------------------------------
638 void 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 //--------------------------------------------------------------------------------------------------------
649 void 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 //--------------------------------------------------------------------------------------------------------
657 void 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 //--------------------------------------------------------------------------------------------------------
690 AliPID::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 //--------------------------------------------------------------------------------------------------------
706 AliRunLoader* 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 }