1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-------------------------------------------------------------------------
17 // Class AliRsnSelectorRL
18 // ------------------------
19 // Reader for conversion of ESD output into the internal format
20 // used for resonance study.
22 // author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
24 // adapted for TSelector compliance
25 // by : R. Vernet (email: renaud.vernet@cern.ch)
26 //-------------------------------------------------------------------------
28 #include <Riostream.h>
32 #include <TParticle.h>
34 #include <TObjString.h>
35 #include <TObjectTable.h>
36 #include <TOrdCollection.h>
40 #include "AliHeader.h"
41 #include "AliESDtrack.h"
42 #include "AliRunLoader.h"
43 #include "AliRsnDaughter.h"
44 #include "AliRsnEvent.h"
45 #include "AliRsnSelectorRL.h"
46 #include "TAlienFile.h"
48 ClassImp(AliRsnSelectorRL)
50 //--------------------------------------------------------------------------------------------------------
51 AliRsnSelectorRL::AliRsnSelectorRL(TTree*) :
66 // Initializes all working parameters to default values:
67 // - ESD particle identification
68 // - rejection of non-ITS-refitted tracks
69 // - maximum distance allowed from primary vertex = 3.0 cm (beam pipe)
72 for (Int_t i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = 1.0;
75 fMaxRadius = 3.0; // the beam pipe
78 AliRsnSelectorRL::~AliRsnSelectorRL() {
81 //--------------------------------------------------------------------------------------------------------
82 AliRsnSelectorRL::AliRsnSelectorRL(const AliRsnSelectorRL& obj) :
83 AliSelectorRL(), // not implemented a copy constructor for AliRsnSelectorRL
84 fOutputPath(obj.fOutputPath),
85 fDebugFlag(obj.fDebugFlag),
86 fStoreKineInfo(obj.fStoreKineInfo),
87 fCheckITSRefit(obj.fCheckITSRefit),
88 fRejectFakes(obj.fRejectFakes),
89 fCopyMomentum(obj.fCopyMomentum),
90 fIsRunLoaderOpen(obj.fIsRunLoaderOpen),
91 fRsnEventTree(obj.fRsnEventTree),
92 fRsnEvent(obj.fRsnEvent),
93 fRsnEventBranch(obj.fRsnEventBranch)
98 fPIDMethod = obj.fPIDMethod;
99 for (Int_t i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = obj.fPrior[i];
100 fPtLimit4PID = obj.fPtLimit4PID;
101 fProbThreshold = obj.fProbThreshold;
102 fMaxRadius = obj.fMaxRadius;
104 //--------------------------------------------------------------------------------------------------------
105 AliRsnSelectorRL& AliRsnSelectorRL::operator=(const AliRsnSelectorRL& obj)
108 // Assignment operator
109 // works in the same way as the copy constructor
112 fDebugFlag = obj.fDebugFlag;
113 fOutputPath = obj.fOutputPath;
114 fStoreKineInfo = obj.fStoreKineInfo;
115 fCheckITSRefit = obj.fCheckITSRefit;
116 fRejectFakes = obj.fRejectFakes;
117 fCopyMomentum = obj.fCopyMomentum;
118 fIsRunLoaderOpen = obj.fIsRunLoaderOpen;
119 fRsnEventTree = obj.fRsnEventTree;
120 fRsnEvent = obj.fRsnEvent;
121 fRsnEventBranch = obj.fRsnEventBranch;
122 fPIDMethod=obj.fPIDMethod;
123 for (Int_t i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = obj.fPrior[i];
124 fPtLimit4PID = obj.fPtLimit4PID;
125 fProbThreshold = obj.fProbThreshold;
126 fMaxRadius = obj.fMaxRadius;
131 //--------------------------------------------------------------------------------------------------------
132 void AliRsnSelectorRL::Clear(Option_t * /*option*/)
138 //--------------------------------------------------------------------------------------------------------
139 Double_t* AliRsnSelectorRL::GetPIDprobabilities(AliRsnDaughter track)
142 // Computes PID probabilites starting from priors and weights
145 Double_t *prob = new Double_t[AliPID::kSPECIES];
147 // step 1 - compute the normalization factor
149 for (i = 0; i < AliPID::kSPECIES; i++) {
150 prob[i] = fPrior[i] * track.GetPIDweight(i);
153 if (sum <= 0.0) return 0;
155 // step 2 - normalize PID weights by the relative prior probability
156 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
162 //--------------------------------------------------------------------------------------------------------
163 void AliRsnSelectorRL::Identify(AliRsnDaughter &track)
166 // Identifies a track according to method selected
168 Bool_t doESDPID = (fPIDMethod == kESDPID);
169 Bool_t keepRecSign = (fPIDMethod == kPerfectPIDwithRecSign);
172 // when doing ESD PID it is supposed that charge sign
173 // comes from ESD track and it is not modified
174 Double_t pt = track.GetPt();
175 if (pt <= fPtLimit4PID) {
176 Double_t *prob = GetPIDprobabilities(track);
177 if (!prob) track.SetPDG(0);
178 Int_t idx[AliPID::kSPECIES];
179 TMath::Sort(AliPID::kSPECIES, prob, idx);
181 Double_t maxprob = prob[imax];
182 if (maxprob >= fProbThreshold) {
183 track.SetPDG((UShort_t)AliPID::ParticleCode(imax));
192 Short_t truePDG = track.GetTruePDG();
193 track.SetPDG((UShort_t)TMath::Abs(truePDG));
195 if (TMath::Abs(truePDG) <= 13) {
196 if (truePDG > 0) track.SetSign(-1); else track.SetSign(1);
199 if (truePDG > 0) track.SetSign(1); else track.SetSign(-1);
204 //--------------------------------------------------------------------------------------------------------
205 void AliRsnSelectorRL::SetPriorProbabilities(Double_t *prior)
208 // Set prior probabilities to be used in case of ESD PID.
212 for (i = 0; i < AliPID::kSPECIES; i++) fPrior[i] = prior[i];
214 //--------------------------------------------------------------------------------------------------------
215 void AliRsnSelectorRL::SetPriorProbability(AliPID::EParticleType type, Double_t value)
218 // Sets prior probability referred to a single particle species.
220 if (type >= AliPID::kElectron && type < AliPID::kPhoton) fPrior[type] = value;
222 //--------------------------------------------------------------------------------------------------------
223 AliPID::EParticleType AliRsnSelectorRL::FindType(Int_t pdg)
226 // Finds the enum value corresponding to a PDG code
228 pdg = TMath::Abs(pdg);
230 case 11: return AliPID::kElectron; break;
231 case 13: return AliPID::kMuon; break;
232 case 211: return AliPID::kPion; break;
233 case 321: return AliPID::kKaon; break;
234 case 2212: return AliPID::kProton; break;
235 default : return AliPID::kPhoton;
238 //--------------------------------------------------------------------------------------------------------
240 //--------------------------------------------------------
241 // The following is Selector stuff
242 //--------------------------------------------------------
244 //--------------------------------------------------------------------------------------------------------
245 void AliRsnSelectorRL::Begin(TTree *)
248 // Implementation of BEGIN method
251 TString option = GetOption();
253 //--------------------------------------------------------------------------------------------------------
254 void AliRsnSelectorRL::SlaveBegin(TTree * tree)
257 // Implementation of secondary BEGIN which
258 // is called by separate process managers
260 Info("SlaveBegin", "");
262 TString option = GetOption();
264 //--------------------------------------------------------------------------------------------------------
265 void AliRsnSelectorRL::Init(TTree *tree)
269 // Connects the selector to a TTree and links the branch
270 // which is used to make translation ESD --> Rsn
277 Info("Init", "fTree=%p fTree->GetCurrentFile()=%p", fTree, fTree->GetCurrentFile());
279 fTree->SetBranchAddress("ESD", &fESD);
280 fRsnEventTree = new TTree("selection", "AliRsnEvents");
281 TTree::SetBranchStyle(1);
282 fRsnEvent = new AliRsnEvent;
283 fRsnEventBranch = fRsnEventTree->Branch("events", "AliRsnEvent", &fRsnEvent, 32000, 1);
284 fRsnEventBranch->SetAutoDelete(kFALSE);
286 //--------------------------------------------------------------------------------------------------------
287 Bool_t AliRsnSelectorRL::Process(Long64_t entry)
290 // Main core of the Selector processing
291 // Reads the ESD input and creates a TTree or AliRsnEvent objects
293 if (fDebugFlag) Info("Process", "Processing event %d", entry);
294 if (!AliSelectorRL::Process(entry)) return kFALSE;
296 AliStack* stack = GetStack();
298 Warning("Process", "NULL stack: cannot get kinematics info");
299 fStoreKineInfo = kFALSE;
302 Int_t ntracks, nSelTracks = 0;
306 vertex[0] = fESD->GetVertex()->GetXv();
307 vertex[1] = fESD->GetVertex()->GetYv();
308 vertex[2] = fESD->GetVertex()->GetZv();
311 ntracks = fESD->GetNumberOfTracks();
313 Warning("Process", "Event %d has no tracks!", entry);
316 if (fDebugFlag) Info("Process", "Number of ESD tracks : %d", ntracks);
318 // create new AliRsnEvent object
319 fRsnEvent->Clear("DELETE");
322 // store tracks from ESD
325 AliRsnDaughter track;
326 for (index = 0; index < ntracks; index++) {
327 if (fDebugFlag) Info("Process","Track # %d", index);
328 AliESDtrack *esdTrack = fESD->GetTrack(index);
330 // check against vertex constraint
332 vtot = (v[0] - vertex[0])*(v[0] - vertex[0]);
333 vtot += (v[1] - vertex[1])*(v[1] - vertex[1]);
334 vtot += (v[2] - vertex[2])*(v[2] - vertex[2]);
335 vtot = TMath::Sqrt(vtot);
336 if (vtot > fMaxRadius) continue;
339 label = esdTrack->GetLabel();
340 if (fRejectFakes && (label < 0)) continue;
342 // copy ESDtrack data into RsnDaughter (and make Bayesian PID)
343 if (!track.Adopt(esdTrack, fCheckITSRefit)) continue;
344 track.SetIndex(index);
346 // if requested, get Kine info
347 if (fStoreKineInfo) {
348 if (fDebugFlag) Info("Process", "Getting part label=%d stack=%p", label, stack);
349 TParticle *part = stack->Particle(TMath::Abs(label));
350 if (fDebugFlag) Info("Process", "part=%p", part);
351 track.SetTruePDG(part->GetPdgCode());
352 Int_t mother = part->GetFirstMother();
353 track.SetMother(mother);
355 TParticle *mum = stack->Particle(mother);
356 track.SetMotherPDG(mum->GetPdgCode());
358 // if requested, reconstructed momentum is replaced with true momentum
359 if (fCopyMomentum) track.SetPxPyPz(part->Px(), part->Py(), part->Pz());
365 // store in TClonesArray
366 if (fDebugFlag) track.Print("SLITVPMNW");
367 fRsnEvent->AddTrack(track);
371 // compute total multiplicity
372 fRsnEvent->ComputeMultiplicity();
374 // link to events tree and fill
375 fRsnEventTree->Fill();
379 //--------------------------------------------------------------------------------------------------------
380 void AliRsnSelectorRL::SlaveTerminate()
384 // Partial termination method
386 Info("SlaveTerminate", "");
388 // Add the histograms to the output on each slave server
389 // fOutput->Add(fRsnEventTree);
391 //--------------------------------------------------------------------------------------------------------
392 void AliRsnSelectorRL::Terminate()
395 // Global termination method
397 Info("Terminate","");
398 // fRsnEventTree = dynamic_cast<TTree*>(fOutput->FindObject("fRsnEventTree"));
400 //AliSelector::Terminate();
401 cout << fOutputPath << endl;
402 Info("Terminate", Form("Saving in: %s", fOutputPath->Data()));
403 TFile* file = TFile::Open(fOutputPath->Data(), "RECREATE");
404 fRsnEventTree->Write();
407 delete fRsnEventTree;