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