]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/Cascades/AliAnalysisTaskExtractPerformanceCascade.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
-> Changes to triggering scheme for p-Pb
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Cascades / AliAnalysisTaskExtractPerformanceCascade.cxx
CommitLineData
76029adc 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//
18// Modified version of AliAnalysisTaskCheckCascade.cxx.
19// This is a 'hybrid' output version, in that it uses a classic TTree
20// ROOT object to store the candidates, plus a couple of histograms filled on
21// a per-event basis for storing variables too numerous to put in a tree.
22//
23// --- Adapted to look for lambdas as well, using code from
24// AliAnalysisTaskCheckPerformanceStrange.cxx
25//
26// --- Algorithm Description
27// 1. Loop over primaries in stack to acquire generated charged Xi
9a8f3aee 28// 2. Loop over stack to find Cascades, fill TH3Fs "PrimRawPt"s for Efficiency
76029adc 29// 3. Perform Physics Selection
30// 4. Perform Primary Vertex |z|<10cm selection
31// 5. Perform Primary Vertex NoTPCOnly vertexing selection (>0 contrib.)
32// 6. Perform Pileup Rejection
33// 7. Analysis Loops:
34// 7a. Fill TH3Fs "PrimAnalysisPt" for control purposes only
76029adc 35//
36// Please Report Any Bugs!
37//
38// --- David Dobrigkeit Chinellato
39// (david.chinellato@gmail.com)
40//
41// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
42
43class TTree;
44class TParticle;
45class TVector3;
46
47//class AliMCEventHandler;
48//class AliMCEvent;
49//class AliStack;
50
51class AliESDVertex;
52class AliAODVertex;
53class AliESDv0;
54class AliAODv0;
55
56#include <Riostream.h>
57#include "TList.h"
58#include "TH1.h"
59#include "TH2.h"
60#include "TH3.h"
61#include "TFile.h"
62#include "THnSparse.h"
63#include "TVector3.h"
64#include "TCanvas.h"
65#include "TMath.h"
66#include "TLegend.h"
67#include "AliLog.h"
68
69#include "AliESDEvent.h"
70#include "AliAODEvent.h"
71#include "AliV0vertexer.h"
72#include "AliCascadeVertexer.h"
73#include "AliESDpid.h"
74#include "AliESDtrack.h"
75#include "AliESDtrackCuts.h"
76#include "AliInputEventHandler.h"
77#include "AliAnalysisManager.h"
78#include "AliMCEventHandler.h"
79#include "AliMCEvent.h"
80#include "AliStack.h"
81
82#include "AliCFContainer.h"
83#include "AliMultiplicity.h"
84#include "AliAODMCParticle.h"
85#include "AliESDcascade.h"
86#include "AliAODcascade.h"
87#include "AliESDUtils.h"
88#include "AliGenEventHeader.h"
89
90#include "AliAnalysisTaskExtractPerformanceCascade.h"
91
92using std::cout;
93using std::endl;
94
95ClassImp(AliAnalysisTaskExtractPerformanceCascade)
96
97AliAnalysisTaskExtractPerformanceCascade::AliAnalysisTaskExtractPerformanceCascade()
98 : AliAnalysisTaskSE(), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0),
99 fkIsNuclear ( kFALSE ),
100 fkLowEnergyPP ( kFALSE ),
101
102//------------------------------------------------
103// HISTOGRAMS
104// --- Filled on an Event-by-event basis
105//------------------------------------------------
106 fHistV0MultiplicityBeforeTrigSel(0),
107 fHistV0MultiplicityForTrigEvt(0),
108 fHistV0MultiplicityForSelEvt(0),
109 fHistV0MultiplicityForSelEvtNoTPCOnly(0),
110 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
111 fHistMultiplicityBeforeTrigSel(0),
112 fHistMultiplicityForTrigEvt(0),
113 fHistMultiplicity(0),
114 fHistMultiplicityNoTPCOnly(0),
115 fHistMultiplicityNoTPCOnlyNoPileup(0),
116
117//------------------------------------------------
118// PARTICLE HISTOGRAMS
119// --- Filled on a Particle-by-Particle basis
120//------------------------------------------------
121 f3dHistGenPtVsYVsMultXiMinus(0),
122 f3dHistGenPtVsYVsMultXiPlus(0),
123 f3dHistGenPtVsYVsMultOmegaMinus(0),
124 f3dHistGenPtVsYVsMultOmegaPlus(0),
125 f3dHistGenSelectedPtVsYVsMultXiMinus(0),
126 f3dHistGenSelectedPtVsYVsMultXiPlus(0),
127 f3dHistGenSelectedPtVsYVsMultOmegaMinus(0),
128 f3dHistGenSelectedPtVsYVsMultOmegaPlus(0),
129 fHistPVx(0),
130 fHistPVy(0),
131 fHistPVz(0),
132 fHistPVxAnalysis(0),
133 fHistPVyAnalysis(0),
134 fHistPVzAnalysis(0)
135{
136 // Dummy Constructor
137}
138
139AliAnalysisTaskExtractPerformanceCascade::AliAnalysisTaskExtractPerformanceCascade(const char *name)
140 : AliAnalysisTaskSE(name), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0),
141 fkIsNuclear ( kFALSE ),
142 fkLowEnergyPP ( kFALSE ),
143
144//------------------------------------------------
145// HISTOGRAMS
146// --- Filled on an Event-by-event basis
147//------------------------------------------------
148 fHistV0MultiplicityBeforeTrigSel(0),
149 fHistV0MultiplicityForTrigEvt(0),
150 fHistV0MultiplicityForSelEvt(0),
151 fHistV0MultiplicityForSelEvtNoTPCOnly(0),
152 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
153 fHistMultiplicityBeforeTrigSel(0),
154 fHistMultiplicityForTrigEvt(0),
155 fHistMultiplicity(0),
156 fHistMultiplicityNoTPCOnly(0),
157 fHistMultiplicityNoTPCOnlyNoPileup(0),
158
159
160//------------------------------------------------
161// PARTICLE HISTOGRAMS
162// --- Filled on a Particle-by-Particle basis
163//------------------------------------------------
164 f3dHistGenPtVsYVsMultXiMinus(0),
165 f3dHistGenPtVsYVsMultXiPlus(0),
166 f3dHistGenPtVsYVsMultOmegaMinus(0),
167 f3dHistGenPtVsYVsMultOmegaPlus(0),
168 f3dHistGenSelectedPtVsYVsMultXiMinus(0),
169 f3dHistGenSelectedPtVsYVsMultXiPlus(0),
170 f3dHistGenSelectedPtVsYVsMultOmegaMinus(0),
171 f3dHistGenSelectedPtVsYVsMultOmegaPlus(0),
172 fHistPVx(0),
173 fHistPVy(0),
174 fHistPVz(0),
175 fHistPVxAnalysis(0),
176 fHistPVyAnalysis(0),
177 fHistPVzAnalysis(0)
178{
179 // Constructor
d8841e95 180
181 //Set Variables for re-running the cascade vertexers (as done for MS paper)
182
183 // New Loose : 1st step for the 7 TeV pp analysis
184
185 fV0Sels[0] = 33. ; // max allowed chi2
186 fV0Sels[1] = 0.02; // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
187 fV0Sels[2] = 0.02; // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
188 fV0Sels[3] = 2.0 ; // max allowed DCA between the daughter tracks (LHC09a4 : 0.5)
189 fV0Sels[4] = 0.95; // min allowed cosine of V0's pointing angle (LHC09a4 : 0.99)
190 fV0Sels[5] = 1.0 ; // min radius of the fiducial volume (LHC09a4 : 0.2)
191 fV0Sels[6] = 100. ; // max radius of the fiducial volume (LHC09a4 : 100.0)
192
193 fCascSels[0] = 33. ; // max allowed chi2 (same as PDC07)
194 fCascSels[1] = 0.05 ; // min allowed V0 impact parameter (PDC07 : 0.05 / LHC09a4 : 0.025 )
195 fCascSels[2] = 0.010; // "window" around the Lambda mass (PDC07 : 0.008 / LHC09a4 : 0.010 )
196 fCascSels[3] = 0.03 ; // min allowed bachelor's impact parameter (PDC07 : 0.035 / LHC09a4 : 0.025 )
197 fCascSels[4] = 2.0 ; // max allowed DCA between the V0 and the bachelor (PDC07 : 0.1 / LHC09a4 : 0.2 )
198 fCascSels[5] = 0.95 ; // min allowed cosine of the cascade pointing angle (PDC07 : 0.9985 / LHC09a4 : 0.998 )
199 fCascSels[6] = 0.4 ; // min radius of the fiducial volume (PDC07 : 0.9 / LHC09a4 : 0.2 )
200 fCascSels[7] = 100. ; // max radius of the fiducial volume (PDC07 : 100 / LHC09a4 : 100 )
201
76029adc 202 // Output slot #0 writes into a TList container (Cascade)
203 DefineOutput(1, TList::Class());
204 DefineOutput(2, TTree::Class());
205}
206
207
208AliAnalysisTaskExtractPerformanceCascade::~AliAnalysisTaskExtractPerformanceCascade()
209{
210//------------------------------------------------
211// DESTRUCTOR
212//------------------------------------------------
213
214 if (fListHist){
215 delete fListHist;
216 fListHist = 0x0;
217 }
218 if (fTreeCascade){
219 delete fTreeCascade;
220 fTreeCascade = 0x0;
221 }
222 //cleanup esd track cuts object too...
223 if (fESDtrackCuts){
224 delete fESDtrackCuts;
225 fESDtrackCuts = 0x0;
226 }
227
228}
229
230//________________________________________________________________________
231void AliAnalysisTaskExtractPerformanceCascade::UserCreateOutputObjects()
232{
233 OpenFile(2);
234 // Called once
235
236//------------------------------------------------
237
238 fTreeCascade = new TTree("fTreeCascade","CascadeCandidates");
239
240//------------------------------------------------
241// fTreeCascade Branch definitions - Cascade Tree
242//------------------------------------------------
243
244//------------------------------------------------
245// fTreeCascade Branch definitions
246//------------------------------------------------
247
248//-----------BASIC-INFO---------------------------
249/* 1*/ fTreeCascade->Branch("fTreeCascVarCharge",&fTreeCascVarCharge,"fTreeCascVarCharge/I");
250/* 2*/ fTreeCascade->Branch("fTreeCascVarMassAsXi",&fTreeCascVarMassAsXi,"fTreeCascVarMassAsXi/F");
251/* 3*/ fTreeCascade->Branch("fTreeCascVarMassAsOmega",&fTreeCascVarMassAsOmega,"fTreeCascVarMassAsOmega/F");
252/* 4*/ fTreeCascade->Branch("fTreeCascVarPt",&fTreeCascVarPt,"fTreeCascVarPt/F");
253/* 5*/ fTreeCascade->Branch("fTreeCascVarPtMC",&fTreeCascVarPtMC,"fTreeCascVarPtMC/F");
254/* 6*/ fTreeCascade->Branch("fTreeCascVarRapXi",&fTreeCascVarRapXi,"fTreeCascVarRapXi/F");
255/* 7*/ fTreeCascade->Branch("fTreeCascVarRapMC",&fTreeCascVarRapMC,"fTreeCascVarRapMC/F");
256/* 8*/ fTreeCascade->Branch("fTreeCascVarRapOmega",&fTreeCascVarRapOmega,"fTreeCascVarRapOmega/F");
257/* 9*/ fTreeCascade->Branch("fTreeCascVarNegEta",&fTreeCascVarNegEta,"fTreeCascVarNegEta/F");
258/*10*/ fTreeCascade->Branch("fTreeCascVarPosEta",&fTreeCascVarPosEta,"fTreeCascVarPosEta/F");
259/*11*/ fTreeCascade->Branch("fTreeCascVarBachEta",&fTreeCascVarBachEta,"fTreeCascVarBachEta/F");
260//-----------INFO-FOR-CUTS------------------------
261/*12*/ fTreeCascade->Branch("fTreeCascVarDCACascDaughters",&fTreeCascVarDCACascDaughters,"fTreeCascVarDCACascDaughters/F");
262/*13*/ fTreeCascade->Branch("fTreeCascVarDCABachToPrimVtx",&fTreeCascVarDCABachToPrimVtx,"fTreeCascVarDCABachToPrimVtx/F");
263/*14*/ fTreeCascade->Branch("fTreeCascVarDCAV0Daughters",&fTreeCascVarDCAV0Daughters,"fTreeCascVarDCAV0Daughters/F");
264/*15*/ fTreeCascade->Branch("fTreeCascVarDCAV0ToPrimVtx",&fTreeCascVarDCAV0ToPrimVtx,"fTreeCascVarDCAV0ToPrimVtx/F");
265/*16*/ fTreeCascade->Branch("fTreeCascVarDCAPosToPrimVtx",&fTreeCascVarDCAPosToPrimVtx,"fTreeCascVarDCAPosToPrimVtx/F");
266/*17*/ fTreeCascade->Branch("fTreeCascVarDCANegToPrimVtx",&fTreeCascVarDCANegToPrimVtx,"fTreeCascVarDCANegToPrimVtx/F");
267/*18*/ fTreeCascade->Branch("fTreeCascVarCascCosPointingAngle",&fTreeCascVarCascCosPointingAngle,"fTreeCascVarCascCosPointingAngle/F");
268/*19*/ fTreeCascade->Branch("fTreeCascVarCascRadius",&fTreeCascVarCascRadius,"fTreeCascVarCascRadius/F");
269/*20*/ fTreeCascade->Branch("fTreeCascVarV0Mass",&fTreeCascVarV0Mass,"fTreeCascVarV0Mass/F");
270/*21*/ fTreeCascade->Branch("fTreeCascVarV0CosPointingAngle",&fTreeCascVarV0CosPointingAngle,"fTreeCascVarV0CosPointingAngle/F");
271/*22*/ fTreeCascade->Branch("fTreeCascVarV0Radius",&fTreeCascVarV0Radius,"fTreeCascVarV0Radius/F");
272/*23*/ fTreeCascade->Branch("fTreeCascVarLeastNbrClusters",&fTreeCascVarLeastNbrClusters,"fTreeCascVarLeastNbrClusters/I");
273//-----------MULTIPLICITY-INFO--------------------
274/*24*/ fTreeCascade->Branch("fTreeCascVarMultiplicity",&fTreeCascVarMultiplicity,"fTreeCascVarMultiplicity/I");
275//-----------DECAY-LENGTH-INFO--------------------
276/*25*/ fTreeCascade->Branch("fTreeCascVarDistOverTotMom",&fTreeCascVarDistOverTotMom,"fTreeCascVarDistOverTotMom/F");
277//-----------MC-PID-------------------------------
278/*26*/ fTreeCascade->Branch("fTreeCascVarPID",&fTreeCascVarPID,"fTreeCascVarPID/I");
279/*27*/ fTreeCascade->Branch("fTreeCascVarPIDBachelor",&fTreeCascVarPIDBachelor,"fTreeCascVarPIDBachelor/I");
280/*28*/ fTreeCascade->Branch("fTreeCascVarPIDNegative",&fTreeCascVarPIDNegative,"fTreeCascVarPIDNegative/I");
281/*29*/ fTreeCascade->Branch("fTreeCascVarPIDPositive",&fTreeCascVarPIDPositive,"fTreeCascVarPIDPositive/I");
282/*30*/ fTreeCascade->Branch("fTreeCascVarPosTransMom",&fTreeCascVarPosTransMom,"fTreeCascVarPosTransMom/F");
283/*31*/ fTreeCascade->Branch("fTreeCascVarNegTransMom",&fTreeCascVarNegTransMom,"fTreeCascVarNegTransMom/F");
284/*32*/ fTreeCascade->Branch("fTreeCascVarPosTransMomMC",&fTreeCascVarPosTransMomMC,"fTreeCascVarPosTransMomMC/F");
285/*33*/ fTreeCascade->Branch("fTreeCascVarNegTransMomMC",&fTreeCascVarNegTransMomMC,"fTreeCascVarNegTransMomMC/F");
286//------------------------------------------------
287/*34*/ fTreeCascade->Branch("fTreeCascVarNegNSigmaPion",&fTreeCascVarNegNSigmaPion,"fTreeCascVarNegNSigmaPion/F");
288/*35*/ fTreeCascade->Branch("fTreeCascVarNegNSigmaProton",&fTreeCascVarNegNSigmaProton,"fTreeCascVarNegNSigmaProton/F");
289/*36*/ fTreeCascade->Branch("fTreeCascVarPosNSigmaPion",&fTreeCascVarPosNSigmaPion,"fTreeCascVarPosNSigmaPion/F");
290/*37*/ fTreeCascade->Branch("fTreeCascVarPosNSigmaProton",&fTreeCascVarPosNSigmaProton,"fTreeCascVarPosNSigmaProton/F");
291/*38*/ fTreeCascade->Branch("fTreeCascVarBachNSigmaPion",&fTreeCascVarBachNSigmaPion,"fTreeCascVarBachNSigmaPion/F");
292/*39*/ fTreeCascade->Branch("fTreeCascVarBachNSigmaKaon",&fTreeCascVarBachNSigmaKaon,"fTreeCascVarBachNSigmaKaon/F");
293
294//------------------------------------------------
295// Particle Identification Setup
296//------------------------------------------------
297
298 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
299 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
300 fPIDResponse = inputHandler->GetPIDResponse();
301
302// Multiplicity
303
304 if(! fESDtrackCuts ){
305 fESDtrackCuts = new AliESDtrackCuts();
306 }
307
308//------------------------------------------------
309// V0 Multiplicity Histograms
310//------------------------------------------------
311
312 // Create histograms
313 OpenFile(1);
314 fListHist = new TList();
315 fListHist->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
316
317
318 if(! fHistV0MultiplicityBeforeTrigSel) {
319 fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel",
320 "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events",
321 25, 0, 25);
322 fListHist->Add(fHistV0MultiplicityBeforeTrigSel);
323 }
324
325 if(! fHistV0MultiplicityForTrigEvt) {
326 fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt",
327 "V0s per event (for triggered evt);Nbr of V0s/Evt;Events",
328 25, 0, 25);
329 fListHist->Add(fHistV0MultiplicityForTrigEvt);
330 }
331
332 if(! fHistV0MultiplicityForSelEvt) {
333 fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt",
334 "V0s per event;Nbr of V0s/Evt;Events",
335 25, 0, 25);
336 fListHist->Add(fHistV0MultiplicityForSelEvt);
337 }
338
339 if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
340 fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly",
341 "V0s per event;Nbr of V0s/Evt;Events",
342 25, 0, 25);
343 fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
344 }
345 if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
346 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup",
347 "V0s per event;Nbr of V0s/Evt;Events",
348 25, 0, 25);
349 fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
350 }
351
352//------------------------------------------------
353// Track Multiplicity Histograms
354//------------------------------------------------
355
356 if(! fHistMultiplicityBeforeTrigSel) {
357 fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel",
358 "Tracks per event;Nbr of Tracks;Events",
359 200, 0, 200);
360 fListHist->Add(fHistMultiplicityBeforeTrigSel);
361 }
362 if(! fHistMultiplicityForTrigEvt) {
363 fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt",
364 "Tracks per event;Nbr of Tracks;Events",
365 200, 0, 200);
366 fListHist->Add(fHistMultiplicityForTrigEvt);
367 }
368 if(! fHistMultiplicity) {
369 fHistMultiplicity = new TH1F("fHistMultiplicity",
370 "Tracks per event;Nbr of Tracks;Events",
371 200, 0, 200);
372 fListHist->Add(fHistMultiplicity);
373 }
374 if(! fHistMultiplicityNoTPCOnly) {
375 fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly",
376 "Tracks per event;Nbr of Tracks;Events",
377 200, 0, 200);
378 fListHist->Add(fHistMultiplicityNoTPCOnly);
379 }
380 if(! fHistMultiplicityNoTPCOnlyNoPileup) {
381 fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup",
382 "Tracks per event;Nbr of Tracks;Events",
383 200, 0, 200);
384 fListHist->Add(fHistMultiplicityNoTPCOnlyNoPileup);
385 }
386
387//------------------------------------------------
388// Generated Particle Histograms
389//------------------------------------------------
390
391 Int_t lCustomNBins = 200;
392 Double_t lCustomPtUpperLimit = 20;
393 Int_t lCustomNBinsMultiplicity = 100;
394
395//----------------------------------
396// Raw Generated (Pre-physics-selection)
397//----------------------------------
398
399//--------------------------------------------------------------------------------------
400//--- 3D Histo (Pt, Y, Multiplicity) for generated XiMinus/Plus, all generated
401
402 if(! f3dHistGenPtVsYVsMultXiMinus) {
403 f3dHistGenPtVsYVsMultXiMinus = new TH3F( "f3dHistGenPtVsYVsMultXiMinus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
404 fListHist->Add(f3dHistGenPtVsYVsMultXiMinus);
405 }
406 if(! f3dHistGenPtVsYVsMultXiPlus) {
407 f3dHistGenPtVsYVsMultXiPlus = new TH3F( "f3dHistGenPtVsYVsMultXiPlus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
408 fListHist->Add(f3dHistGenPtVsYVsMultXiPlus);
409 }
410//--- 3D Histo (Pt, Y, Multiplicity) for generated OmegaMinus/Plus
411
412 if(! f3dHistGenPtVsYVsMultOmegaMinus) {
413 f3dHistGenPtVsYVsMultOmegaMinus = new TH3F( "f3dHistGenPtVsYVsMultOmegaMinus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
414 fListHist->Add(f3dHistGenPtVsYVsMultOmegaMinus);
415 }
416 if(! f3dHistGenPtVsYVsMultOmegaPlus) {
417 f3dHistGenPtVsYVsMultOmegaPlus = new TH3F( "f3dHistGenPtVsYVsMultOmegaPlus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
418 fListHist->Add(f3dHistGenPtVsYVsMultOmegaPlus);
419 }
420
421//--------------------------------------------------------------------------------------
422//--- 3D Histo (Pt, Y, Multiplicity) for generated XiMinus/Plus, at selected analysis evts
423
424 if(! f3dHistGenSelectedPtVsYVsMultXiMinus) {
425 f3dHistGenSelectedPtVsYVsMultXiMinus = new TH3F( "f3dHistGenSelectedPtVsYVsMultXiMinus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
426 fListHist->Add(f3dHistGenSelectedPtVsYVsMultXiMinus);
427 }
428 if(! f3dHistGenSelectedPtVsYVsMultXiPlus) {
429 f3dHistGenSelectedPtVsYVsMultXiPlus = new TH3F( "f3dHistGenSelectedPtVsYVsMultXiPlus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
430 fListHist->Add(f3dHistGenSelectedPtVsYVsMultXiPlus);
431 }
432//--- 3D Histo (Pt, Y, Multiplicity) for generated OmegaMinus/Plus
433
434 if(! f3dHistGenSelectedPtVsYVsMultOmegaMinus) {
435 f3dHistGenSelectedPtVsYVsMultOmegaMinus = new TH3F( "f3dHistGenSelectedPtVsYVsMultOmegaMinus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
436 fListHist->Add(f3dHistGenSelectedPtVsYVsMultOmegaMinus);
437 }
438 if(! f3dHistGenSelectedPtVsYVsMultOmegaPlus) {
439 f3dHistGenSelectedPtVsYVsMultOmegaPlus = new TH3F( "f3dHistGenSelectedPtVsYVsMultOmegaPlus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
440 fListHist->Add(f3dHistGenSelectedPtVsYVsMultOmegaPlus);
441 }
442
443//----------------------------------
444// Primary Vertex Position Histos
445//----------------------------------
446
447 if(! fHistPVx) {
448 fHistPVx = new TH1F("fHistPVx",
449 "PV x position;Nbr of Evts;x",
450 2000, -0.5, 0.5);
451 fListHist->Add(fHistPVx);
452 }
453 if(! fHistPVy) {
454 fHistPVy = new TH1F("fHistPVy",
455 "PV y position;Nbr of Evts;y",
456 2000, -0.5, 0.5);
457 fListHist->Add(fHistPVy);
458 }
459 if(! fHistPVz) {
460 fHistPVz = new TH1F("fHistPVz",
461 "PV z position;Nbr of Evts;z",
462 400, -20, 20);
463 fListHist->Add(fHistPVz);
464 }
465
466 if(! fHistPVxAnalysis) {
467 fHistPVxAnalysis = new TH1F("fHistPVxAnalysis",
468 "PV x position;Nbr of Evts;x",
469 2000, -0.5, 0.5);
470 fListHist->Add(fHistPVxAnalysis);
471 }
472 if(! fHistPVyAnalysis) {
473 fHistPVyAnalysis = new TH1F("fHistPVyAnalysis",
474 "PV y position;Nbr of Evts;y",
475 2000, -0.5, 0.5);
476 fListHist->Add(fHistPVyAnalysis);
477 }
478 if(! fHistPVzAnalysis) {
479 fHistPVzAnalysis = new TH1F("fHistPVzAnalysis",
480 "PV z position;Nbr of Evts;z",
481 400, -20, 20);
482 fListHist->Add(fHistPVzAnalysis);
483 }
484
485 //List of Histograms: Normal
486 PostData(1, fListHist);
487
488 //TTree Object: Saved to base directory. Should cache to disk while saving.
489 //(Important to avoid excessive memory usage, particularly when merging)
490 PostData(2, fTreeCascade);
491
492}// end UserCreateOutputObjects
493
494
495//________________________________________________________________________
496void AliAnalysisTaskExtractPerformanceCascade::UserExec(Option_t *)
497{
498 // Main loop
499 // Called for each event
500
501 AliESDEvent *lESDevent = 0x0;
502 AliMCEvent *lMCevent = 0x0;
503 AliStack *lMCstack = 0x0;
504
505 Int_t lNumberOfV0s = -1;
506 Double_t lTrkgPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
507 Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
508 Double_t lMagneticField = -10.;
509
510 // Connect to the InputEvent
511 // After these lines, we should have an ESD/AOD event + the number of V0s in it.
512
513 // Appropriate for ESD analysis!
514
515 lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
516 if (!lESDevent) {
517 AliWarning("ERROR: lESDevent not available \n");
518 return;
519 }
520
521/* --- Acquisition of exact event ID
522 fTreeVariableRunNumber = lESDevent->GetRunNumber();
523 fTreeVariableEventNumber =
524 ( ( ((ULong64_t)lESDevent->GetPeriodNumber() ) << 36 ) |
525 ( ((ULong64_t)lESDevent->GetOrbitNumber () ) << 12 ) |
526 ((ULong64_t)lESDevent->GetBunchCrossNumber() ) );
527*/
528 lMCevent = MCEvent();
529 if (!lMCevent) {
530 Printf("ERROR: Could not retrieve MC event \n");
531 cout << "Name of the file with pb :" << fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
532 return;
533 }
534
535 lMCstack = lMCevent->Stack();
536 if (!lMCstack) {
537 Printf("ERROR: Could not retrieve MC stack \n");
538 cout << "Name of the file with pb :" << fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
539 return;
540 }
541 TArrayF mcPrimaryVtx;
542 AliGenEventHeader* mcHeader=lMCevent->GenEventHeader();
543 if(!mcHeader) return;
544 mcHeader->PrimaryVertex(mcPrimaryVtx);
545
546//------------------------------------------------
547// Multiplicity Information Acquistion
548//------------------------------------------------
549
550 //REVISED multiplicity estimator after 'multiplicity day' (2011)
551 Int_t lMultiplicity = -100;
552
553 //testing purposes
554 if(fkIsNuclear == kFALSE) lMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
555
556 //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
557 //---> Warning: Experimental
558 if(fkIsNuclear == kTRUE){
559 AliCentrality* centrality;
560 centrality = lESDevent->GetCentrality();
561 lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( "V0M" ) ) );
562 if (centrality->GetQuality()>1) {
563 PostData(1, fListHist);
564 PostData(2, fTreeCascade);
565 return;
566 }
567 }
568
569 //Set variable for filling tree afterwards!
570 //---> pp case......: GetReferenceMultiplicity
571 //---> Pb-Pb case...: Centrality by V0M
572
573 fTreeCascVarMultiplicity = lMultiplicity;
574
575 fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
576 fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
577
578//------------------------------------------------
579// MC Information Acquistion
580//------------------------------------------------
581
582 Int_t iNumberOfPrimaries = -1;
583 iNumberOfPrimaries = lMCstack->GetNprimary();
584 if(iNumberOfPrimaries < 1) return;
585
586//------------------------------------------------
587// Variable Definition
588//------------------------------------------------
589
590 Int_t lNbMCPrimary = 0;
591
76029adc 592 Double_t lPtCurrentPart = 0;
593
594 //Int_t lComeFromSigma = 0;
595
596 // current mc particle 's mother
597 //Int_t iCurrentMother = 0;
598 lNbMCPrimary = lMCstack->GetNprimary();
599
600//------------------------------------------------
601// Pre-Physics Selection
602//------------------------------------------------
603
604//----- Loop on primary Xi, Omega --------------------------------------------------------------
605 for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++)
606 {// This is the begining of the loop on primaries
607
608 TParticle* lCurrentParticlePrimary = 0x0;
609 lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
610 if(!lCurrentParticlePrimary){
611 Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
612 continue;
613 }
614 if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 || TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3334 ) {
615 Double_t lRapXiMCPrimary = -100;
616 if( (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) != 0 ) {
617 if ( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) !=0 ){
618 lRapXiMCPrimary = 0.5*TMath::Log( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) );
619 }
620 }
621
622 //=================================================================================
623 // Xi Histograms
624 if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){
625 lPtCurrentPart = lCurrentParticlePrimary->Pt();
626 f3dHistGenPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
627 }
628 if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){
629 lPtCurrentPart = lCurrentParticlePrimary->Pt();
630 f3dHistGenPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
631 }
632 // Omega Histograms
633 if( lCurrentParticlePrimary->GetPdgCode() == 3334 ){
634 lPtCurrentPart = lCurrentParticlePrimary->Pt();
635 f3dHistGenPtVsYVsMultOmegaMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
636 }
637 if( lCurrentParticlePrimary->GetPdgCode() == -3334 ){
638 lPtCurrentPart = lCurrentParticlePrimary->Pt();
639 f3dHistGenPtVsYVsMultOmegaPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
640 }
641 }
642 }
643//----- End Loop on primary Xi, Omega ----------------------------------------------------------
644
645// ---> Set Variables to Zero again
646// ---> Variable Definition
647
76029adc 648 lPtCurrentPart = 0;
649
650//------------------------------------------------
651// Physics Selection
652//------------------------------------------------
653
654 UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
655 Bool_t isSelected = 0;
656 isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
657
658 //pp at 2.76TeV: special case, ignore FastOnly
659 if ( (fkLowEnergyPP == kTRUE) && (maskIsSelected& AliVEvent::kFastOnly) ){
660 PostData(1, fListHist);
661 PostData(2, fTreeCascade);
662 return;
663 }
664 //Standard Min-Bias Selection
665 if ( ! isSelected ) {
666 PostData(1, fListHist);
667 PostData(2, fTreeCascade);
668 return;
669 }
670
d8841e95 671//------------------------------------------------
672// Rerun cascade vertexer!
673//------------------------------------------------
674
675 lESDevent->ResetCascades();
676 lESDevent->ResetV0s();
677
678 AliV0vertexer lV0vtxer;
679 AliCascadeVertexer lCascVtxer;
680
681 lV0vtxer.SetDefaultCuts(fV0Sels);
682 lCascVtxer.SetDefaultCuts(fCascSels);
683
684 lV0vtxer.Tracks2V0vertices(lESDevent);
685 lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
686
76029adc 687//------------------------------------------------
688// After Trigger Selection
689//------------------------------------------------
690
691 lNumberOfV0s = lESDevent->GetNumberOfV0s();
692
693 //Set variable for filling tree afterwards!
694 fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
695 fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
696
697//------------------------------------------------
698// Getting: Primary Vertex + MagField Info
699//------------------------------------------------
700
701 const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
702 // get the vtx stored in ESD found with tracks
703 lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
704
705 const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
706 // get the best primary vertex available for the event
707 // As done in AliCascadeVertexer, we keep the one which is the best one available.
708 // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
709 // This one will be used for next calculations (DCA essentially)
710 lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
711
712 Double_t lPrimaryVtxPosition[3];
713 const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
714 lPrimaryVtxPosition[0] = primaryVtx->GetX();
715 lPrimaryVtxPosition[1] = primaryVtx->GetY();
716 lPrimaryVtxPosition[2] = primaryVtx->GetZ();
717 fHistPVx->Fill( lPrimaryVtxPosition[0] );
718 fHistPVy->Fill( lPrimaryVtxPosition[1] );
719 fHistPVz->Fill( lPrimaryVtxPosition[2] );
720
721//------------------------------------------------
722// Primary Vertex Z position: SKIP
723//------------------------------------------------
724
725 if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) {
726 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
727 PostData(1, fListHist);
728 PostData(2, fTreeCascade);
729 return;
730 }
731
732 lMagneticField = lESDevent->GetMagneticField( );
733 fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
734 fHistMultiplicity->Fill(lMultiplicity);
735
736//------------------------------------------------
737// SKIP: Events with well-established PVtx
738//------------------------------------------------
739
740 const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
741 const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
742 if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
743 AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
744 PostData(1, fListHist);
745 PostData(2, fTreeCascade);
746 return;
747 }
748 fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
749 fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
750
751//------------------------------------------------
752// Pileup Rejection Studies
753//------------------------------------------------
754
755 // FIXME : quality selection regarding pile-up rejection
756 if(lESDevent->IsPileupFromSPD() && !fkIsNuclear ){// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5. -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
757 AliWarning("Pb / This is tagged as Pileup from SPD... return !");
758 PostData(1, fListHist);
759 PostData(2, fTreeCascade);
760 return;
761 }
762 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
763 fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
764
765 //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
766 if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() ) ){
767 fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
768 fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
769 fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
770 }
771
772//------------------------------------------------
773// stack loop starts here
774//------------------------------------------------
775
776//----- Loop on primary Xi, Omega --------------------------------------------------------------
777 for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++)
778 {// This is the begining of the loop on primaries
779
780 TParticle* lCurrentParticlePrimary = 0x0;
781 lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
782 if(!lCurrentParticlePrimary){
783 Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
784 continue;
785 }
786 if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 || TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3334 ) {
787 Double_t lRapXiMCPrimary = -100;
788 if( (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) != 0 ) {
789 if ( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) !=0 ){
790 lRapXiMCPrimary = 0.5*TMath::Log( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) );
791 }
792 }
793
794 //=================================================================================
795 // Xi Histograms
796 if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){
797 lPtCurrentPart = lCurrentParticlePrimary->Pt();
798 f3dHistGenSelectedPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
799 }
800 if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){
801 lPtCurrentPart = lCurrentParticlePrimary->Pt();
802 f3dHistGenSelectedPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
803 }
804 // Omega Histograms
805 if( lCurrentParticlePrimary->GetPdgCode() == 3334 ){
806 lPtCurrentPart = lCurrentParticlePrimary->Pt();
807 f3dHistGenSelectedPtVsYVsMultOmegaMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
808 }
809 if( lCurrentParticlePrimary->GetPdgCode() == -3334 ){
810 lPtCurrentPart = lCurrentParticlePrimary->Pt();
811 f3dHistGenSelectedPtVsYVsMultOmegaPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
812 }
813 }
814 }
815//----- End Loop on primary Xi, Omega ----------------------------------------------------------
816
817//------------------------------------------------
818// MAIN CASCADE LOOP STARTS HERE
819//------------------------------------------------
820// Code Credit: Antonin Maire (thanks^100)
821// ---> This is an adaptation
822
823 Long_t ncascades = 0;
824 ncascades = lESDevent->GetNumberOfCascades();
825
826
827 for (Int_t iXi = 0; iXi < ncascades; iXi++){
828 //------------------------------------------------
829 // Initializations
830 //------------------------------------------------
831 //Double_t lTrkgPrimaryVtxRadius3D = -500.0;
832 //Double_t lBestPrimaryVtxRadius3D = -500.0;
833
834 // - 1st part of initialisation : variables needed to store AliESDCascade data members
835 Double_t lEffMassXi = 0. ;
836 //Double_t lChi2Xi = -1. ;
837 Double_t lDcaXiDaughters = -1. ;
838 Double_t lXiCosineOfPointingAngle = -1. ;
839 Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
840 Double_t lXiRadius = -1000. ;
841
842 // - 2nd part of initialisation : Nbr of clusters within TPC for the 3 daughter cascade tracks
843 Int_t lPosTPCClusters = -1; // For ESD only ...//FIXME : wait for availability in AOD
844 Int_t lNegTPCClusters = -1; // For ESD only ...
845 Int_t lBachTPCClusters = -1; // For ESD only ...
846
847 // - 3rd part of initialisation : about V0 part in cascades
848 Double_t lInvMassLambdaAsCascDghter = 0.;
849 //Double_t lV0Chi2Xi = -1. ;
850 Double_t lDcaV0DaughtersXi = -1.;
851
852 Double_t lDcaBachToPrimVertexXi = -1., lDcaV0ToPrimVertexXi = -1.;
853 Double_t lDcaPosToPrimVertexXi = -1.;
854 Double_t lDcaNegToPrimVertexXi = -1.;
855 Double_t lV0CosineOfPointingAngleXi = -1. ;
856 Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
857 Double_t lV0RadiusXi = -1000.0;
858 Double_t lV0quality = 0.;
859
860 // - 4th part of initialisation : Effective masses
861 Double_t lInvMassXiMinus = 0.;
862 Double_t lInvMassXiPlus = 0.;
863 Double_t lInvMassOmegaMinus = 0.;
864 Double_t lInvMassOmegaPlus = 0.;
865
866 // - 6th part of initialisation : extra info for QA
867 Double_t lXiMomX = 0. , lXiMomY = 0., lXiMomZ = 0.;
868 Double_t lXiTransvMom = 0. ;
869 Double_t lXiTransvMomMC= 0. ;
870 Double_t lXiTotMom = 0. ;
871
872 Double_t lBachMomX = 0., lBachMomY = 0., lBachMomZ = 0.;
873 //Double_t lBachTransvMom = 0.;
874 //Double_t lBachTotMom = 0.;
875
876 fTreeCascVarNegNSigmaPion = -100;
877 fTreeCascVarNegNSigmaProton = -100;
878 fTreeCascVarPosNSigmaPion = -100;
879 fTreeCascVarPosNSigmaProton = -100;
880 fTreeCascVarBachNSigmaPion = -100;
881 fTreeCascVarBachNSigmaKaon = -100;
882
883 Short_t lChargeXi = -2;
884 //Double_t lV0toXiCosineOfPointingAngle = 0. ;
885
886 Double_t lRapXi = -20.0, lRapOmega = -20.0, lRapMC = -20.0; // lEta = -20.0, lTheta = 360., lPhi = 720. ;
887 //Double_t lAlphaXi = -200., lPtArmXi = -200.0;
888
889 // -------------------------------------
890 // II.ESD - Calculation Part dedicated to Xi vertices (ESD)
891
892 AliESDcascade *xi = lESDevent->GetCascade(iXi);
893 if (!xi) continue;
894
895
896 // - II.Step 1 : around primary vertex
897 //-------------
898 //lTrkgPrimaryVtxRadius3D = TMath::Sqrt( lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
899 // lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
900 // lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
901
902 //lBestPrimaryVtxRadius3D = TMath::Sqrt( lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
903 // lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
904 // lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
905
906 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)
907 //-------------
908 lV0quality = 0.;
909 xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
910
911 lEffMassXi = xi->GetEffMassXi();
912 //lChi2Xi = xi->GetChi2Xi();
913 lDcaXiDaughters = xi->GetDcaXiDaughters();
914 lXiCosineOfPointingAngle = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
915 lBestPrimaryVtxPos[1],
916 lBestPrimaryVtxPos[2] );
917 // Take care : the best available vertex should be used (like in AliCascadeVertexer)
918
919 xi->GetXYZcascade( lPosXi[0], lPosXi[1], lPosXi[2] );
920 lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0] + lPosXi[1]*lPosXi[1] );
921
922 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
923 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
924 //-------------
925
926 UInt_t lIdxPosXi = (UInt_t) TMath::Abs( xi->GetPindex() );
927 UInt_t lIdxNegXi = (UInt_t) TMath::Abs( xi->GetNindex() );
928 UInt_t lBachIdx = (UInt_t) TMath::Abs( xi->GetBindex() );
929 // Care track label can be negative in MC production (linked with the track quality)
930 // However = normally, not the case for track index ...
931
932 // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
933 if(lBachIdx == lIdxNegXi) {
934 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
935 }
936 if(lBachIdx == lIdxPosXi) {
937 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
938 }
939
940 AliESDtrack *pTrackXi = lESDevent->GetTrack( lIdxPosXi );
941 AliESDtrack *nTrackXi = lESDevent->GetTrack( lIdxNegXi );
942 AliESDtrack *bachTrackXi = lESDevent->GetTrack( lBachIdx );
943
944 if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
945 AliWarning("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
946 continue;
947 }
948
949 fTreeCascVarPosEta = pTrackXi->Eta();
950 fTreeCascVarNegEta = nTrackXi->Eta();
951 fTreeCascVarBachEta = bachTrackXi->Eta();
952
953 //------------------------------------------------
954 // TPC dEdx information
955 //------------------------------------------------
956 fTreeCascVarNegNSigmaPion = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kPion );
957 fTreeCascVarNegNSigmaProton = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kProton );
958 fTreeCascVarPosNSigmaPion = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kPion );
959 fTreeCascVarPosNSigmaProton = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kProton );
960 fTreeCascVarBachNSigmaPion = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kPion );
961 fTreeCascVarBachNSigmaKaon = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kKaon );
962
963 //------------------------------------------------
964 // TPC Number of clusters info
965 // --- modified to save the smallest number
966 // --- of TPC clusters for the 3 tracks
967 //------------------------------------------------
968
969 lPosTPCClusters = pTrackXi->GetTPCNcls();
970 lNegTPCClusters = nTrackXi->GetTPCNcls();
971 lBachTPCClusters = bachTrackXi->GetTPCNcls();
972
973 // 1 - Poor quality related to TPCrefit
974 ULong_t pStatus = pTrackXi->GetStatus();
975 ULong_t nStatus = nTrackXi->GetStatus();
976 ULong_t bachStatus = bachTrackXi->GetStatus();
977 if ((pStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
978 if ((nStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
979 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach. track has no TPCrefit ... continue!"); continue; }
980 // 2 - Poor quality related to TPC clusters: lowest cut of 70 clusters
981 if(lPosTPCClusters < 70) { AliWarning("Pb / V0 Pos. track has less than 70 TPC clusters ... continue!"); continue; }
982 if(lNegTPCClusters < 70) { AliWarning("Pb / V0 Neg. track has less than 70 TPC clusters ... continue!"); continue; }
983 if(lBachTPCClusters < 70) { AliWarning("Pb / Bach. track has less than 70 TPC clusters ... continue!"); continue; }
984 Int_t leastnumberofclusters = 1000;
985 if( lPosTPCClusters < leastnumberofclusters ) leastnumberofclusters = lPosTPCClusters;
986 if( lNegTPCClusters < leastnumberofclusters ) leastnumberofclusters = lNegTPCClusters;
987 if( lBachTPCClusters < leastnumberofclusters ) leastnumberofclusters = lBachTPCClusters;
988
989 lInvMassLambdaAsCascDghter = xi->GetEffMass();
990 // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
991 lDcaV0DaughtersXi = xi->GetDcaV0Daughters();
992 //lV0Chi2Xi = xi->GetChi2V0();
993
994 lV0CosineOfPointingAngleXi = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
995 lBestPrimaryVtxPos[1],
996 lBestPrimaryVtxPos[2] );
997
998 lDcaV0ToPrimVertexXi = xi->GetD( lBestPrimaryVtxPos[0],
999 lBestPrimaryVtxPos[1],
1000 lBestPrimaryVtxPos[2] );
1001
1002 lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
1003 lBestPrimaryVtxPos[1],
1004 lMagneticField ) );
1005 // Note : AliExternalTrackParam::GetD returns an algebraic value ...
1006
1007 xi->GetXYZ( lPosV0Xi[0], lPosV0Xi[1], lPosV0Xi[2] );
1008 lV0RadiusXi = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1] );
1009
1010 lDcaPosToPrimVertexXi = TMath::Abs( pTrackXi ->GetD( lBestPrimaryVtxPos[0],
1011 lBestPrimaryVtxPos[1],
1012 lMagneticField ) );
1013
1014 lDcaNegToPrimVertexXi = TMath::Abs( nTrackXi ->GetD( lBestPrimaryVtxPos[0],
1015 lBestPrimaryVtxPos[1],
1016 lMagneticField ) );
1017
1018 // - II.Step 4 : around effective masses (ESD)
1019 // ~ change mass hypotheses to cover all the possibilities : Xi-/+, Omega -/+
1020
1021 if( bachTrackXi->Charge() < 0 ) {
1022 lV0quality = 0.;
1023 xi->ChangeMassHypothesis(lV0quality , 3312);
1024 // Calculate the effective mass of the Xi- candidate.
1025 // pdg code 3312 = Xi-
1026 lInvMassXiMinus = xi->GetEffMassXi();
1027
1028 lV0quality = 0.;
1029 xi->ChangeMassHypothesis(lV0quality , 3334);
1030 // Calculate the effective mass of the Xi- candidate.
1031 // pdg code 3334 = Omega-
1032 lInvMassOmegaMinus = xi->GetEffMassXi();
1033
1034 lV0quality = 0.;
1035 xi->ChangeMassHypothesis(lV0quality , 3312); // Back to default hyp.
1036 }// end if negative bachelor
1037
1038
1039 if( bachTrackXi->Charge() > 0 ){
1040 lV0quality = 0.;
1041 xi->ChangeMassHypothesis(lV0quality , -3312);
1042 // Calculate the effective mass of the Xi+ candidate.
1043 // pdg code -3312 = Xi+
1044 lInvMassXiPlus = xi->GetEffMassXi();
1045
1046 lV0quality = 0.;
1047 xi->ChangeMassHypothesis(lV0quality , -3334);
1048 // Calculate the effective mass of the Xi+ candidate.
1049 // pdg code -3334 = Omega+
1050 lInvMassOmegaPlus = xi->GetEffMassXi();
1051
1052 lV0quality = 0.;
1053 xi->ChangeMassHypothesis(lV0quality , -3312); // Back to "default" hyp.
1054 }// end if positive bachelor
1055 // - II.Step 6 : extra info for QA (ESD)
1056 // miscellaneous pieces of info that may help regarding data quality assessment.
1057 //-------------
1058
1059 xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
1060 lXiTransvMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY );
1061 lXiTotMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY + lXiMomZ*lXiMomZ );
1062
1063 xi->GetBPxPyPz( lBachMomX, lBachMomY, lBachMomZ );
1064 //lBachTransvMom = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY );
1065 //lBachTotMom = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY + lBachMomZ*lBachMomZ );
1066
1067 lChargeXi = xi->Charge();
1068
1069 //lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
1070
1071 lRapXi = xi->RapXi();
1072 lRapOmega = xi->RapOmega();
1073 //lEta = xi->Eta();
1074 //lTheta = xi->Theta() *180.0/TMath::Pi();
1075 //lPhi = xi->Phi() *180.0/TMath::Pi();
1076 //lAlphaXi = xi->AlphaXi();
1077 //lPtArmXi = xi->PtArmXi();
1078
1079//------------------------------------------------
1080// Associate Cascade Candidates to Monte Carlo!
1081//------------------------------------------------
1082
1083//Warning: Not using Continues... Need to fill tree later!
1084
1085 Int_t lPDGCodeCascade = 0;
1086
1087 Int_t lPID_BachMother = 0;
1088 Int_t lPID_NegMother = 0;
1089 Int_t lPID_PosMother = 0;
1090
1091
1092 fTreeCascVarPIDPositive = 0;
1093 fTreeCascVarPIDNegative = 0;
1094 fTreeCascVarPIDBachelor = 0;
1095
1096
1097 if(fDebug > 5)
1098 cout << "MC EventNumber : " << lMCevent->Header()->GetEvent()
1099 << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
1100
1101 // - Step 4.1 : level of the V0 daughters
1102
1103//----------------------------------------
1104// Regular MC ASSOCIATION STARTS HERE
1105//----------------------------------------
1106
1107 Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
1108 // Abs value = needed ! question of quality track association ...
1109 Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
1110 Int_t lblBach = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
1111
1112 TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
1113 TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
1114 TParticle* mcBach = lMCstack->Particle( lblBach );
1115
1116 fTreeCascVarPosTransMomMC = mcPosV0Dghter->Pt();
1117 fTreeCascVarNegTransMomMC = mcNegV0Dghter->Pt();
1118
1119 fTreeCascVarPIDPositive = mcPosV0Dghter -> GetPdgCode();
1120 fTreeCascVarPIDNegative = mcNegV0Dghter -> GetPdgCode();
1121 fTreeCascVarPIDBachelor = mcBach->GetPdgCode();
1122
1123 // - Step 4.2 : level of the Xi daughters
1124
1125 Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ;
1126 Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
1127
1128 //Rather uncivilized: Open brackets for each 'continue'
1129 if(! (lblMotherPosV0Dghter != lblMotherNegV0Dghter) ) { // same mother
1130 if(! (lblMotherPosV0Dghter < 0) ) { // mother != primary (!= -1)
1131 if(! (lblMotherNegV0Dghter < 0) ) {
1132
1133 // mothers = Lambda candidate ... a priori
1134
1135 TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
1136 TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );
1137
1138 // - Step 4.3 : level of Xi candidate
1139
1140 Int_t lblGdMotherPosV0Dghter = mcMotherPosV0Dghter->GetFirstMother() ;
1141 Int_t lblGdMotherNegV0Dghter = mcMotherNegV0Dghter->GetFirstMother() ;
1142
1143 if(! (lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter) ) {
1144 if(! (lblGdMotherPosV0Dghter < 0) ) { // primary lambda ...
1145 if(! (lblGdMotherNegV0Dghter < 0) ) { // primary lambda ...
1146
1147 // Gd mothers = Xi candidate ... a priori
1148
1149 TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
1150 TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );
1151
1152 Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother() );
1153
1154 // if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
1155 if(!(lblMotherBach != lblGdMotherPosV0Dghter)) { //same mother for bach and V0 daughters
1156
1157 TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
1158
1159 lPID_BachMother = mcMotherBach->GetPdgCode();
1160 lPID_NegMother = mcGdMotherPosV0Dghter->GetPdgCode();
1161 lPID_PosMother = mcGdMotherNegV0Dghter->GetPdgCode();
1162
1163 if(lPID_BachMother==lPID_NegMother && lPID_BachMother==lPID_PosMother){
1164 lPDGCodeCascade = lPID_BachMother;
1165 lXiTransvMomMC = mcMotherBach->Pt();
1166 if ( (mcMotherBach->Energy() + mcMotherBach->Pz()) / (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) !=0 ){
1167 lRapMC = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) / (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
1168 }
1169 }
1170
1171 }}}}}}} //Ends all conditionals above...
1172
1173 //----------------------------------------
1174 // Regular MC ASSOCIATION ENDS HERE
1175 //----------------------------------------
1176
1177 //------------------------------------------------
1178 // Set Variables for adding to tree
1179 //------------------------------------------------
1180
1181/* 1*/ fTreeCascVarCharge = lChargeXi;
1182/* 2*/ if(lInvMassXiMinus!=0) fTreeCascVarMassAsXi = lInvMassXiMinus;
1183/* 2*/ if(lInvMassXiPlus!=0) fTreeCascVarMassAsXi = lInvMassXiPlus;
1184/* 3*/ if(lInvMassOmegaMinus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaMinus;
1185/* 3*/ if(lInvMassOmegaPlus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaPlus;
1186/* 4*/ fTreeCascVarPt = lXiTransvMom;
1187/* 4*/ fTreeCascVarPtMC = lXiTransvMomMC;
1188/* 5*/ fTreeCascVarRapXi = lRapXi ;
1189/* 5*/ fTreeCascVarRapMC = lRapMC ;
1190/* 6*/ fTreeCascVarRapOmega = lRapOmega ;
1191/* 7*/ fTreeCascVarDCACascDaughters = lDcaXiDaughters;
1192/* 8*/ fTreeCascVarDCABachToPrimVtx = lDcaBachToPrimVertexXi;
1193/* 9*/ fTreeCascVarDCAV0Daughters = lDcaV0DaughtersXi;
1194/*10*/ fTreeCascVarDCAV0ToPrimVtx = lDcaV0ToPrimVertexXi;
1195/*11*/ fTreeCascVarDCAPosToPrimVtx = lDcaPosToPrimVertexXi;
1196/*12*/ fTreeCascVarDCANegToPrimVtx = lDcaNegToPrimVertexXi;
1197/*13*/ fTreeCascVarCascCosPointingAngle = lXiCosineOfPointingAngle;
1198/*14*/ fTreeCascVarCascRadius = lXiRadius;
1199/*15*/ fTreeCascVarV0Mass = lInvMassLambdaAsCascDghter;
1200/*16*/ fTreeCascVarV0CosPointingAngle = lV0CosineOfPointingAngleXi;
1201/*17*/ fTreeCascVarV0Radius = lV0RadiusXi;
1202/*20*/ fTreeCascVarLeastNbrClusters = leastnumberofclusters;
1203/*21*/ fTreeCascVarMultiplicity = lMultiplicity; //multiplicity, whatever that may be
1204
1205/*23*/ fTreeCascVarDistOverTotMom = TMath::Sqrt(
1206 TMath::Power( lPosXi[0] - lBestPrimaryVtxPos[0] , 2) +
1207 TMath::Power( lPosXi[1] - lBestPrimaryVtxPos[1] , 2) +
1208 TMath::Power( lPosXi[2] - lBestPrimaryVtxPos[2] , 2)
1209 );
1210/*23*/ fTreeCascVarDistOverTotMom /= (lXiTotMom+1e-13);
1211/*24*/ //Not specified here, it has been set already (TRunNumber)
1212/*25*/ fTreeCascVarPID = lPDGCodeCascade;
1213
1214//------------------------------------------------
1215// Fill Tree!
1216//------------------------------------------------
1217
1218// The conditional is meant to decrease excessive
1219// memory usage! Be careful when loosening the
1220// cut!
1221
1222 //Xi Mass window: 150MeV wide
1223 //Omega mass window: 150MeV wide
1224
1225 if( (fTreeCascVarMassAsXi<1.32+0.075&&fTreeCascVarMassAsXi>1.32-0.075) ||
1226 (fTreeCascVarMassAsOmega<1.68+0.075&&fTreeCascVarMassAsOmega>1.68-0.075) ){
1227 fTreeCascade->Fill();
1228 }
1229
1230//------------------------------------------------
1231// Fill tree over.
1232//------------------------------------------------
1233
1234 }// end of the Cascade loop (ESD or AOD)
1235
1236 // Post output data.
1237 PostData(1, fListHist);
1238 PostData(2, fTreeCascade);
1239}
1240
1241//________________________________________________________________________
1242void AliAnalysisTaskExtractPerformanceCascade::Terminate(Option_t *)
1243{
1244 // Draw result to the screen
1245 // Called once at the end of the query
1246
1247 TList *cRetrievedList = 0x0;
1248 cRetrievedList = (TList*)GetOutputData(1);
1249 if(!cRetrievedList){
1250 Printf("ERROR - AliAnalysisTaskExtractCascade : ouput data container list not available\n");
1251 return;
1252 }
1253
1254 fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt") );
1255 if (!fHistV0MultiplicityForTrigEvt) {
1256 Printf("ERROR - AliAnalysisTaskExtractCascade : fHistV0MultiplicityForTrigEvt not available");
1257 return;
1258 }
1259
1260 TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractCascade","V0 Multiplicity",10,10,510,510);
1261 canCheck->cd(1)->SetLogy();
1262
1263 fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1264 fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1265}
1266
1267//----------------------------------------------------------------------------
1268
1269Double_t AliAnalysisTaskExtractPerformanceCascade::MyRapidity(Double_t rE, Double_t rPz) const
1270{
1271 // Local calculation for rapidity
1272 Double_t ReturnValue = -100;
1273 if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){
1274 ReturnValue = 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1275 }
1276 return ReturnValue;
1277}