PWG2/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 }