]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/GammaConv/AliV0Reader.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Transition PWG4 --> PWGGA
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliV0Reader.cxx
CommitLineData
a0b94e5c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2eedd4ed 3 * *
4 * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *
5 * Version 1.0 *
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. *
a0b94e5c 14 **************************************************************************/
15
16////////////////////////////////////////////////
17//---------------------------------------------
18// Class used to do analysis on conversion pairs
19//---------------------------------------------
20////////////////////////////////////////////////
21
22// --- ROOT system ---
23#include <TMath.h>
24
25//---- ANALYSIS system ----
26#include "AliV0Reader.h"
27#include "AliAnalysisManager.h"
28#include "AliESDInputHandler.h"
67381a40 29#include "AliPID.h"
a0b94e5c 30#include "AliESDtrack.h"
31#include "AliMCEvent.h"
32#include "AliKFVertex.h"
3b77b2d1 33#include "AliKFParticle.h"
a0b94e5c 34#include "AliStack.h"
35#include "AliMCEventHandler.h"
10d100d4 36#include "AliESDpid.h"
5e55d806 37#include "AliGammaConversionBGHandler.h"
037dc2db 38#include "AliESDtrackCuts.h"
2ea0dd4a 39#include "TRandom3.h"
3b77b2d1 40#include "AliGenCocktailEventHeader.h"
41#include "TList.h"
a0b94e5c 42
43class iostream;
44class AliESDv0;
45class TFormula;
2ea0dd4a 46class TRandom3;
a0b94e5c 47
48using namespace std;
49
50ClassImp(AliV0Reader)
51
52
9c1cb6f7 53AliESDpid* AliV0Reader::fgESDpid = 0x0;
a0b94e5c 54
55AliV0Reader::AliV0Reader() :
2eedd4ed 56 TObject(),
57 fMCStack(NULL),
58 // fMCTruth(NULL),
59 fMCEvent(NULL), // for CF
60 fChain(NULL),
61 // fESDHandler(NULL),
62 fESDEvent(NULL),
63 fCFManager(NULL),
64 //fESDpid(NULL),
65 fHistograms(NULL),
66 fCurrentV0IndexNumber(0),
67 fCurrentV0(NULL),
68 fCurrentNegativeKFParticle(NULL),
69 fCurrentPositiveKFParticle(NULL),
70 fCurrentMotherKFCandidate(NULL),
71 fCurrentNegativeESDTrack(NULL),
72 fCurrentPositiveESDTrack(NULL),
73 fNegativeTrackLorentzVector(NULL),
74 fPositiveTrackLorentzVector(NULL),
75 fMotherCandidateLorentzVector(NULL),
76 fCurrentXValue(0),
77 fCurrentYValue(0),
78 fCurrentZValue(0),
79 fPositiveTrackPID(0),
80 fNegativeTrackPID(0),
81 fNegativeMCParticle(NULL),
82 fPositiveMCParticle(NULL),
83 fMotherMCParticle(NULL),
84 fMotherCandidateKFMass(0),
85 fMotherCandidateKFWidth(0),
86 fUseKFParticle(kTRUE),
87 fUseESDTrack(kFALSE),
88 fDoMC(kFALSE),
89 fMaxVertexZ(100.),// 100 cm(from the 0)
90 fMaxR(10000),// 100 meter(outside of ALICE)
91 fMinR(0),// 100 meter(outside of ALICE)
92 fEtaCut(0.),
93 fEtaCutMin(-0.1),
94 fRapidityMesonCut(0.),
95 fPtCut(0.),
96 fSinglePtCut(0.),
97 fMaxZ(0.),
98 fMinClsTPC(0.),
99 fMinClsTPCToF(0.),
100 fLineCutZRSlope(0.),
101 fLineCutZValue(0.),
102 fLineCutZRSlopeMin(0.),
103 fLineCutZValueMin(0.),
104 fChi2CutConversion(0.),
105 fChi2CutMeson(0.),
106 fAlphaCutMeson(1.),
107 fAlphaMinCutMeson(0.),
108 fPIDProbabilityCutNegativeParticle(0),
109 fPIDProbabilityCutPositiveParticle(0),
110 fDodEdxSigmaCut(kFALSE),
111 fDoTOFsigmaCut(kFALSE), // RRnewTOF
112 fPIDnSigmaAboveElectronLine(100),
113 fPIDnSigmaBelowElectronLine(-100),
114 fTofPIDnSigmaAboveElectronLine(100), // RRnewTOF
115 fTofPIDnSigmaBelowElectronLine(-100), // RRnewTOF
116 fPIDnSigmaAbovePionLine(-100),
117 fPIDnSigmaAbovePionLineHighPt(-100),
118 fPIDMinPnSigmaAbovePionLine(100),
119 fPIDMaxPnSigmaAbovePionLine(100),
120 fDoKaonRejectionLowP(kFALSE),
121 fDoProtonRejectionLowP(kFALSE),
122 fDoPionRejectionLowP(kFALSE),
123 fPIDnSigmaAtLowPAroundKaonLine(0),
124 fPIDnSigmaAtLowPAroundProtonLine(0),
125 fPIDnSigmaAtLowPAroundPionLine(0),
126 fPIDMinPKaonRejectionLowP(0),
127 fPIDMinPProtonRejectionLowP(0),
128 fPIDMinPPionRejectionLowP(0),
129 fDoQtGammaSelection(kFALSE),
130 fDoHighPtQtGammaSelection(kFALSE), // RRnew
131 fQtMax(100.),
132 fHighPtQtMax(100.), // RRnew
133 fPtBorderForQt(100.), // RRnew
134 fXVertexCut(0.),
135 fYVertexCut(0.),
136 fZVertexCut(0.),
3b77b2d1 137 fPsiPairCut(0.),
138 fCosinePointCut(0.),
2eedd4ed 139 fNSigmaMass(0.),
140 fUseImprovedVertex(kFALSE),
141 fUseOwnXYZCalculation(kFALSE),
142 fUseConstructGamma(kFALSE),
143 fDoCF(kFALSE),
144 fUseEtaMinCut(kFALSE),
145 fUseOnFlyV0Finder(kTRUE),
146 fUpdateV0AlreadyCalled(kFALSE),
147 fCurrentEventGoodV0s(NULL),
148 fV0Pindex(),
149 fV0Nindex(),
150// fPreviousEventGoodV0s(),
151 fCalculateBackground(kFALSE),
152 fBGEventHandler(NULL),
153 fBGEventInitialized(kFALSE),
154 fEsdTrackCuts(NULL),
155 fNumberOfESDTracks(0),
156 fNEventsForBGCalculation(20),
157 fUseChargedTrackMultiplicityForBG(kTRUE),
158 fNumberOfGoodV0s(0),
159 fIsHeavyIon(0),
160 fUseCorrectedTPCClsInfo(kFALSE),
161 fUseMCPSmearing(kTRUE),
162 fPBremSmearing(1.),
163 fPSigSmearing(0.),
164 fPSigSmearingCte(0.),
165 fRandom(0),
166 fBrem(NULL),
167 fDoPhotonAsymmetryCut(0),
3b77b2d1 168 fdoESDQtCut(0),
2eedd4ed 169 fMinPPhotonAsymmetryCut(100.),
3b77b2d1 170 fMinPhotonAsymmetry(0.),
171 fExcludeBackgroundEventForGammaCorrection(0.),
172 fNumberOfPrimerisFromHijingAndPythia(0)
a0b94e5c 173{
2eedd4ed 174 //fESDpid = new AliESDpid;
a0b94e5c 175}
176
177
178AliV0Reader::AliV0Reader(const AliV0Reader & original) :
2eedd4ed 179 TObject(original),
180 fMCStack(original.fMCStack),
181 // fMCTruth(original.fMCTruth),
182 fMCEvent(original.fMCEvent), // for CF
183 fChain(original.fChain),
184 // fESDHandler(original.fESDHandler),
185 fESDEvent(original.fESDEvent),
186 fCFManager(original.fCFManager),
187 // fESDpid(original.fESDpid),
188 fHistograms(original.fHistograms),
189 fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),
190 fCurrentV0(original.fCurrentV0),
191 fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),
192 fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),
193 fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),
194 fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),
195 fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),
196 fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),
197 fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),
198 fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),
199 fCurrentXValue(original.fCurrentXValue),
200 fCurrentYValue(original.fCurrentYValue),
201 fCurrentZValue(original.fCurrentZValue),
202 fPositiveTrackPID(original.fPositiveTrackPID),
203 fNegativeTrackPID(original.fNegativeTrackPID),
204 fNegativeMCParticle(original.fNegativeMCParticle),
205 fPositiveMCParticle(original.fPositiveMCParticle),
206 fMotherMCParticle(original.fMotherMCParticle),
207 fMotherCandidateKFMass(original.fMotherCandidateKFMass),
208 fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),
209 fUseKFParticle(kTRUE),
210 fUseESDTrack(kFALSE),
211 fDoMC(kFALSE),
212 fMaxVertexZ(original.fMaxVertexZ),
213 fMaxR(original.fMaxR),
214 fMinR(original.fMinR),
215 fEtaCut(original.fEtaCut),
216 fEtaCutMin(original.fEtaCutMin),
217 fRapidityMesonCut(original.fRapidityMesonCut),
218 fPtCut(original.fPtCut),
219 fSinglePtCut(original.fSinglePtCut),
220 fMaxZ(original.fMaxZ),
221 fMinClsTPC(original.fMinClsTPC),
222 fMinClsTPCToF(original.fMinClsTPCToF),
223 fLineCutZRSlope(original.fLineCutZRSlope),
224 fLineCutZValue(original.fLineCutZValue),
225 fLineCutZRSlopeMin(original.fLineCutZRSlopeMin),
226 fLineCutZValueMin(original.fLineCutZValueMin),
227 fChi2CutConversion(original.fChi2CutConversion),
228 fChi2CutMeson(original.fChi2CutMeson),
229 fAlphaCutMeson(original.fAlphaCutMeson),
230 fAlphaMinCutMeson(original.fAlphaMinCutMeson),
231 fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),
232 fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),
233 fDodEdxSigmaCut(original.fDodEdxSigmaCut),
234 fDoTOFsigmaCut(original.fDoTOFsigmaCut), // RRnewTOF
235 fPIDnSigmaAboveElectronLine(original.fPIDnSigmaAboveElectronLine),
236 fPIDnSigmaBelowElectronLine(original.fPIDnSigmaBelowElectronLine),
237 fTofPIDnSigmaAboveElectronLine(original.fTofPIDnSigmaAboveElectronLine), // RRnewTOF
238 fTofPIDnSigmaBelowElectronLine(original.fTofPIDnSigmaBelowElectronLine), // RRnewTOF
239 fPIDnSigmaAbovePionLine(original.fPIDnSigmaAbovePionLine),
240 fPIDnSigmaAbovePionLineHighPt(original.fPIDnSigmaAbovePionLineHighPt),
241 fPIDMinPnSigmaAbovePionLine(original.fPIDMinPnSigmaAbovePionLine),
242 fPIDMaxPnSigmaAbovePionLine(original.fPIDMaxPnSigmaAbovePionLine),
243 fDoKaonRejectionLowP(original.fDoKaonRejectionLowP),
244 fDoProtonRejectionLowP(original.fDoProtonRejectionLowP),
245 fDoPionRejectionLowP(original.fDoPionRejectionLowP),
246 fPIDnSigmaAtLowPAroundKaonLine(original.fPIDnSigmaAtLowPAroundKaonLine),
247 fPIDnSigmaAtLowPAroundProtonLine(original.fPIDnSigmaAtLowPAroundProtonLine),
248 fPIDnSigmaAtLowPAroundPionLine(original.fPIDnSigmaAtLowPAroundPionLine),
249 fPIDMinPKaonRejectionLowP(original.fPIDMinPKaonRejectionLowP),
250 fPIDMinPProtonRejectionLowP(original.fPIDMinPProtonRejectionLowP),
251 fPIDMinPPionRejectionLowP(original.fPIDMinPPionRejectionLowP),
252 fDoQtGammaSelection(original.fDoQtGammaSelection),
253 fDoHighPtQtGammaSelection(original.fDoHighPtQtGammaSelection), // RRnew
254 fQtMax(original.fQtMax),
255 fHighPtQtMax(original.fHighPtQtMax), // RRnew
256 fPtBorderForQt(original.fPtBorderForQt), // RRnew
257 fXVertexCut(original.fXVertexCut),
258 fYVertexCut(original.fYVertexCut),
259 fZVertexCut(original.fZVertexCut),
3b77b2d1 260 fPsiPairCut(original.fPsiPairCut),
261 fCosinePointCut(original.fCosinePointCut),
2eedd4ed 262 fNSigmaMass(original.fNSigmaMass),
263 fUseImprovedVertex(original.fUseImprovedVertex),
264 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
265 fUseConstructGamma(original.fUseConstructGamma),
266 fDoCF(original.fDoCF),
267 fUseEtaMinCut(original.fUseEtaMinCut),
268 fUseOnFlyV0Finder(original.fUseOnFlyV0Finder),
269 fUpdateV0AlreadyCalled(original.fUpdateV0AlreadyCalled),
270 fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
271 fV0Pindex(original.fV0Pindex),
272 fV0Nindex(original.fV0Nindex),
273 // fPreviousEventGoodV0s(original.fPreviousEventGoodV0s),
274 fCalculateBackground(original.fCalculateBackground),
275 fBGEventHandler(original.fBGEventHandler),
276 fBGEventInitialized(original.fBGEventInitialized),
277 fEsdTrackCuts(original.fEsdTrackCuts),
278 fNumberOfESDTracks(original.fNumberOfESDTracks),
279 fNEventsForBGCalculation(original.fNEventsForBGCalculation),
280 fUseChargedTrackMultiplicityForBG(original.fUseChargedTrackMultiplicityForBG),
281 fNumberOfGoodV0s(original.fNumberOfGoodV0s),
282 fIsHeavyIon(original.fIsHeavyIon),
283 fUseCorrectedTPCClsInfo(original.fUseCorrectedTPCClsInfo),
284 fUseMCPSmearing(original.fUseMCPSmearing),
285 fPBremSmearing(original.fPBremSmearing),
286 fPSigSmearing(original.fPSigSmearing),
287 fPSigSmearingCte(original.fPSigSmearingCte),
288 fRandom(original.fRandom),
289 fBrem(original.fBrem),
290 fDoPhotonAsymmetryCut(original.fDoPhotonAsymmetryCut),
3b77b2d1 291 fdoESDQtCut(original.fdoESDQtCut),
2eedd4ed 292 fMinPPhotonAsymmetryCut(original.fMinPPhotonAsymmetryCut),
3b77b2d1 293 fMinPhotonAsymmetry(original.fMinPhotonAsymmetry),
294 fExcludeBackgroundEventForGammaCorrection(original.fExcludeBackgroundEventForGammaCorrection),
295 fNumberOfPrimerisFromHijingAndPythia(original.fNumberOfPrimerisFromHijingAndPythia)
a0b94e5c 296{
297
298}
299
300
301AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)
302{
2eedd4ed 303 // assignment operator
304 return *this;
a0b94e5c 305}
9640a3d1 306AliV0Reader::~AliV0Reader()
307{
2eedd4ed 308 // if(fESDpid){
309 // delete fESDpid;
310 //}
9640a3d1 311}
a0b94e5c 312
48682642 313//____________________________________________________________________________
314void AliV0Reader::SetInputAndMCEvent(AliVEvent* esd, AliMCEvent* mc) {
2eedd4ed 315 // Connect the data pointers
48682642 316
2eedd4ed 317 SetInputEvent(esd);
318 SetMC(mc);
48682642 319
320}
321
322
a0b94e5c 323void AliV0Reader::Initialize(){
2eedd4ed 324 //see header file for documentation
325
326 fUpdateV0AlreadyCalled = kFALSE;
327
328 /*
329 // Get the input handler from the manager
330 fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
331 if(fESDHandler == NULL){
332 //print warning here
333 }
334
335 // Get pointer to esd event from input handler
336 fESDEvent = fESDHandler->GetEvent();
337 if(fESDEvent == NULL){
338 //print warning here
339 }
340
341 //Get pointer to MCTruth
342 fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
343 */
344
345
346
347 // fMCTruth = mcH->MCEvent();
348 // fMC = mcH->MCEvent();
349 // stack = fMC->Stack();
350
351
352 //if(fMCTruth == NULL){
353 //print warning here
354 // fDoMC = kFALSE;
355 //}
356
357 if(fMCEvent == NULL){
358 fDoMC = kFALSE;
359 }
360
361 //Get pointer to the mc stack
362 // if(fMCTruth){
363 if(fMCEvent){
364 fMCStack = fMCEvent->Stack();
365 if(fMCStack == NULL){
366 //print warning here
367 }
368 // Better parameters for MonteCarlo from A. Kalweit 2010/01/8
369// fESDpid->GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
370// 1.75295e+01,
371// 3.40030e-09,
372// 1.96178e+00,
373// 3.91720e+00);
374 }
375 else{
376 // Better parameters for data from A. Kalweit 2010/01/8
377 // fESDpid->GetTPCResponse().SetBetheBlochParameters(0.0283086,
378// 2.63394e+01,
379// 5.04114e-11,
380// 2.12543e+00,
381// 4.88663e+00);
382 }
a0b94e5c 383
2eedd4ed 384 // for CF
385 //Get pointer to the mc event
386 if(fDoCF && fDoMC){
387 //fMCEvent = fMCTruth->MCEvent();
388 if(fMCEvent == NULL){
389 //print warning here
390 fDoCF = kFALSE;
391 }
392 }
a0b94e5c 393
2eedd4ed 394 fUseEtaMinCut = kFALSE;
395 if ( fEtaCutMin != -0.1) {
396 fUseEtaMinCut = kTRUE;
397 }
398
399
400 AliKFParticle::SetField(fESDEvent->GetMagneticField());
401
402 // fCurrentEventGoodV0s = new TClonesArray("TClonesArray", 0);
403 if(fCurrentEventGoodV0s == NULL){
404 fCurrentEventGoodV0s = new TClonesArray("AliKFParticle", 0);
405 }
406
407 fV0Pindex.clear();
408 fV0Nindex.clear();
409
410 if(gRandom != NULL){
411 delete gRandom;
412 gRandom= new TRandom3(0);
413 }
414
415
416 if (fBrem == NULL){
417 fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
418 // tests done with 1.0e-14
419 fBrem->SetParameter(0,fPBremSmearing);
420 fBrem->SetNpx(100000);
421 }
422
423 if(fCalculateBackground == kTRUE){
424 if(fBGEventInitialized == kFALSE){
425
3b77b2d1 426 Double_t zBinLimitsArray[9];
2eedd4ed 427 zBinLimitsArray[0] = -50.00;
428 zBinLimitsArray[1] = -3.375;
429 zBinLimitsArray[2] = -1.605;
430 zBinLimitsArray[3] = -0.225;
431 zBinLimitsArray[4] = 1.065;
432 zBinLimitsArray[5] = 2.445;
433 zBinLimitsArray[6] = 4.245;
434 zBinLimitsArray[7] = 50.00;
435 zBinLimitsArray[8] = 1000.00;
436
3b77b2d1 437 Double_t multiplicityBinLimitsArray[6];
2eedd4ed 438 if(fUseChargedTrackMultiplicityForBG == kTRUE){
439 multiplicityBinLimitsArray[0] = 0;
440 multiplicityBinLimitsArray[1] = 8.5;
441 multiplicityBinLimitsArray[2] = 16.5;
442 multiplicityBinLimitsArray[3] = 27.5;
443 multiplicityBinLimitsArray[4] = 41.5;
444 multiplicityBinLimitsArray[5] = 100.;
445 if(fIsHeavyIon){
446 multiplicityBinLimitsArray[0] = 0;
447 multiplicityBinLimitsArray[1] = 200.;
448 multiplicityBinLimitsArray[2] = 500.;
449 multiplicityBinLimitsArray[3] = 1000.;
450 multiplicityBinLimitsArray[4] = 1500.;
451 multiplicityBinLimitsArray[5] = 3000.;
452 }
453 fBGEventHandler = new AliGammaConversionBGHandler(9,6,fNEventsForBGCalculation);
454 } else {
455 multiplicityBinLimitsArray[0] = 2;
456 multiplicityBinLimitsArray[1] = 3;
457 multiplicityBinLimitsArray[2] = 4;
458 multiplicityBinLimitsArray[3] = 5;
459 multiplicityBinLimitsArray[4] = 9999;
460 if(fIsHeavyIon){
461 multiplicityBinLimitsArray[0] = 2;
462 multiplicityBinLimitsArray[1] = 10;
463 multiplicityBinLimitsArray[2] = 30;
464 multiplicityBinLimitsArray[3] = 50;
465 multiplicityBinLimitsArray[4] = 9999;
466 }
467
468 fBGEventHandler = new AliGammaConversionBGHandler(9,5,fNEventsForBGCalculation);
469 }
470
471
472
473 /*
474 // ---------------------------------
475 Double_t *zBinLimitsArray = new Double_t[1];
476 zBinLimitsArray[0] = 999999.00;
477
478 Double_t *multiplicityBinLimitsArray= new Double_t[1];
479 multiplicityBinLimitsArray[0] = 99999999.00;
480 fBGEventHandler = new AliGammaConversionBGHandler(1,1,10);
481 // ---------------------------------
482 */
483 fBGEventHandler->Initialize(zBinLimitsArray, multiplicityBinLimitsArray);
484 fBGEventInitialized = kTRUE;
485 }
486 }
3b77b2d1 487
488 if(fDoMC && fExcludeBackgroundEventForGammaCorrection){
489 fNumberOfPrimerisFromHijingAndPythia = GetNumberOfHijingPlusPythiaPrimeries();
490 }
a0b94e5c 491}
492
493AliESDv0* AliV0Reader::GetV0(Int_t index){
2eedd4ed 494 //see header file for documentation
495 fCurrentV0 = fESDEvent->GetV0(index);
496 UpdateV0Information();
497 return fCurrentV0;
a0b94e5c 498}
499
c8206114 500Int_t AliV0Reader::GetNumberOfContributorsVtx(){
2eedd4ed 501 // returns number of contributors to the vertex
502 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
503 return fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
504 }
505
506 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
507 // return 0;
508 //-AM test pi0s without SPD only vertex
509 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
510 return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
511
512 }
513 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
514 cout<<"number of contributors from bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
515 return 0;
516 }
517 }
518 return 0;
c8206114 519}
2eedd4ed 520
a0b94e5c 521Bool_t AliV0Reader::CheckForPrimaryVertex(){
2eedd4ed 522 //see headerfile for documentation
523 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
524 return 1;
525 }
526 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
527 // SPD vertex
528 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
529 // return 0;
530 //-AM test pi0s without SPD only vertex
531 //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
532 return 1;
533 }
534 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
535 // cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
536 return 0;
537 }
538 }
539 return 0;
540 // return fESDEvent->GetPrimaryVertex()->GetNContributors()>0;
a0b94e5c 541}
542
10e3319b 543Bool_t AliV0Reader::CheckForPrimaryVertexZ(){
2eedd4ed 544 //see headerfile for documentation
545 if(TMath::Abs(fESDEvent->GetPrimaryVertex()->GetZ())<GetMaxVertexZ()){
546 return kTRUE;
547 }else{
548 return kFALSE;
549 }
550 return kTRUE;
10e3319b 551}
77880bd8 552
3b77b2d1 553
554Bool_t AliV0Reader::CheckV0FinderStatus(AliESDv0 *v0){
2eedd4ed 555 // see headerfile for documentation
556 if(fUseOnFlyV0Finder){
3b77b2d1 557 if(!v0->GetOnFlyStatus()){
2eedd4ed 558 return kFALSE;
559 }
560 }
561 if(!fUseOnFlyV0Finder){
3b77b2d1 562 if(v0->GetOnFlyStatus()){
2eedd4ed 563 return kFALSE;
564 }
565 }
566 return kTRUE;
77880bd8 567}
568
569
570
3b77b2d1 571
a0b94e5c 572Bool_t AliV0Reader::NextV0(){
2eedd4ed 573 //see header file for documentation
574 Bool_t iResult=kFALSE;
575 while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){
576 fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);
037dc2db 577
2eedd4ed 578 fUpdateV0AlreadyCalled=kFALSE;
3b77b2d1 579
2eedd4ed 580 if(fHistograms != NULL){
581 fHistograms->FillHistogram("ESD_AllV0s_InvMass",GetMotherCandidateMass());
3b77b2d1 582 fHistograms->FillHistogram("ESD_AllV0s_Pt",GetMotherCandidatePt());
2eedd4ed 583 }
a0b94e5c 584
2eedd4ed 585 // moved it up here so that the correction framework can access pt and eta information
586 if(UpdateV0Information() == kFALSE){
587 fCurrentV0IndexNumber++;
588 continue;
589 }
cb90a330 590
3b77b2d1 591
592 if(fDoMC && fExcludeBackgroundEventForGammaCorrection){ // Remove all V0s from BGEvent
593 Bool_t isFromBGEvent = kFALSE;
594 isFromBGEvent = IsParticleFromBGEvent(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
595 if(isFromBGEvent){
596 fHistograms->FillHistogram("ESD_CutMCBgEvent_InvMass",GetMotherCandidateMass());
597 fHistograms->FillHistogram("ESD_CutMCBgEvent_Pt",GetMotherCandidatePt());
598 fCurrentV0IndexNumber++;
599 continue;
600 }
601 isFromBGEvent = IsParticleFromBGEvent(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
602 if(isFromBGEvent){
603 fHistograms->FillHistogram("ESD_CutMCBgEvent_InvMass",GetMotherCandidateMass());
604 fHistograms->FillHistogram("ESD_CutMCBgEvent_Pt",GetMotherCandidatePt());
605 fCurrentV0IndexNumber++;
606 continue;
607 }
608 }
609
610
611
2eedd4ed 612 Double_t containerInput[3];
613 if(fDoCF){
614 containerInput[0] = GetMotherCandidatePt();
615 containerInput[1] = GetMotherCandidateEta();
616 containerInput[2] = GetMotherCandidateMass();
617 }
618 /*
619 if(fDoCF){
620 containerInput[0] = GetMotherCandidatePt();
621 containerInput[1] = GetMotherCandidateEta();
622 containerInput[2] = GetMotherCandidateMass();
623
624 fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign); // for CF
625 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit); // for CF
626 fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks); // for CF
627 }
628 */
629
630 //checks if on the fly mode is set
3b77b2d1 631 if ( !CheckV0FinderStatus(fCurrentV0) ){
632
2eedd4ed 633 if(fHistograms != NULL){
634 fHistograms->FillHistogram("ESD_CutGetOnFly_InvMass",GetMotherCandidateMass());
3b77b2d1 635 fHistograms->FillHistogram("ESD_CutGetOnFly_Pt",GetMotherCandidatePt());
2eedd4ed 636 }
637 fCurrentV0IndexNumber++;
638 continue;
639 }
640 if(fDoCF){
641 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly); // for CF
642 }
643
644 if(fHistograms != NULL){
645 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_InvMass",GetMotherCandidateMass());
3b77b2d1 646 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_Pt",GetMotherCandidatePt());
2eedd4ed 647 }
48682642 648
3b77b2d1 649 Double_t armenterosQtAlpha[2] = {0,0};
650 Double_t armenterosQtAlphaKF[2] = {0,0};
651 Double_t armenterosQtAlphaESD[2] = {0,0};
652 Double_t armenterosQtAlphaKFNew[2] = {0,0};
653 Double_t armenterosQtAlphaESDMC[2] = {0,0};
654 Double_t armenterosQtAlphaMC[2] = {0,0};
655
656 GetArmenterosQtAlpha(GetNegativeKFParticle(), // old KF way calculating Qt Alpha
657 GetPositiveKFParticle(),
658 GetMotherCandidateKFCombination(),
659 armenterosQtAlphaKF);
660 GetArmenterosQtAlpha(fCurrentV0,armenterosQtAlphaESD); // ESD way calculating Qt Alpha
661 GetArmenterosQtAlpha(GetNegativeKFParticle(), // new KF way calculating Qt Alpha
662 GetPositiveKFParticle(),
663 armenterosQtAlphaKFNew,fdoESDQtCut);
664
665 fHistograms->FillHistogram("ESD_AllV0sCurrentFinderKF_alfa_qt",armenterosQtAlphaKF[1],armenterosQtAlphaKF[0]);
666 fHistograms->FillHistogram("ESD_AllV0sCurrentFinderESD_alfa_qt",armenterosQtAlphaESD[1],armenterosQtAlphaESD[0]);
667 fHistograms->FillHistogram("ESD_AllV0sCurrentFinderKFNew_alfa_qt",armenterosQtAlphaKFNew[1],armenterosQtAlphaKFNew[0]);
668
669 if(fdoESDQtCut == 0){
670 armenterosQtAlpha[0] = armenterosQtAlphaKF[0];
671 armenterosQtAlpha[1] = armenterosQtAlphaKF[1];
672 }
673 else if(fdoESDQtCut == 1){
674 armenterosQtAlpha[0] = armenterosQtAlphaESD[0];
675 armenterosQtAlpha[1] = armenterosQtAlphaESD[1];
676 }
677
678 else if(fdoESDQtCut == 2 || fdoESDQtCut == 3){
679 armenterosQtAlpha[0] = armenterosQtAlphaKFNew[0];
680 armenterosQtAlpha[1] = armenterosQtAlphaKFNew[1];
681 }
682
2eedd4ed 683 if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){ // avoid like sign
684 // iResult=kFALSE;
685 if(fHistograms != NULL ){
686 fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
3b77b2d1 687 fHistograms->FillHistogram("ESD_CutLikeSign_Pt",GetMotherCandidatePt());
2eedd4ed 688 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
689 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
690 // fUpdateV0AlreadyCalled = kTRUE;
691 }
692 fCurrentV0IndexNumber++;
693 continue;
694 }
695 if(fDoCF){
696 fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign); // for CF
697 }
ebcfaa7e 698
2eedd4ed 699 if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ||
700 !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
701 // if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) ||
702 // !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
703 // iResult=kFALSE;
704 if(fHistograms != NULL){
705 fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
3b77b2d1 706 fHistograms->FillHistogram("ESD_CutRefit_Pt",GetMotherCandidatePt());
2eedd4ed 707 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
708 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
709 //fUpdateV0AlreadyCalled = kTRUE;
710 }
711 fCurrentV0IndexNumber++;
712 continue;
713 }
714 if(fDoCF){
715 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit); // for CF
716 }
ebcfaa7e 717
718
719
2eedd4ed 720 if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
721 fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
722 //iResult=kFALSE;
723 if(fHistograms != NULL ){
724 fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
3b77b2d1 725 fHistograms->FillHistogram("ESD_CutKink_Pt",GetMotherCandidatePt());
2eedd4ed 726 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
727 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
728 //fUpdateV0AlreadyCalled = kTRUE;
729 }
730 fCurrentV0IndexNumber++;
731 continue;
732 }
733
734 if(fDoCF){
735 fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks); // for CF
736 }
ebcfaa7e 737
2eedd4ed 738 if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
739 if(fHistograms != NULL){
740 fHistograms->FillHistogram("ESD_CutR_InvMass",GetMotherCandidateMass());
3b77b2d1 741 fHistograms->FillHistogram("ESD_CutR_Pt",GetMotherCandidatePt());
2eedd4ed 742 }
743 fCurrentV0IndexNumber++;
744 continue;
745 }
746 if(fDoCF){
747 fCFManager->GetParticleContainer()->Fill(containerInput,kStepR); // for CF
748 }
749 if(GetXYRadius()<fMinR){ // cuts on distance from collision point
750 if(fHistograms != NULL){
751 fHistograms->FillHistogram("ESD_CutMinR_InvMass",GetMotherCandidateMass());
3b77b2d1 752 fHistograms->FillHistogram("ESD_CutMinR_Pt",GetMotherCandidatePt());
2eedd4ed 753 }
754 fCurrentV0IndexNumber++;
755 continue;
756 }
757
758 //if((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue > GetXYRadius() ) { // cuts out regions where we do not reconstruct
759 if( GetXYRadius() <= ((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue)){
760 if(fHistograms != NULL){
761 fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
3b77b2d1 762 fHistograms->FillHistogram("ESD_CutLine_Pt",GetMotherCandidatePt());
2eedd4ed 763 }
764 fCurrentV0IndexNumber++;
765 continue;
766 } else if (fUseEtaMinCut && GetXYRadius() >= ((TMath::Abs(fCurrentZValue)*fLineCutZRSlopeMin)-fLineCutZValueMin )){
767 if(fHistograms != NULL){
768 fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
3b77b2d1 769 fHistograms->FillHistogram("ESD_CutLine_Pt",GetMotherCandidatePt());
2eedd4ed 770 }
771 fCurrentV0IndexNumber++;
772 continue;
773 }
cc5a88a2 774
2eedd4ed 775 if(fDoCF){
776 fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine); // for CF
777 }
cc5a88a2 778
2eedd4ed 779 if(TMath::Abs(fCurrentZValue) > fMaxZ ){ // cuts out regions where we do not reconstruct
780 if(fHistograms != NULL){
781 fHistograms->FillHistogram("ESD_CutZ_InvMass",GetMotherCandidateMass());
3b77b2d1 782 fHistograms->FillHistogram("ESD_CutZ_Pt",GetMotherCandidatePt());
2eedd4ed 783 }
784 fCurrentV0IndexNumber++;
785 continue;
786 }
787 if(fDoCF){
788 fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ); // for CF
789 }
cc5a88a2 790
2eedd4ed 791 if(fUseKFParticle){
792 if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut || TMath::Abs(fMotherCandidateLorentzVector->Eta())< fEtaCutMin){
793 if(fHistograms != NULL){
794 fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
3b77b2d1 795 fHistograms->FillHistogram("ESD_CutEta_Pt",GetMotherCandidatePt());
2eedd4ed 796 }
797 fCurrentV0IndexNumber++;
798 continue;
799 }
cc5a88a2 800
2eedd4ed 801 if(TMath::Abs(fCurrentNegativeKFParticle->GetEta())> fEtaCut || TMath::Abs(fCurrentNegativeKFParticle->GetEta())< fEtaCutMin){
802 if(fHistograms != NULL){
803 fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
3b77b2d1 804 fHistograms->FillHistogram("ESD_CutEta_Pt",GetMotherCandidatePt());
2eedd4ed 805 }
806 fCurrentV0IndexNumber++;
807 continue;
808 }
cc5a88a2 809
2eedd4ed 810 if(TMath::Abs(fCurrentPositiveKFParticle->GetEta())> fEtaCut || TMath::Abs(fCurrentPositiveKFParticle->GetEta())< fEtaCutMin){
811 if(fHistograms != NULL){
812 fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
3b77b2d1 813 fHistograms->FillHistogram("ESD_CutEta_Pt",GetMotherCandidatePt());
2eedd4ed 814 }
815 fCurrentV0IndexNumber++;
816 continue;
817 }
818 }
e0155998 819
3b77b2d1 820 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_Pt_Qt",GetMotherCandidatePt(),armenterosQtAlpha[0]);
e0155998 821
2eedd4ed 822 if(fDoMC){
823 if ( HasSameMCMother() == kTRUE){
3b77b2d1 824 GetArmenterosQtAlpha(fNegativeMCParticle,
2eedd4ed 825 fPositiveMCParticle,
826 fMotherMCParticle,
827 armenterosQtAlphaMC);
828 }
3b77b2d1 829 GetArmenterosQtAlpha(fNegativeMCParticle,
2eedd4ed 830 fPositiveMCParticle,
831 GetMotherCandidateKFCombination(),
832 armenterosQtAlphaESDMC );
833 }
3b77b2d1 834
835 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_goodtracks_alfa_qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 836 if( fCurrentNegativeKFParticle->GetPt()> 0.150 && fCurrentPositiveKFParticle->GetPt()> 0.150){
3b77b2d1 837 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_minPt_GT_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 838 }
839 if(fDoMC){
840 fHistograms->FillHistogram("ESD_TrueConvAllV0s_ESDMother_Alpha_Qt",armenterosQtAlphaESDMC[1],armenterosQtAlphaESDMC[0]);
841 if ( HasSameMCMother() == kTRUE){
842 fHistograms->FillHistogram("ESD_TrueConvSameMother_ESDMother_Alpha_Qt",armenterosQtAlphaESDMC[1],armenterosQtAlphaESDMC[0]);
843 fHistograms->FillHistogram("ESD_TrueConvSameMother_MCMother_Alpha_Qt",armenterosQtAlphaMC[1],armenterosQtAlphaMC[0]);
844 if (fMotherMCParticle->GetPdgCode() == 22 ){
3b77b2d1 845 fHistograms->FillHistogram("ESD_TrueConvGamma_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
846 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Qt",GetMotherCandidatePt(),armenterosQtAlpha[0]);
2eedd4ed 847 } else if ( fMotherMCParticle->GetPdgCode() == 310 ){
3b77b2d1 848 fHistograms->FillHistogram("ESD_TrueConvK0s_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 849 } else if ( fMotherMCParticle->GetPdgCode() == 113 ){
3b77b2d1 850 fHistograms->FillHistogram("ESD_TrueConvRho0_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 851 } else if ( fMotherMCParticle->GetPdgCode() == 333 ){
3b77b2d1 852 fHistograms->FillHistogram("ESD_TrueConvPhi_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 853 } else if ( (fMotherMCParticle->GetPdgCode() == 3122 || fMotherMCParticle->GetPdgCode() == -3122) ){
3b77b2d1 854 fHistograms->FillHistogram("ESD_TrueConvLambda_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 855 } else if ( (fMotherMCParticle->GetPdgCode() == 2114 || fMotherMCParticle->GetPdgCode() == -2114) ){
3b77b2d1 856 fHistograms->FillHistogram("ESD_TrueConvDelta_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 857 } else if ( (fMotherMCParticle->GetPdgCode() == 313 ||
858 fMotherMCParticle->GetPdgCode() == 323 ||
859 fMotherMCParticle->GetPdgCode() == -323 ) ){
3b77b2d1 860 fHistograms->FillHistogram("ESD_TrueConvKStar_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 861 } else {
3b77b2d1 862 fHistograms->FillHistogram("ESD_TrueConvUnknown_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 863 fHistograms->FillHistogram("ESD_TrueConvUnknown_Qt_PDG",fMotherMCParticle->GetPdgCode());
864 // cout << "unidentfied mother: pdg-C mother " << fMotherMCParticle->GetPdgCode() << " daughters " << fNegativeMCParticle->GetPdgCode() << "\t" << fPositiveMCParticle->GetPdgCode() << endl;
865 }
866 } else {
3b77b2d1 867 fHistograms->FillHistogram("ESD_TrueConvComb_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
2eedd4ed 868 }
869 }
870
871 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_E_dEdxP",fCurrentNegativeESDTrack->P(),fCurrentNegativeESDTrack->GetTPCsignal());
872 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_P_dEdxP",fCurrentPositiveESDTrack->P(),fCurrentPositiveESDTrack->GetTPCsignal());
873 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_E_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron));
874 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_P_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron));
875 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PiPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion));
876 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PiMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion));
877 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_KPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon));
878 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_KMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon));
879 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kProton));
880 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton));
881 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_MuPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kMuon));
882 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_MuMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kMuon));
883
884 if(fDodEdxSigmaCut == kTRUE){
885 if( fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
886 fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ||
887 fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
888 fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
889 //iResult=kFALSE;
890 if(fHistograms != NULL ){
891 fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
3b77b2d1 892 fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_Pt",GetMotherCandidatePt());
2eedd4ed 893 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
894 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
895 //fUpdateV0AlreadyCalled = kTRUE;
14d7f2f9 896 }
897 fCurrentV0IndexNumber++;
898 continue;
899 }
2eedd4ed 900 if(fDoCF){
901 fCFManager->GetParticleContainer()->Fill(containerInput,kStepdEdxElectronselection); // for CF
902 }
903
904 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PiPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion));
905 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PiMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion));
906 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_MuPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kMuon));
907 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_MuMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kMuon));
908
909
910 if( fCurrentPositiveESDTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentPositiveESDTrack->P()<fPIDMaxPnSigmaAbovePionLine ){
911 if(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
912 fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
913 fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
914 // iResult=kFALSE;
915 if(fHistograms != NULL){
916 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
3b77b2d1 917 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
2eedd4ed 918 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
919 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
920 //fUpdateV0AlreadyCalled = kTRUE;
921 }
922 fCurrentV0IndexNumber++;
923 continue;
924 }
925 }
926
927 if( fCurrentNegativeESDTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentNegativeESDTrack->P()<fPIDMaxPnSigmaAbovePionLine){
928 if(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
929 fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
930 fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
931 // iResult=kFALSE;
932 if(fHistograms != NULL){
933 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
3b77b2d1 934 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
2eedd4ed 935 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
936 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
937 //fUpdateV0AlreadyCalled = kTRUE;
938 }
939 fCurrentV0IndexNumber++;
940 continue;
941 }
942 }
943
944 // High Pt
945 if( fCurrentPositiveESDTrack->P()>fPIDMaxPnSigmaAbovePionLine ){
946 if(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
947 fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
948 fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
949 // iResult=kFALSE;
950 if(fHistograms != NULL){
951 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
3b77b2d1 952 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
2eedd4ed 953 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
954 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
955 //fUpdateV0AlreadyCalled = kTRUE;
956 }
957 fCurrentV0IndexNumber++;
958 continue;
959 }
960 }
961
962 if( fCurrentNegativeESDTrack->P()>fPIDMaxPnSigmaAbovePionLine){
963 if(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
964 fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
965 fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
966 // iResult=kFALSE;
967 if(fHistograms != NULL){
968 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
3b77b2d1 969 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
2eedd4ed 970 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
971 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
972 //fUpdateV0AlreadyCalled = kTRUE;
973 }
974 fCurrentV0IndexNumber++;
975 continue;
976 }
977 }
978
979
980 if(fDoCF){
981 fCFManager->GetParticleContainer()->Fill(containerInput,kStepdEdxPionrejection); // for CF
982 }
983
14d7f2f9 984 }
2eedd4ed 985
986 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_KPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon));
987 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_KMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon));
988
989 if(fDoKaonRejectionLowP == kTRUE){
990 if( fCurrentNegativeESDTrack->P()<fPIDMinPKaonRejectionLowP ){
991 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
992 if(fHistograms != NULL){
993 fHistograms->FillHistogram("ESD_CutKaonRejectionLowP_InvMass",GetMotherCandidateMass());
994 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
995 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
996 //fUpdateV0AlreadyCalled = kTRUE;
997 }
998 fCurrentV0IndexNumber++;
999 continue;
1000 }
1001 }
1002 if( fCurrentPositiveESDTrack->P()<fPIDMinPKaonRejectionLowP ){
1003 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
1004 if(fHistograms != NULL){
1005 fHistograms->FillHistogram("ESD_CutKaonRejectionLowP_InvMass",GetMotherCandidateMass());
1006 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1007 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1008 //fUpdateV0AlreadyCalled = kTRUE;
1009 }
1010 fCurrentV0IndexNumber++;
1011 continue;
1012 }
1013 }
1014 }
1015
1016 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton));
1017 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon));
1018
1019 if(fDoProtonRejectionLowP == kTRUE){
1020 if( fCurrentNegativeESDTrack->P()<fPIDMinPProtonRejectionLowP){
1021 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
1022 if(fHistograms != NULL){
1023 fHistograms->FillHistogram("ESD_CutProtonRejectionLowP_InvMass",GetMotherCandidateMass());
1024 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1025 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1026 //fUpdateV0AlreadyCalled = kTRUE;
1027 }
1028 fCurrentV0IndexNumber++;
1029 continue;
1030 }
1031 }
1032 if( fCurrentPositiveESDTrack->P()<fPIDMinPProtonRejectionLowP ){
1033 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
1034 if(fHistograms != NULL){
1035 fHistograms->FillHistogram("ESD_CutProtonRejectionLowP_InvMass",GetMotherCandidateMass());
1036 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1037 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1038 //fUpdateV0AlreadyCalled = kTRUE;
1039 }
1040 fCurrentV0IndexNumber++;
1041 continue;
1042 }
1043 }
1044 }
1045
1046 if(fDoPionRejectionLowP == kTRUE){
1047 if( fCurrentNegativeESDTrack->P()<fPIDMinPPionRejectionLowP ){
1048 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
1049 if(fHistograms != NULL){
1050 fHistograms->FillHistogram("ESD_CutPionRejectionLowP_InvMass",GetMotherCandidateMass());
1051 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1052 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1053 //fUpdateV0AlreadyCalled = kTRUE;
1054 }
1055 fCurrentV0IndexNumber++;
1056 continue;
1057 }
1058 }
1059 if( fCurrentPositiveESDTrack->P()<fPIDMinPPionRejectionLowP ){
1060 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
1061 if(fHistograms != NULL){
1062 fHistograms->FillHistogram("ESD_CutPionRejectionLowP_InvMass",GetMotherCandidateMass());
1063 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1064 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1065 //fUpdateV0AlreadyCalled = kTRUE;
1066 }
1067 fCurrentV0IndexNumber++;
1068 continue;
1069 }
1070 }
1071 }
1072
1073
3b77b2d1 1074 Double_t psiPair = -1;
1075 psiPair = GetPsiPair(fCurrentV0);
1076
1077 if(psiPair > fPsiPairCut){
1078 if(fHistograms != NULL ){
1079 fHistograms->FillHistogram("ESD_CutPsiPair_InvMass",GetMotherCandidateMass());
1080 fHistograms->FillHistogram("ESD_CutPsiPair_Pt",GetMotherCandidatePt());
1081 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1082 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1083 // fUpdateV0AlreadyCalled = kTRUE;
1084 }
1085
1086 fCurrentV0IndexNumber++;
1087 continue;
1088 }
1089
1090
1091 Double_t cosineOfPointingAngle = -1;
1092 cosineOfPointingAngle = GetV0CosineOfPointingAngle(fCurrentXValue,fCurrentYValue,fCurrentZValue);
1093 if(cosineOfPointingAngle < fCosinePointCut){
1094 if(fHistograms != NULL ){
1095 fHistograms->FillHistogram("ESD_CutCosinePoint_InvMass",GetMotherCandidateMass());
1096 fHistograms->FillHistogram("ESD_CutCosinePoint_Pt",GetMotherCandidatePt());
1097 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1098 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1099 // fUpdateV0AlreadyCalled = kTRUE;
1100 }
1101 fCurrentV0IndexNumber++;
1102 continue;
1103 }
1104
1105
2eedd4ed 1106 if( fDoTOFsigmaCut == kTRUE ){ // RRnewTOF start /////////////////////////////////////////////////////////////////////////////
1107 Bool_t PosTrackNotTOFelec = kFALSE;
1108 Bool_t NegTrackNotTOFelec = kFALSE;
1109 if( (fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTOFpid) && !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTOFmismatch) ){
1110 Double_t t0pos = fgESDpid->GetTOFResponse().GetStartTime(fCurrentPositiveESDTrack->P());
1111 Double_t fnSigmaPos = fgESDpid->NumberOfSigmasTOF(fCurrentPositiveESDTrack, AliPID::kElectron, t0pos);
1112 if( (fnSigmaPos>fTofPIDnSigmaAboveElectronLine) || (fnSigmaPos<fTofPIDnSigmaBelowElectronLine) ) PosTrackNotTOFelec = kTRUE;
1113 }
1114 if( (fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTOFpid) && !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTOFmismatch) ){
1115 Double_t t0neg = fgESDpid->GetTOFResponse().GetStartTime(fCurrentNegativeESDTrack->P());
1116 Double_t fnSigmaNeg = fgESDpid->NumberOfSigmasTOF(fCurrentNegativeESDTrack, AliPID::kElectron, t0neg);
1117 if( (fnSigmaNeg>fTofPIDnSigmaAboveElectronLine) || (fnSigmaNeg<fTofPIDnSigmaBelowElectronLine) ) NegTrackNotTOFelec = kTRUE;
1118 }
1119 if( (PosTrackNotTOFelec==kTRUE) || (NegTrackNotTOFelec==kTRUE) ){
14d7f2f9 1120 if(fHistograms != NULL){
2eedd4ed 1121 fHistograms->FillHistogram("ESD_CutTOFsigmaElec_InvMass",GetMotherCandidateMass());
3b77b2d1 1122 fHistograms->FillHistogram("ESD_CutTOFsigmaElec_Pt",GetMotherCandidatePt());
14d7f2f9 1123 }
1124 fCurrentV0IndexNumber++;
1125 continue;
1126 }
2eedd4ed 1127 } /////////////////////////////// RRnewTOF end ///////////////////////////////////////////////////////////////////////////////
1128
1129
1130 // Gamma selection based on QT from Armenteros
1131 if(fDoQtGammaSelection == kTRUE){ // RRnew start : apply different qT-cut above/below
1132 if(fDoHighPtQtGammaSelection){
1133 if(GetMotherCandidatePt() < fPtBorderForQt){
3b77b2d1 1134 if(armenterosQtAlpha[0]>fQtMax){
2eedd4ed 1135 if(fHistograms != NULL){
1136 fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
3b77b2d1 1137 fHistograms->FillHistogram("ESD_CutQt_Pt",GetMotherCandidatePt());
2eedd4ed 1138 }
1139 fCurrentV0IndexNumber++;
1140 continue;
1141 }
1142 } else {
3b77b2d1 1143 if(armenterosQtAlpha[0]>fHighPtQtMax) {
2eedd4ed 1144 if(fHistograms != NULL){
1145 fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
3b77b2d1 1146 fHistograms->FillHistogram("ESD_CutQt_Pt",GetMotherCandidatePt());
2eedd4ed 1147 }
1148 fCurrentV0IndexNumber++;
1149 continue;
1150 }
1151 }
1152 } else {
3b77b2d1 1153 if(armenterosQtAlpha[0]>fQtMax){
2eedd4ed 1154 if(fHistograms != NULL){
1155 fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
3b77b2d1 1156 fHistograms->FillHistogram("ESD_CutQt_Pt",GetMotherCandidatePt());
2eedd4ed 1157 }
1158 fCurrentV0IndexNumber++;
1159 continue;
1160 }
1161 }
1162 } // RRnew end
1163
1164 if(fDoPhotonAsymmetryCut == kTRUE){
1165 if( fNegativeTrackLorentzVector->P()>fMinPPhotonAsymmetryCut ){
1166 Double_t trackNegAsy=0;
1167 if (fCurrentMotherKFCandidate->GetP()!=0.){
1168 trackNegAsy= fNegativeTrackLorentzVector->P()/fMotherCandidateLorentzVector->P();
1169 }
1170 if( trackNegAsy<fMinPhotonAsymmetry ||trackNegAsy>(1.- fMinPhotonAsymmetry)){
1171 if(fHistograms != NULL){
1172 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_InvMass",GetMotherCandidateMass());
3b77b2d1 1173 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_Pt",GetMotherCandidatePt());
2eedd4ed 1174 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1175 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1176 //fUpdateV0AlreadyCalled = kTRUE;
1177 }
1178 fCurrentV0IndexNumber++;
1179 continue;
1180 }
1181 }
1182
1183 if( fPositiveTrackLorentzVector->P()>fMinPPhotonAsymmetryCut ){
1184 Double_t trackPosAsy=0;
1185 if (fCurrentMotherKFCandidate->GetP()!=0.){
1186 trackPosAsy= fPositiveTrackLorentzVector->P()/fMotherCandidateLorentzVector->P();
1187 }
1188 if( trackPosAsy<fMinPhotonAsymmetry ||trackPosAsy>(1.- fMinPhotonAsymmetry)){
1189 if(fHistograms != NULL){
1190 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_InvMass",GetMotherCandidateMass());
3b77b2d1 1191 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_Pt",GetMotherCandidatePt());
2eedd4ed 1192 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1193 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1194 //fUpdateV0AlreadyCalled = kTRUE;
1195 }
1196 fCurrentV0IndexNumber++;
1197 continue;
1198 }
1199 }
14d7f2f9 1200 }
2eedd4ed 1201 //checks if we have a prim vertex
1202 //if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
1203 if(GetNumberOfContributorsVtx()<=0) {
14d7f2f9 1204 if(fHistograms != NULL){
2eedd4ed 1205 fHistograms->FillHistogram("ESD_CutNContributors_InvMass",GetMotherCandidateMass());
3b77b2d1 1206 fHistograms->FillHistogram("ESD_CutNContributors_Pt",GetMotherCandidatePt());
14d7f2f9 1207 }
1208 fCurrentV0IndexNumber++;
1209 continue;
2eedd4ed 1210 }
1211 if(fDoCF){
1212 fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors); // for CF
1213 }
a0b94e5c 1214
2eedd4ed 1215 //Check the pid probability
1216 if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
1217 if(fHistograms != NULL){
1218 fHistograms->FillHistogram("ESD_CutPIDProb_InvMass",GetMotherCandidateMass());
3b77b2d1 1219 fHistograms->FillHistogram("ESD_CutPIDProb_Pt",GetMotherCandidatePt());
2eedd4ed 1220 }
1221 fCurrentV0IndexNumber++;
1222 continue;
1223 }
1224 if(fDoCF){
1225 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID); // for CF
1226 }
a0b94e5c 1227
a0b94e5c 1228
2eedd4ed 1229 /* Moved further up so corr framework can work
1230 if(UpdateV0Information() == kFALSE){
1231 fCurrentV0IndexNumber++;
1232 continue;
1233 }
1234 */
1235 if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC || fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
1236 if(fHistograms != NULL){
1237 fHistograms->FillHistogram("ESD_CutMinNClsTPC_InvMass",GetMotherCandidateMass());
3b77b2d1 1238 fHistograms->FillHistogram("ESD_CutMinNClsTPC_Pt",GetMotherCandidatePt());
2eedd4ed 1239 }
1240 fCurrentV0IndexNumber++;
1241 continue;
1242 }
1243 if(fDoCF){
1244 fCFManager->GetParticleContainer()->Fill(containerInput,kStepMinClsTPC); // for CF
1245 }
1246 Double_t negclsToF = 0.;
1247 if (!fUseCorrectedTPCClsInfo ){
1248 if(fCurrentNegativeESDTrack->GetTPCNclsF()!=0 ){
1249 negclsToF = (Double_t)fCurrentNegativeESDTrack->GetNcls(1)/(Double_t)fCurrentNegativeESDTrack->GetTPCNclsF();
1250 }
1251 } else {
1252 negclsToF = fCurrentNegativeESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
1253 }
1254
1255 Double_t posclsToF = 0.;
1256 if (!fUseCorrectedTPCClsInfo ){
1257 if(fCurrentPositiveESDTrack->GetTPCNclsF()!=0 ){
1258 posclsToF = (Double_t)fCurrentPositiveESDTrack->GetNcls(1)/(Double_t)fCurrentPositiveESDTrack->GetTPCNclsF();
1259 }
1260 }else{
1261 posclsToF = fCurrentPositiveESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
1262 }
1263
1264 if( negclsToF < fMinClsTPCToF || posclsToF < fMinClsTPCToF ){
1265 if(fHistograms != NULL){
1266 fHistograms->FillHistogram("ESD_CutMinNClsTPCToF_InvMass",GetMotherCandidateMass());
3b77b2d1 1267 fHistograms->FillHistogram("ESD_CutMinNClsTPCToF_Pt",GetMotherCandidatePt());
2eedd4ed 1268 }
1269 fCurrentV0IndexNumber++;
1270 continue;
1271 }
bd6d9fa3 1272
1273
9640a3d1 1274
a0b94e5c 1275
2eedd4ed 1276 if(fUseKFParticle){
48682642 1277
1278
2eedd4ed 1279 if( fCurrentNegativeKFParticle->GetPt()< fSinglePtCut || fCurrentPositiveKFParticle->GetPt()< fSinglePtCut){
1280 if(fHistograms != NULL){
1281 fHistograms->FillHistogram("ESD_CutSinglePt_InvMass",GetMotherCandidateMass());
3b77b2d1 1282 fHistograms->FillHistogram("ESD_CutSinglePt_Pt",GetMotherCandidatePt());
2eedd4ed 1283 }
1284 fCurrentV0IndexNumber++;
1285 continue;
1286 }
1287 if(fDoCF){
1288 fCFManager->GetParticleContainer()->Fill(containerInput,kStepSinglePt); // for CF
1289 }
48682642 1290
1291
2eedd4ed 1292 if(fCurrentMotherKFCandidate->GetNDF()<=0){
1293 if(fHistograms != NULL){
1294 fHistograms->FillHistogram("ESD_CutNDF_InvMass",GetMotherCandidateMass());
3b77b2d1 1295 fHistograms->FillHistogram("ESD_CutNDF_Pt",GetMotherCandidatePt());
2eedd4ed 1296 }
1297 fCurrentV0IndexNumber++;
1298 continue;
1299 }
1300 if(fDoCF){
1301 fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF); // for CF
1302 }
a0b94e5c 1303
2eedd4ed 1304 Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
1305 if(chi2V0 > fChi2CutConversion || chi2V0 <=0){
1306 if(fHistograms != NULL){
1307 fHistograms->FillHistogram("ESD_CutChi2_InvMass",GetMotherCandidateMass());
3b77b2d1 1308 fHistograms->FillHistogram("ESD_CutChi2_Pt",GetMotherCandidatePt());
2eedd4ed 1309 }
1310 fCurrentV0IndexNumber++;
1311 continue;
1312 }
1313 if(fDoCF){
1314 fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2); // for CF
1315 }
1316
1317 if(fDoCF){
1318 fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta); // for CF
1319 }
a0b94e5c 1320
2eedd4ed 1321 if(fMotherCandidateLorentzVector->Pt()<fPtCut){
1322 if(fHistograms != NULL){
1323 fHistograms->FillHistogram("ESD_CutPt_InvMass",GetMotherCandidateMass());
3b77b2d1 1324 fHistograms->FillHistogram("ESD_CutPt_Pt",GetMotherCandidatePt());
2eedd4ed 1325 }
1326 fCurrentV0IndexNumber++;
1327 continue;
1328 }
1329 if(fDoCF){
1330 fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt); // for CF
1331 }
a0b94e5c 1332
2eedd4ed 1333 }
1334 else if(fUseESDTrack){
1335 //TODO
1336 }
9640a3d1 1337
2eedd4ed 1338 if(fHistograms != NULL){
1339 fHistograms->FillHistogram("ESD_GoodV0s_InvMass",GetMotherCandidateMass());
1340 }
9640a3d1 1341
2eedd4ed 1342 // fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
5e55d806 1343
2eedd4ed 1344 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1345 fCurrentMotherKFCandidate->E()=fCurrentMotherKFCandidate->GetP();
1346 }
2ea0dd4a 1347
2eedd4ed 1348 if(fDoMC&& fUseMCPSmearing>0){
1349 SmearKFParticle(fCurrentMotherKFCandidate);
2ea0dd4a 1350
2eedd4ed 1351 }
2ea0dd4a 1352
2eedd4ed 1353 new((*fCurrentEventGoodV0s)[fCurrentEventGoodV0s->GetEntriesFast()]) AliKFParticle(*fCurrentMotherKFCandidate);
1354 fV0Pindex.push_back(fCurrentV0->GetPindex());
1355 fV0Nindex.push_back(fCurrentV0->GetNindex());
5e55d806 1356
2eedd4ed 1357 iResult=kTRUE;//means we have a v0 who survived all the cuts applied
a0b94e5c 1358
2eedd4ed 1359 fNumberOfGoodV0s++;
5ce758b0 1360
2eedd4ed 1361 fCurrentV0IndexNumber++;
a0b94e5c 1362
2eedd4ed 1363 break;
1364 }
1365 return iResult;
a0b94e5c 1366}
1367
1368Bool_t AliV0Reader::UpdateV0Information(){
2eedd4ed 1369 //see header file for documentation
a0b94e5c 1370
3b77b2d1 1371 const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0);
1372 const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0);
1373
2eedd4ed 1374 Bool_t iResult=kTRUE; // for taking out not refitted, kinks and like sign tracks
a0b94e5c 1375
2eedd4ed 1376 Bool_t switchTracks = kFALSE;
a0b94e5c 1377
2eedd4ed 1378 fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1379 fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1380
1381
1382 if(fCurrentPositiveESDTrack->GetSign() == -1 && fCurrentNegativeESDTrack->GetSign() == 1){ // switch wrong signed tracks
1383 fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1384 fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1385 switchTracks = kTRUE;
1386 }
1387
a0b94e5c 1388
2eedd4ed 1389 if(fCurrentNegativeKFParticle != NULL){
1390 delete fCurrentNegativeKFParticle;
1391 }
1392 if(switchTracks == kFALSE){
1393 fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);
1394 }
1395 else{
1396 fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fNegativeTrackPID);
1397 }
a0b94e5c 1398
2eedd4ed 1399 if(fCurrentPositiveKFParticle != NULL){
1400 delete fCurrentPositiveKFParticle;
1401 }
1402 if(switchTracks == kFALSE){
1403 fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);
1404 }
1405 else{
1406 fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fPositiveTrackPID);
1407 }
1408
1409 if(fCurrentMotherKFCandidate != NULL){
1410 delete fCurrentMotherKFCandidate;
1411 }
1412
1413 if(fUseConstructGamma==kTRUE){
1414 fCurrentMotherKFCandidate = new AliKFParticle;//(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1415 fCurrentMotherKFCandidate->ConstructGamma(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1416 }else{
1417 fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1418 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1419 fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
1420 }
1421 }
1422 if(fUseImprovedVertex == kTRUE){
1423 AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
1424 primaryVertexImproved+=*fCurrentMotherKFCandidate;
1425 fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
1426 }
a0b94e5c 1427
2eedd4ed 1428 fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);
ebcfaa7e 1429
2eedd4ed 1430 if(fNegativeTrackLorentzVector != NULL){
1431 delete fNegativeTrackLorentzVector;
1432 }
1433 if(fUseKFParticle){
1434 fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());
1435 }
1436 else { //if(fUseESDTrack){
1437 fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());
1438 }
a0b94e5c 1439
2eedd4ed 1440 if(fPositiveTrackLorentzVector != NULL){
1441 delete fPositiveTrackLorentzVector;
1442 }
1443 if(fUseKFParticle){
1444 fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());
1445 }
1446 else { // if(fUseESDTrack){ fPositiveTrackLorentzVector must be reinitialized, so assuming use ESD if not kfparticle Svein.
1447 fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());
1448 }
a0b94e5c 1449
2eedd4ed 1450 if(fMotherCandidateLorentzVector != NULL){
1451 delete fMotherCandidateLorentzVector;
1452 }
239a56d1 1453
2eedd4ed 1454 fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
a0b94e5c 1455
2eedd4ed 1456 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1457 fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.);
1458 }
1459
a0b94e5c 1460
2eedd4ed 1461 if(fDoMC == kTRUE){
1462 fMotherMCParticle= NULL;
1463 if(switchTracks == kFALSE){
1464 fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
1465 fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
1466 }else{
1467 fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
1468 fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
1469 }
1470
1471 if(fPositiveMCParticle->GetMother(0)>-1){
1472 fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));
1473 }
1474 }
1475
1476
1477
a0b94e5c 1478
a0b94e5c 1479
2eedd4ed 1480 // for CF
1481// Double_t containerInput[3];
1482// if(fDoCF){
1483// containerInput[0] = GetMotherCandidatePt();
1484// containerInput[1] = GetMotherCandidateEta();
1485// containerInput[2] = GetMotherCandidateMass();
1486
1487// fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign); // for CF
1488// fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit); // for CF
1489// fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks); // for CF
1490// }
1491
cb90a330 1492
2eedd4ed 1493 if(fUseOwnXYZCalculation == kFALSE){
1494 if(fUseConstructGamma == kFALSE){
1495 fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
1496 }else{
1497 fCurrentXValue=GetMotherCandidateKFCombination()->GetX();
1498 fCurrentYValue=GetMotherCandidateKFCombination()->GetY();
1499 fCurrentZValue=GetMotherCandidateKFCombination()->GetZ();
1500 }
1501 }
1502 else{
3b77b2d1 1503 Double_t convpos[3]={0,0,0};
1504 GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos);
1505// fCurrentMotherKF->SetConversionPoint(convpos);
1506
1507// Double_t convpos[2];
1508// convpos[0]=0;
1509// convpos[1]=0;
1510//
1511// GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
1512//
2eedd4ed 1513 fCurrentXValue = convpos[0];
1514 fCurrentYValue = convpos[1];
3b77b2d1 1515 fCurrentZValue = convpos[2];
2eedd4ed 1516 }
2eedd4ed 1517 fUpdateV0AlreadyCalled = kTRUE;
1518
1519 return iResult;
a0b94e5c 1520}
1521
1522
1523
1524Bool_t AliV0Reader::HasSameMCMother(){
2eedd4ed 1525 //see header file for documentation
a0b94e5c 1526
2eedd4ed 1527 Bool_t iResult = kFALSE;
1528 if(fDoMC == kTRUE){
1529 if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){
1530 if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))
a0b94e5c 1531 if(fMotherMCParticle){
2eedd4ed 1532 iResult = kTRUE;
1533 }
1534 }
a0b94e5c 1535 }
2eedd4ed 1536 return iResult;
a0b94e5c 1537}
1538
1539Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){
2eedd4ed 1540 //see header file for documentation
a0b94e5c 1541
2eedd4ed 1542 Bool_t iResult=kFALSE;
a0b94e5c 1543
2eedd4ed 1544 // Double_t *posProbArray = new Double_t[10];
1545 // Double_t *negProbArray = new Double_t[10];
1546 //-AM The TPCpid method expects an array of length kSPECIES that is 5 not 10
1547
1548 Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1549 Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1550
1551 AliESDtrack* negTrack = GetNegativeESDTrack();
1552 AliESDtrack* posTrack = GetPositiveESDTrack();
1553 //fESDEvent->GetTrack(fCurrentV0->GetNindex());
1554 //fESDEvent->GetTrack(fCurrentV0->GetPindex());
1555 //-AM for switchtracks==true the above is a bug
1556
1557 if(negProbArray && posProbArray){
1558
1559 negTrack->GetTPCpid(negProbArray);
1560 posTrack->GetTPCpid(posProbArray);
1561
1562 // if(negProbArray != NULL && posProbArray != NULL){ // this is not allowed anymore for some reason(RC19)
1563 if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){
1564 iResult=kTRUE;
1565 }
1566 }
1567 delete [] posProbArray;
1568 delete [] negProbArray;
1569 return iResult;
a0b94e5c 1570}
1571
1572void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){
2eedd4ed 1573 // see header file for documentation
1574
1575 //Double_t *posProbArray = new Double_t[10];
1576 // Double_t *negProbArray = new Double_t[10];
1577 //-AM The TPCpid method expects an array of length kSPECIES that is 5 not 10
1578 Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1579 Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1580
1581// AliESDtrack* negTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1582// AliESDtrack* posTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1583 //-AM for switchtracks the above is a bug
1584 AliESDtrack* negTrack = GetNegativeESDTrack();
1585 AliESDtrack* posTrack = GetPositiveESDTrack();
1586
1587 if(negProbArray && posProbArray){
1588 negTrack->GetTPCpid(negProbArray);
1589 posTrack->GetTPCpid(posProbArray);
1590
1591 // if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
1592 negPIDProb = negProbArray[GetSpeciesIndex(-1)];
1593 posPIDProb = posProbArray[GetSpeciesIndex(1)];
1594 }
1595 delete [] posProbArray;
1596 delete [] negProbArray;
a0b94e5c 1597}
239a56d1 1598
9c1cb6f7 1599void AliV0Reader::GetPIDProbabilityMuonPion(Double_t &negPIDProb,Double_t & posPIDProb){
2eedd4ed 1600 // see header file for documentation
9c1cb6f7 1601
1602
2eedd4ed 1603 Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1604 Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
9c1cb6f7 1605
2eedd4ed 1606 // AliESDtrack* negTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1607 // AliESDtrack* posTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1608 //-AM for switchtracks the above is a bug
9c1cb6f7 1609
2eedd4ed 1610 AliESDtrack* negTrack = GetNegativeESDTrack();
1611 AliESDtrack* posTrack = GetPositiveESDTrack();
9c1cb6f7 1612
2eedd4ed 1613 if(negProbArray && posProbArray){
1614 negTrack->GetTPCpid(negProbArray);
1615 posTrack->GetTPCpid(posProbArray);
1616
1617 // if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
239a56d1 1618
2eedd4ed 1619 negPIDProb = negProbArray[1]+negProbArray[2];
1620 posPIDProb = posProbArray[1]+posProbArray[2];
1621 }
1622 delete [] posProbArray;
1623 delete [] negProbArray;
9c1cb6f7 1624}
a0b94e5c 1625
1626void AliV0Reader::UpdateEventByEventData(){
2eedd4ed 1627 //see header file for documentation
1628 if(fCurrentEventGoodV0s->GetEntriesFast() >0 ){
1629 if(fCalculateBackground){
1630 if(fUseChargedTrackMultiplicityForBG == kTRUE){
5ce758b0 1631 fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
1632 //filling z and multiplicity histograms
1633 fHistograms->FillHistogram("ESD_Z_distribution",fESDEvent->GetPrimaryVertex()->GetZ());
1634 fHistograms->FillHistogram("ESD_multiplicity_distribution",CountESDTracks());
1635 fHistograms->FillHistogram("ESD_ZvsMultiplicity",fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
2eedd4ed 1636 }
1637 else{ // means we use #V0s for multiplicity
5ce758b0 1638 fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
1639 //filling z and multiplicity histograms
1640 fHistograms->FillHistogram("ESD_Z_distribution",fESDEvent->GetPrimaryVertex()->GetZ());
1641 fHistograms->FillHistogram("ESD_multiplicity_distribution",fNumberOfGoodV0s);
1642 fHistograms->FillHistogram("ESD_ZvsMultiplicity",fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
2eedd4ed 1643 }
1644 }
1645 }
1646 fCurrentEventGoodV0s->Delete();
1647 fCurrentV0IndexNumber=0;
1648 fNumberOfESDTracks=0;
1649
1650 fV0Pindex.clear();
1651 fV0Nindex.clear();
1652
7f3c7cc6 1653
1654
2eedd4ed 1655 // fBGEventHandler->PrintBGArray(); // for debugging
a0b94e5c 1656}
1657
1658
1659Double_t AliV0Reader::GetNegativeTrackPhi() const{
2eedd4ed 1660 //see header file for documentation
a0b94e5c 1661
2eedd4ed 1662 Double_t offset=0;
1663 if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
1664 offset = -2*TMath::Pi();
1665 }
1666 return fNegativeTrackLorentzVector->Phi()+offset;
a0b94e5c 1667}
1668
1669Double_t AliV0Reader::GetPositiveTrackPhi() const{
2eedd4ed 1670 //see header file for documentation
a0b94e5c 1671
2eedd4ed 1672 Double_t offset=0;
1673 if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
1674 offset = -2*TMath::Pi();
1675 }
1676 return fPositiveTrackLorentzVector->Phi()+offset;
a0b94e5c 1677}
1678
1679Double_t AliV0Reader::GetMotherCandidatePhi() const{
2eedd4ed 1680 //see header file for documentation
a0b94e5c 1681
2eedd4ed 1682 Double_t offset=0;
1683 if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
1684 offset = -2*TMath::Pi();
1685 }
1686 return fMotherCandidateLorentzVector->Phi()+offset;
a0b94e5c 1687}
1688
1689
1690Double_t AliV0Reader::GetMotherCandidateRapidity() const{
2eedd4ed 1691 //see header file for documentation
a0b94e5c 1692
2eedd4ed 1693 Double_t rapidity=0;
1694 if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
1695 else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
1696 return rapidity;
a0b94e5c 1697
1698}
1699
1700
1701
1702
1703
1704Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
2eedd4ed 1705 //see header file for documentation
a0b94e5c 1706
2eedd4ed 1707 Int_t iResult = 10; // Unknown particle
a0b94e5c 1708
2eedd4ed 1709 if(chargeOfTrack==-1){ //negative track
1710 switch(abs(fNegativeTrackPID)){
1711 case 11: //electron
1712 iResult = 0;
1713 break;
1714 case 13: //muon
1715 iResult = 1;
1716 break;
1717 case 211: //pion
1718 iResult = 2;
1719 break;
1720 case 321: //kaon
1721 iResult = 3;
1722 break;
1723 case 2212: //proton
1724 iResult = 4;
1725 break;
1726 case 22: //photon
1727 iResult = 5;
1728 break;
1729 case 111: //pi0
1730 iResult = 6;
1731 break;
1732 case 2112: //neutron
1733 iResult = 7;
1734 break;
1735 case 311: //K0
1736 iResult = 8;
1737 break;
a0b94e5c 1738
2eedd4ed 1739 //Put in here for kSPECIES::kEleCon ????
1740 }
1741 }
1742 else if(chargeOfTrack==1){ //positive track
1743 switch(abs(fPositiveTrackPID)){
1744 case 11: //electron
1745 iResult = 0;
1746 break;
1747 case 13: //muon
1748 iResult = 1;
1749 break;
1750 case 211: //pion
1751 iResult = 2;
1752 break;
1753 case 321: //kaon
1754 iResult = 3;
1755 break;
1756 case 2212: //proton
1757 iResult = 4;
1758 break;
1759 case 22: //photon
1760 iResult = 5;
1761 break;
1762 case 111: //pi0
1763 iResult = 6;
1764 break;
1765 case 2112: //neutron
1766 iResult = 7;
1767 break;
1768 case 311: //K0
1769 iResult = 8;
1770 break;
a0b94e5c 1771
2eedd4ed 1772 //Put in here for kSPECIES::kEleCon ????
1773 }
1774 }
1775 else{
1776 //Wrong parameter.. Print warning
1777 }
1778 return iResult;
a0b94e5c 1779}
1780
3b77b2d1 1781//_____________________________________________________________________________________________________________________
1782Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
2eedd4ed 1783 // see header file for documentation
2eedd4ed 1784
1785 Double_t helix[6];
1786 track->GetHelixParameters(helix,b);
1787
1788 Double_t xpos = helix[5];
1789 Double_t ypos = helix[0];
1790 Double_t radius = TMath::Abs(1./helix[4]);
1791 Double_t phi = helix[2];
1792
1793 if(phi < 0){
3b77b2d1 1794 phi = phi + 2*TMath::Pi();
2eedd4ed 1795 }
1796
3b77b2d1 1797 phi -= TMath::Pi()/2.;
2eedd4ed 1798 Double_t xpoint = radius * TMath::Cos(phi);
1799 Double_t ypoint = radius * TMath::Sin(phi);
1800
1801 if(b<0){
3b77b2d1 1802 if(charge > 0){
1803 xpoint = - xpoint;
1804 ypoint = - ypoint;
1805 }
1806
1807 if(charge < 0){
1808 xpoint = xpoint;
1809 ypoint = ypoint;
1810 }
2eedd4ed 1811 }
3b77b2d1 1812 if(b>0){
1813 if(charge > 0){
1814 xpoint = xpoint;
1815 ypoint = ypoint;
1816 }
2eedd4ed 1817
3b77b2d1 1818 if(charge < 0){
1819 xpoint = - xpoint;
1820 ypoint = - ypoint;
1821 }
1822 }
2eedd4ed 1823 center[0] = xpos + xpoint;
1824 center[1] = ypos + ypoint;
1825
1826 return 1;
a0b94e5c 1827}
1828
3b77b2d1 1829//_________________________________________________________________________________________________________
1830// Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
1831// // Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1832// // see header file for documentation
1833//
1834// Double_t pi = 3.14159265358979323846;
1835//
1836// Double_t helix[6];
1837// track->GetHelixParameters(helix,b);
1838//
1839// Double_t xpos = helix[5];
1840// Double_t ypos = helix[0];
1841// Double_t radius = TMath::Abs(1./helix[4]);
1842// Double_t phi = helix[2];
1843//
1844// if(phi < 0){
1845// phi = phi + 2*pi;
1846// }
1847//
1848// phi -= pi/2.;
1849// Double_t xpoint = radius * TMath::Cos(phi);
1850// Double_t ypoint = radius * TMath::Sin(phi);
1851//
1852// if(b<0){
1853// if(charge > 0){
1854// xpoint = - xpoint;
1855// ypoint = - ypoint;
1856// }
1857//
1858// if(charge < 0){
1859// xpoint = xpoint;
1860// ypoint = ypoint;
1861// }
1862// }
1863// if(b>0){
1864// if(charge > 0){
1865// xpoint = xpoint;
1866// ypoint = ypoint;
1867// }
1868//
1869// if(charge < 0){
1870// xpoint = - xpoint;
1871// ypoint = - ypoint;
1872// }
1873// }
1874// center[0] = xpos + xpoint;
1875// center[1] = ypos + ypoint;
1876//
1877// return 1;
1878// }
1879
1880//____________________________________________________________________________________________________________________________
1881Bool_t AliV0Reader::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3]){
1882
1883 if(!pparam||!nparam)return kFALSE;
2eedd4ed 1884
1885 Double_t helixcenterpos[2];
1886 GetHelixCenter(pparam,GetMagneticField(),pparam->Charge(),helixcenterpos);
3b77b2d1 1887
2eedd4ed 1888 Double_t helixcenterneg[2];
1889 GetHelixCenter(nparam,GetMagneticField(),nparam->Charge(),helixcenterneg);
3b77b2d1 1890
2eedd4ed 1891 Double_t helixpos[6];
1892 pparam->GetHelixParameters(helixpos,GetMagneticField());
1893 Double_t posradius = TMath::Abs(1./helixpos[4]);
3b77b2d1 1894
2eedd4ed 1895 Double_t helixneg[6];
1896 nparam->GetHelixParameters(helixneg,GetMagneticField());
1897 Double_t negradius = TMath::Abs(1./helixneg[4]);
3b77b2d1 1898
1899 // Calculate xy-position
1900
2eedd4ed 1901 Double_t xpos = helixcenterpos[0];
1902 Double_t ypos = helixcenterpos[1];
1903 Double_t xneg = helixcenterneg[0];
1904 Double_t yneg = helixcenterneg[1];
3b77b2d1 1905
2eedd4ed 1906 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
3b77b2d1 1907 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
1908
1909
2eedd4ed 1910 // Calculate z-position
2eedd4ed 1911
3b77b2d1 1912 Double_t deltaXPos = convpos[0] - xpos;
1913 Double_t deltaYPos = convpos[1] - ypos;
2eedd4ed 1914
3b77b2d1 1915 Double_t deltaXNeg = convpos[0] - xneg;
1916 Double_t deltaYNeg = convpos[1] - yneg;
2eedd4ed 1917
3b77b2d1 1918 Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
1919 Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
2eedd4ed 1920
3b77b2d1 1921 Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
1922 Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
1923
1924 Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
1925 Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
1926
1927 Double_t b = fESDEvent->GetMagneticField();
1928
1929 AliExternalTrackParam p(*pparam);
1930 AliExternalTrackParam n(*nparam);
1931
1932 TVector2 vertexPos(vertexXPos,vertexYPos);
1933 TVector2 vertexNeg(vertexXNeg,vertexYNeg);
1934
1935 // Convert to local coordinate system
1936 vertexPos=vertexPos.Rotate(-p.GetAlpha());
1937 vertexNeg=vertexNeg.Rotate(-p.GetAlpha());
1938
1939 // Propagate Track Params to Vertex
1940 p.PropagateTo(vertexPos.X(),b);
1941 n.PropagateTo(vertexNeg.X(),b);
2eedd4ed 1942
3b77b2d1 1943 convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
2eedd4ed 1944
3b77b2d1 1945 return kTRUE;
1946}
1947//
1948// //__________________________________________________________________________________________________________
4a6157dc 1949Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
2eedd4ed 1950 //see header file for documentation
a0b94e5c 1951
2eedd4ed 1952 Double_t helixcenterpos[2];
1953 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
a0b94e5c 1954
2eedd4ed 1955 Double_t helixcenterneg[2];
1956 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
a0b94e5c 1957
2eedd4ed 1958 Double_t poshelix[6];
1959 ptrack->GetHelixParameters(poshelix,b);
1960 Double_t posradius = TMath::Abs(1./poshelix[4]);
a0b94e5c 1961
2eedd4ed 1962 Double_t neghelix[6];
1963 ntrack->GetHelixParameters(neghelix,b);
1964 Double_t negradius = TMath::Abs(1./neghelix[4]);
a0b94e5c 1965
2eedd4ed 1966 Double_t xpos = helixcenterpos[0];
1967 Double_t ypos = helixcenterpos[1];
1968 Double_t xneg = helixcenterneg[0];
1969 Double_t yneg = helixcenterneg[1];
a0b94e5c 1970
2eedd4ed 1971 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1972 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
a0b94e5c 1973
2eedd4ed 1974 return 1;
a0b94e5c 1975}
3b77b2d1 1976//
1977//
1978//
4a6157dc 1979Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
2eedd4ed 1980 //see header file for documentation
a0b94e5c 1981
2eedd4ed 1982 Double_t helixpos[6];
1983 ptrack->GetHelixParameters(helixpos,b);
a0b94e5c 1984
2eedd4ed 1985 Double_t helixneg[6];
1986 ntrack->GetHelixParameters(helixneg,b);
a0b94e5c 1987
2eedd4ed 1988 Double_t negtrackradius = TMath::Abs(1./helixneg[4]);
1989 Double_t postrackradius = TMath::Abs(1./helixpos[4]);
a0b94e5c 1990
2eedd4ed 1991 Double_t pi = 3.14159265358979323846;
a0b94e5c 1992
2eedd4ed 1993 Double_t convpos[2];
1994 GetConvPosXY(ptrack,ntrack,b,convpos);
a0b94e5c 1995
2eedd4ed 1996 Double_t convposx = convpos[0];
1997 Double_t convposy = convpos[1];
a0b94e5c 1998
2eedd4ed 1999 Double_t helixcenterpos[2];
2000 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
a0b94e5c 2001
2eedd4ed 2002 Double_t helixcenterneg[2];
2003 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
a0b94e5c 2004
2eedd4ed 2005 Double_t xpos = helixcenterpos[0];
2006 Double_t ypos = helixcenterpos[1];
2007 Double_t xneg = helixcenterneg[0];
2008 Double_t yneg = helixcenterneg[1];
a0b94e5c 2009
2eedd4ed 2010 Double_t deltaXPos = convposx - xpos;
2011 Double_t deltaYPos = convposy - ypos;
a0b94e5c 2012
2eedd4ed 2013 Double_t deltaXNeg = convposx - xneg;
2014 Double_t deltaYNeg = convposy - yneg;
a0b94e5c 2015
2eedd4ed 2016 Double_t alphaPos = pi + TMath::ATan2(-deltaYPos,-deltaXPos);
2017 Double_t alphaNeg = pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
a0b94e5c 2018
2eedd4ed 2019 Double_t vertexXNeg = xneg + TMath::Abs(negtrackradius)*
2020 TMath::Cos(alphaNeg);
2021 Double_t vertexYNeg = yneg + TMath::Abs(negtrackradius)*
2022 TMath::Sin(alphaNeg);
a0b94e5c 2023
2eedd4ed 2024 Double_t vertexXPos = xpos + TMath::Abs(postrackradius)*
2025 TMath::Cos(alphaPos);
2026 Double_t vertexYPos = ypos + TMath::Abs(postrackradius)*
2027 TMath::Sin(alphaPos);
a0b94e5c 2028
2eedd4ed 2029 Double_t x0neg = helixneg[5];
2030 Double_t y0neg = helixneg[0];
a0b94e5c 2031
2eedd4ed 2032 Double_t x0pos = helixpos[5];
2033 Double_t y0pos = helixpos[0];
a0b94e5c 2034
2eedd4ed 2035 Double_t dNeg = TMath::Sqrt((vertexXNeg - x0neg)*(vertexXNeg - x0neg)
2036 +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
a0b94e5c 2037
2eedd4ed 2038 Double_t dPos = TMath::Sqrt((vertexXPos - x0pos)*(vertexXPos - x0pos)
2039 +(vertexYPos - y0pos)*(vertexYPos - y0pos));
a0b94e5c 2040
2eedd4ed 2041 Double_t rNeg = TMath::Sqrt(negtrackradius*negtrackradius -
2042 dNeg*dNeg/4.);
a0b94e5c 2043
2eedd4ed 2044 Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
2045 dPos*dPos/4.);
a0b94e5c 2046
2eedd4ed 2047 Double_t deltabetaNeg = 2*(pi + TMath::ATan2(-dNeg/2.,-rNeg));
2048 Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
a0b94e5c 2049
2eedd4ed 2050 Double_t deltaUNeg = negtrackradius*deltabetaNeg;
2051 Double_t deltaUPos = postrackradius*deltabetaPos;
a0b94e5c 2052
2eedd4ed 2053 Double_t zphaseNeg = ntrack->GetZ() + deltaUNeg * ntrack->GetTgl();
2054 Double_t zphasePos = ptrack->GetZ() + deltaUPos * ptrack->GetTgl();
a0b94e5c 2055
2eedd4ed 2056 Double_t convposz = (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
a0b94e5c 2057
2eedd4ed 2058 return convposz;
a0b94e5c 2059}
3b77b2d1 2060//
f14266e7 2061AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t /*event*/) const{
2eedd4ed 2062 /*
2063 if(fUseChargedTrackMultiplicityForBG == kTRUE){
2064 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
2065 }
2066 else{ // means we use #v0s as multiplicity
2067 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
2068 }
2069 */
2070 return NULL;
037dc2db 2071}
2072
2073Int_t AliV0Reader::CountESDTracks(){
2eedd4ed 2074 // see header file for documentation
2075 if(fNumberOfESDTracks == 0){ // count the good esd tracks
2076 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2077 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2078 if(!curTrack){
037dc2db 2079 continue;
2eedd4ed 2080 }
2081 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
037dc2db 2082 fNumberOfESDTracks++;
2eedd4ed 2083 }
2084 }
2085 }
037dc2db 2086
2eedd4ed 2087 return fNumberOfESDTracks;
037dc2db 2088}
2089
2090Bool_t AliV0Reader::CheckIfPi0IsMother(Int_t label){
2eedd4ed 2091 // see headerfile for documentation
2092 Bool_t iResult=kFALSE;
2093 // cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2094 if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 111){
2095 iResult=kTRUE;
2096 }
2097 return iResult;
5e55d806 2098}
48682642 2099
10e3319b 2100Bool_t AliV0Reader::CheckIfEtaIsMother(Int_t label){
2eedd4ed 2101 // see headerfile for documentation
2102 Bool_t iResult=kFALSE;
2103 // cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2104 if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 221){
2105 iResult=kTRUE;
2106 }
2107 return iResult;
10e3319b 2108}
2109
48682642 2110
3b77b2d1 2111
2112
2113Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, Double_t armenterosQtAlpha[2], Int_t kfProductionMethod){
2114
2115 AliKFParticle PosParticle = *positiveKFParticle;
2116 AliKFParticle NegParticle = *negativeKFParticle;
2117 AliKFParticle Gamma;
2118 if(kfProductionMethod < 3)
2119 Gamma.ConstructGamma(PosParticle, NegParticle);
2120 else if(kfProductionMethod == 3){
2121 Gamma += PosParticle;
2122 Gamma += NegParticle;
2123 }
2124
2125 Double_t VertexGamma[3] = {Gamma.GetX(), Gamma.GetY(), Gamma.GetZ()};
2126 PosParticle.TransportToPoint(VertexGamma);
2127 NegParticle.TransportToPoint(VertexGamma);
2128
2129 AliKFParticle::GetArmenterosPodolanski(PosParticle, NegParticle, armenterosQtAlpha);
2130
2131 return 1;
2132}
2133
2134
2135
2136
2137Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2eedd4ed 2138 //see header file for documentation
48682642 2139
2eedd4ed 2140 TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
2141 TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
2142 TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
48682642 2143
2eedd4ed 2144 Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2145 Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2146
2147 Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2148 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2149
2150
2151 Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2152
3b77b2d1 2153 armenterosQtAlpha[0]=qt;
2154 armenterosQtAlpha[1]=alfa;
2eedd4ed 2155
2156 return 1;
2157
2158}
2159
3b77b2d1 2160Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliESDv0* v0, Double_t armenterosQtAlpha[2])
2161{ //see header file for documentation
2162
2163 Double_t mn[3] = {0,0,0};
2164 Double_t mp[3] = {0,0,0};
2165 Double_t mm[3] = {0,0,0};
2166
2167
2168 Int_t pIndex = 0, nIndex = 0;
2169 pIndex = v0->GetPindex();
2170 nIndex = v0->GetNindex();
2171
2172 AliESDtrack* d[2];
2173 d[0] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(pIndex));
2174 d[1] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(nIndex));
2175
2176 Int_t sign[2];
2177 sign[0] = (int)d[0]->GetSign();
2178 sign[1] = (int)d[1]->GetSign();
2179
2180 Bool_t correct = kFALSE;
2181
2182
2183 if(-1 == sign[0] && 1 == sign[1]){
2184 correct = kFALSE;
2185 }
2186 else{
2187 correct = kTRUE;
2188 }
2189
2190
2191 if(correct){
2192 v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2193 v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2194 }
2195 else{
2196 v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2197 v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2198 }
2199 v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
2200
2201 TVector3 vecN(mn[0],mn[1],mn[2]);
2202 TVector3 vecP(mp[0],mp[1],mp[2]);
2203 TVector3 vecM(mm[0],mm[1],mm[2]);
2204
2205 Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
2206 Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
2207
2208 Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
2209 ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
2210 Double_t qt = vecP.Mag()*sin(thetaP);
2211
2212 armenterosQtAlpha[0]=qt;
2213 armenterosQtAlpha[1]=alfa;
2214
2215 return 1;
2216
2217}
2218
2219
2220
2221Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2eedd4ed 2222 //see header file for documentation
2223
2224 TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2225 TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2226 TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
48682642 2227
2eedd4ed 2228 Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2229 Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2230
2231 Float_t alfa;
2232 Float_t qt;
2233 if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2234 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2235 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2236 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2237 } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2238 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2239 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2240 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2241 } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2242 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2243 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2244 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2245 } else {
2246 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2247 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2248 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2249 }
2250
3b77b2d1 2251 armenterosQtAlpha[0]=qt;
2252 armenterosQtAlpha[1]=alfa;
2eedd4ed 2253 return 1;
2254
2255}
2256
3b77b2d1 2257Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const TParticle * gammaCandidate, Double_t armenterosQtAlpha[2] ){
2eedd4ed 2258 //see header file for documentation
48682642 2259
2eedd4ed 2260 TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2261 TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2262 TVector3 vecV0(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
2263
2264 Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2265 Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2266
2267 Float_t alfa;
2268 Float_t qt;
2269 if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2270 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2271 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2272 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2273 } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2274 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2275 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2276 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2277 } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2278 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2279 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2280 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2281 } else {
2282 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2283 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2284 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2285 }
2286
3b77b2d1 2287 armenterosQtAlpha[0]=qt;
2288 armenterosQtAlpha[1]=alfa;
2eedd4ed 2289 return 1;
48682642 2290
2291}
2292
2eedd4ed 2293
40051a3e 2294Int_t AliV0Reader::GetFirstTPCRow(Double_t radius){
2295
2296
2eedd4ed 2297 Int_t firstTPCRow=0;
2298 Double_t radiusI = 84.8;
2299 Double_t radiusO = 134.6;
2300 Double_t radiusOB = 198.;
2301 Double_t rSizeI = 0.75;
2302 Double_t rSizeO = 1.;
2303 Double_t rSizeOB = 1.5;
2304 Int_t nClsI=63;
2305 Int_t nClsIO=127;
2306
2307 if(radius <= radiusI){
2308 return firstTPCRow;
2309 }
2310 if(radius>radiusI && radius<=radiusO){
2311 firstTPCRow = (Int_t)((radius-radiusI)/rSizeI);
2312 }
2313 if(radius>radiusO && radius<=radiusOB){
2314 firstTPCRow = (Int_t)(nClsI+(radius-radiusO)/rSizeO);
2315 }
2316
2317 if(radius>radiusOB){
2318 firstTPCRow =(Int_t)(nClsIO+(radius-radiusOB)/rSizeOB);
2319 }
2320
2321
2322 return firstTPCRow;
40051a3e 2323}
2ea0dd4a 2324void AliV0Reader::SmearKFParticle(AliKFParticle * kfParticle)
2325{
2eedd4ed 2326 Double_t facPBrem = 1.;
2327 Double_t facPSig = 0.;
2328
2329 Double_t phi=0.;
2330 Double_t theta=0.;
2331 Double_t P=0.;
2332
2333 P=kfParticle->GetP();
2334 phi=kfParticle->GetPhi();
2335 if( kfParticle->GetP()!=0){
2336 theta=acos( kfParticle->Pz()/ kfParticle->GetP());
2337 }
2338
2339 if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
2340 facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
2341 }
2342
2343 if( fPBremSmearing != 1.){
2344 if(fBrem!=NULL){
2345 facPBrem = fBrem->GetRandom();
2346 }
2347 }
2348
2349 kfParticle->Px() = facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) ;
2350 kfParticle->Py() = facPBrem* (1+facPSig)* P*sin(theta)*sin(phi) ;
2351 kfParticle->Pz() = facPBrem* (1+facPSig)* P*cos(theta) ;
2352 kfParticle->E() = kfParticle->GetP();
2ea0dd4a 2353}
3b77b2d1 2354
2355///________________________________________________________________________
2356const AliExternalTrackParam *AliV0Reader::GetExternalTrackParam(AliESDv0 *v0,Int_t charge){
2357
2358 if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
2359
2360 Int_t label;
2361 if(charge>0)label=0;
2362 else label=1;
2363 // Check for sign flip
2364
2365 if(v0){
2366 if(!v0->GetParamN()||!v0->GetParamP())return 0x0;
2367 if(!GetTrack(v0->GetNindex())||!fESDEvent->GetTrack(v0->GetPindex()))return 0x0;
2368 if((GetTrack(v0->GetPindex()))->Charge()==charge){
2369// fCurrentTrackLabels[label]=v0->GetPindex();
2370 return v0->GetParamP();}
2371 if((GetTrack(v0->GetNindex()))->Charge()==charge){
2372// fCurrentTrackLabels[label]=v0->GetNindex();
2373 return v0->GetParamN();}
2374 }
2375 return 0x0;
2376}
2377
2378///________________________________________________________________________
2379AliVTrack *AliV0Reader::GetTrack(Int_t label){
2380 if(fESDEvent){
2381 return (AliESDtrack*)fESDEvent->GetTrack(label);
2382 }
2383// if(fAODEvent)return (AliAODTrack*)GetAODTrack(label);
2384 return 0x0;
2385}
2386
2387Double_t AliV0Reader::GetV0CosineOfPointingAngle(Double_t V0PointX, Double_t V0PointY, Double_t V0PointZ){
2388 // calculates the pointing angle of the recalculated V0
2389
2390 Double_t momV0[3]; //momentum of the V0
2391 fCurrentV0->GetPxPyPz(momV0[0],momV0[1],momV0[2]);
2392
2393 Double_t PosV0[3]; //Recalculated V0 Position vector
2394
2395 PosV0[0] = V0PointX - fESDEvent->GetPrimaryVertex()->GetX();
2396 PosV0[1] = V0PointY - fESDEvent->GetPrimaryVertex()->GetY();
2397 PosV0[2] = V0PointZ - fESDEvent->GetPrimaryVertex()->GetZ();
2398
2399 Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
2400 Double_t PosV02 = PosV0[0]*PosV0[0] + PosV0[1]*PosV0[1] + PosV0[2]*PosV0[2];
2401
2402 Double_t cosinePointingAngle = (PosV0[0]*momV0[0] + PosV0[1]*momV0[1] + PosV0[2]*momV0[2] ) / TMath::Sqrt(momV02 * PosV02);
2403
2404 return cosinePointingAngle;
2405}
2406
2407
2408Double_t AliV0Reader::GetPsiPair(AliESDv0* v0)
2409{
2410 //
2411 // Angle between daughter momentum plane and plane
2412 //
2413 Float_t magField = fESDEvent->GetMagneticField();
2414
2415 Double_t xyz[3] = {0.,0.,0.};
2416 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
2417
2418 Double_t mn[3] = {0,0,0};
2419 Double_t mp[3] = {0,0,0};
2420
2421 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
2422 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
2423
2424 Double_t deltat = 1.;
2425 deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
2426 Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
2427
2428 Double_t momPosProp[3] = {0,0,0};
2429 Double_t momNegProp[3] = {0,0,0};
2430
2431 AliExternalTrackParam nt = *(v0->GetParamN());
2432 AliExternalTrackParam pt = *(v0->GetParamP());
2433
2434 Double_t psiPair = 4.;
2435 if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2436
2437 if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2438
2439 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
2440 nt.GetPxPyPz(momNegProp);
2441
2442 Double_t pEle =
2443 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
2444
2445 Double_t pPos =
2446 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
2447
2448 Double_t scalarproduct =
2449 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
2450
2451 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
2452
2453 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
2454
2455 return psiPair;
2456}
2457
2458Int_t AliV0Reader::GetNumberOfHijingPlusPythiaPrimeries(){
2459
2460 // Calculate NPrimaries for LHC11a10b_*
2461
2462 Int_t nproduced = 0;
2463 AliGenCocktailEventHeader *cHeader = dynamic_cast<AliGenCocktailEventHeader*>(fMCEvent->GenEventHeader());
2464 if(cHeader){
2465 TList *genHeaders = cHeader->GetHeaders();
2466 AliGenEventHeader* gh = 0;
2467 for(Int_t i = 0; i<genHeaders->GetEntries();i++){
2468 gh = (AliGenEventHeader*)genHeaders->At(i);
2469 TString GeneratorName = gh->GetName();
2470 // cout<<i<<" "<<GeneratorName<<" "<<gh->NProduced()<<endl;
2471 if(GeneratorName.CompareTo("Hijing") == 0){
2472 nproduced = nproduced + gh->NProduced();
2473 }
2474 else if(GeneratorName.CompareTo("Pythia") == 0){
2475 nproduced = nproduced + gh->NProduced();
2476 }
2477 }
2478 }
2479 if(!cHeader){
2480 nproduced = fMCStack->GetNprimary();
2481 }
2482
2483 // cout<<fMCStack->GetNprimary()-nproduced<<endl;
2484
2485 return nproduced;
2486}
2487
2488
2489Bool_t AliV0Reader::IsParticleFromBGEvent(Int_t index){
2490
2491 Bool_t particleFromBG = kFALSE;
2492
2493 if(index == -1) return kFALSE;
2494 if(index > fNumberOfPrimerisFromHijingAndPythia && index < fMCStack->GetNprimary()){
2495 // cout<<index<<" "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2496 return kTRUE;
2497 }
2498 // else cout<<"Passt Noch "<<index<<" "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2499 // cout<<fMCEvent->IsFromBGEvent(index)<<endl;
2500 TParticle *BGParticle = fMCStack->Particle(index);
2501 if(BGParticle->GetMother(0) > -1) return kFALSE;
2502 Int_t indexMother = fMCStack->Particle(index)->GetMother(0);
2503 particleFromBG = IsParticleFromBGEvent(indexMother);
2504
2505 return kFALSE;
2506}