]> 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
sync
[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.),
7033a2a4 171 fExcludeBackgroundEventForGammaCorrection(0),
3b77b2d1 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();
efbc23fc 365 // if(fMCStack == NULL){
366 // //print warning here
367 // }
2eedd4ed 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
92efd725 410// if(gRandom != NULL){
411// delete gRandom;
412// gRandom= new TRandom3(0);
413// }
2eedd4ed 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
92efd725 426 Double_t *zBinLimitsArray = new Double_t[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
92efd725 437 Double_t *multiplicityBinLimitsArray = new Double_t[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){
7033a2a4 489 fNumberOfPrimerisFromHijingAndPythia = GetNumberOfHijingPlusPythiaPrimeries(fExcludeBackgroundEventForGammaCorrection);
3b77b2d1 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
92efd725 1371 const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0);
1372 const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0);
3b77b2d1 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){
92efd725 1390 delete fCurrentNegativeKFParticle;
2eedd4ed 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){
92efd725 1414 fCurrentMotherKFCandidate = new AliKFParticle();//(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
2eedd4ed 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
92efd725 1653 delete fCurrentEventGoodV0s;
1654 fCurrentEventGoodV0s = NULL;
1655 delete fBrem;
1656 fBrem = NULL;
1657 delete fCurrentNegativeKFParticle;
1658 fCurrentNegativeKFParticle = NULL;
1659 delete fCurrentPositiveKFParticle;
1660 fCurrentPositiveKFParticle = NULL;
1661 delete fCurrentMotherKFCandidate;
1662 fCurrentMotherKFCandidate = NULL;
1663 delete fNegativeTrackLorentzVector;
1664 fNegativeTrackLorentzVector = NULL;
1665 delete fPositiveTrackLorentzVector;
1666 fPositiveTrackLorentzVector = NULL;
1667 delete fMotherCandidateLorentzVector;
1668 fMotherCandidateLorentzVector = NULL;
1669
1670 // fBGEventHandler->PrintBGArray(); // for debugging
a0b94e5c 1671}
1672
1673
1674Double_t AliV0Reader::GetNegativeTrackPhi() const{
2eedd4ed 1675 //see header file for documentation
a0b94e5c 1676
2eedd4ed 1677 Double_t offset=0;
1678 if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
1679 offset = -2*TMath::Pi();
1680 }
1681 return fNegativeTrackLorentzVector->Phi()+offset;
a0b94e5c 1682}
1683
1684Double_t AliV0Reader::GetPositiveTrackPhi() const{
2eedd4ed 1685 //see header file for documentation
a0b94e5c 1686
2eedd4ed 1687 Double_t offset=0;
1688 if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
1689 offset = -2*TMath::Pi();
1690 }
1691 return fPositiveTrackLorentzVector->Phi()+offset;
a0b94e5c 1692}
1693
1694Double_t AliV0Reader::GetMotherCandidatePhi() const{
2eedd4ed 1695 //see header file for documentation
a0b94e5c 1696
2eedd4ed 1697 Double_t offset=0;
1698 if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
1699 offset = -2*TMath::Pi();
1700 }
1701 return fMotherCandidateLorentzVector->Phi()+offset;
a0b94e5c 1702}
1703
1704
1705Double_t AliV0Reader::GetMotherCandidateRapidity() const{
2eedd4ed 1706 //see header file for documentation
a0b94e5c 1707
2eedd4ed 1708 Double_t rapidity=0;
1709 if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
1710 else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
1711 return rapidity;
a0b94e5c 1712
1713}
1714
1715
1716
1717
1718
1719Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
2eedd4ed 1720 //see header file for documentation
a0b94e5c 1721
2eedd4ed 1722 Int_t iResult = 10; // Unknown particle
a0b94e5c 1723
2eedd4ed 1724 if(chargeOfTrack==-1){ //negative track
1725 switch(abs(fNegativeTrackPID)){
1726 case 11: //electron
1727 iResult = 0;
1728 break;
1729 case 13: //muon
1730 iResult = 1;
1731 break;
1732 case 211: //pion
1733 iResult = 2;
1734 break;
1735 case 321: //kaon
1736 iResult = 3;
1737 break;
1738 case 2212: //proton
1739 iResult = 4;
1740 break;
1741 case 22: //photon
1742 iResult = 5;
1743 break;
1744 case 111: //pi0
1745 iResult = 6;
1746 break;
1747 case 2112: //neutron
1748 iResult = 7;
1749 break;
1750 case 311: //K0
1751 iResult = 8;
1752 break;
a0b94e5c 1753
2eedd4ed 1754 //Put in here for kSPECIES::kEleCon ????
1755 }
1756 }
1757 else if(chargeOfTrack==1){ //positive track
1758 switch(abs(fPositiveTrackPID)){
1759 case 11: //electron
1760 iResult = 0;
1761 break;
1762 case 13: //muon
1763 iResult = 1;
1764 break;
1765 case 211: //pion
1766 iResult = 2;
1767 break;
1768 case 321: //kaon
1769 iResult = 3;
1770 break;
1771 case 2212: //proton
1772 iResult = 4;
1773 break;
1774 case 22: //photon
1775 iResult = 5;
1776 break;
1777 case 111: //pi0
1778 iResult = 6;
1779 break;
1780 case 2112: //neutron
1781 iResult = 7;
1782 break;
1783 case 311: //K0
1784 iResult = 8;
1785 break;
a0b94e5c 1786
2eedd4ed 1787 //Put in here for kSPECIES::kEleCon ????
1788 }
1789 }
1790 else{
1791 //Wrong parameter.. Print warning
1792 }
1793 return iResult;
a0b94e5c 1794}
1795
3b77b2d1 1796//_____________________________________________________________________________________________________________________
1797Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
2eedd4ed 1798 // see header file for documentation
2eedd4ed 1799
1800 Double_t helix[6];
1801 track->GetHelixParameters(helix,b);
1802
1803 Double_t xpos = helix[5];
1804 Double_t ypos = helix[0];
1805 Double_t radius = TMath::Abs(1./helix[4]);
1806 Double_t phi = helix[2];
1807
1808 if(phi < 0){
3b77b2d1 1809 phi = phi + 2*TMath::Pi();
2eedd4ed 1810 }
1811
3b77b2d1 1812 phi -= TMath::Pi()/2.;
2eedd4ed 1813 Double_t xpoint = radius * TMath::Cos(phi);
1814 Double_t ypoint = radius * TMath::Sin(phi);
1815
1816 if(b<0){
3b77b2d1 1817 if(charge > 0){
1818 xpoint = - xpoint;
1819 ypoint = - ypoint;
1820 }
1821
1822 if(charge < 0){
1823 xpoint = xpoint;
1824 ypoint = ypoint;
1825 }
2eedd4ed 1826 }
3b77b2d1 1827 if(b>0){
1828 if(charge > 0){
1829 xpoint = xpoint;
1830 ypoint = ypoint;
1831 }
2eedd4ed 1832
3b77b2d1 1833 if(charge < 0){
1834 xpoint = - xpoint;
1835 ypoint = - ypoint;
1836 }
1837 }
2eedd4ed 1838 center[0] = xpos + xpoint;
1839 center[1] = ypos + ypoint;
1840
1841 return 1;
a0b94e5c 1842}
1843
3b77b2d1 1844//_________________________________________________________________________________________________________
1845// Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
1846// // Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1847// // see header file for documentation
1848//
1849// Double_t pi = 3.14159265358979323846;
1850//
1851// Double_t helix[6];
1852// track->GetHelixParameters(helix,b);
1853//
1854// Double_t xpos = helix[5];
1855// Double_t ypos = helix[0];
1856// Double_t radius = TMath::Abs(1./helix[4]);
1857// Double_t phi = helix[2];
1858//
1859// if(phi < 0){
1860// phi = phi + 2*pi;
1861// }
1862//
1863// phi -= pi/2.;
1864// Double_t xpoint = radius * TMath::Cos(phi);
1865// Double_t ypoint = radius * TMath::Sin(phi);
1866//
1867// if(b<0){
1868// if(charge > 0){
1869// xpoint = - xpoint;
1870// ypoint = - ypoint;
1871// }
1872//
1873// if(charge < 0){
1874// xpoint = xpoint;
1875// ypoint = ypoint;
1876// }
1877// }
1878// if(b>0){
1879// if(charge > 0){
1880// xpoint = xpoint;
1881// ypoint = ypoint;
1882// }
1883//
1884// if(charge < 0){
1885// xpoint = - xpoint;
1886// ypoint = - ypoint;
1887// }
1888// }
1889// center[0] = xpos + xpoint;
1890// center[1] = ypos + ypoint;
1891//
1892// return 1;
1893// }
1894
1895//____________________________________________________________________________________________________________________________
1896Bool_t AliV0Reader::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3]){
1897
1898 if(!pparam||!nparam)return kFALSE;
2eedd4ed 1899
1900 Double_t helixcenterpos[2];
1901 GetHelixCenter(pparam,GetMagneticField(),pparam->Charge(),helixcenterpos);
3b77b2d1 1902
2eedd4ed 1903 Double_t helixcenterneg[2];
1904 GetHelixCenter(nparam,GetMagneticField(),nparam->Charge(),helixcenterneg);
3b77b2d1 1905
2eedd4ed 1906 Double_t helixpos[6];
1907 pparam->GetHelixParameters(helixpos,GetMagneticField());
1908 Double_t posradius = TMath::Abs(1./helixpos[4]);
3b77b2d1 1909
2eedd4ed 1910 Double_t helixneg[6];
1911 nparam->GetHelixParameters(helixneg,GetMagneticField());
1912 Double_t negradius = TMath::Abs(1./helixneg[4]);
3b77b2d1 1913
1914 // Calculate xy-position
1915
2eedd4ed 1916 Double_t xpos = helixcenterpos[0];
1917 Double_t ypos = helixcenterpos[1];
1918 Double_t xneg = helixcenterneg[0];
1919 Double_t yneg = helixcenterneg[1];
3b77b2d1 1920
2eedd4ed 1921 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
3b77b2d1 1922 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
1923
1924
2eedd4ed 1925 // Calculate z-position
2eedd4ed 1926
3b77b2d1 1927 Double_t deltaXPos = convpos[0] - xpos;
1928 Double_t deltaYPos = convpos[1] - ypos;
2eedd4ed 1929
3b77b2d1 1930 Double_t deltaXNeg = convpos[0] - xneg;
1931 Double_t deltaYNeg = convpos[1] - yneg;
2eedd4ed 1932
3b77b2d1 1933 Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
1934 Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
2eedd4ed 1935
3b77b2d1 1936 Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
1937 Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
1938
1939 Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
1940 Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
1941
1942 Double_t b = fESDEvent->GetMagneticField();
1943
1944 AliExternalTrackParam p(*pparam);
1945 AliExternalTrackParam n(*nparam);
1946
1947 TVector2 vertexPos(vertexXPos,vertexYPos);
1948 TVector2 vertexNeg(vertexXNeg,vertexYNeg);
1949
1950 // Convert to local coordinate system
1951 vertexPos=vertexPos.Rotate(-p.GetAlpha());
1952 vertexNeg=vertexNeg.Rotate(-p.GetAlpha());
1953
1954 // Propagate Track Params to Vertex
1955 p.PropagateTo(vertexPos.X(),b);
1956 n.PropagateTo(vertexNeg.X(),b);
2eedd4ed 1957
3b77b2d1 1958 convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
2eedd4ed 1959
3b77b2d1 1960 return kTRUE;
1961}
1962//
1963// //__________________________________________________________________________________________________________
4a6157dc 1964Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
2eedd4ed 1965 //see header file for documentation
a0b94e5c 1966
2eedd4ed 1967 Double_t helixcenterpos[2];
1968 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
a0b94e5c 1969
2eedd4ed 1970 Double_t helixcenterneg[2];
1971 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
a0b94e5c 1972
2eedd4ed 1973 Double_t poshelix[6];
1974 ptrack->GetHelixParameters(poshelix,b);
1975 Double_t posradius = TMath::Abs(1./poshelix[4]);
a0b94e5c 1976
2eedd4ed 1977 Double_t neghelix[6];
1978 ntrack->GetHelixParameters(neghelix,b);
1979 Double_t negradius = TMath::Abs(1./neghelix[4]);
a0b94e5c 1980
2eedd4ed 1981 Double_t xpos = helixcenterpos[0];
1982 Double_t ypos = helixcenterpos[1];
1983 Double_t xneg = helixcenterneg[0];
1984 Double_t yneg = helixcenterneg[1];
a0b94e5c 1985
2eedd4ed 1986 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1987 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
a0b94e5c 1988
2eedd4ed 1989 return 1;
a0b94e5c 1990}
3b77b2d1 1991//
1992//
1993//
4a6157dc 1994Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
2eedd4ed 1995 //see header file for documentation
a0b94e5c 1996
2eedd4ed 1997 Double_t helixpos[6];
1998 ptrack->GetHelixParameters(helixpos,b);
a0b94e5c 1999
2eedd4ed 2000 Double_t helixneg[6];
2001 ntrack->GetHelixParameters(helixneg,b);
a0b94e5c 2002
2eedd4ed 2003 Double_t negtrackradius = TMath::Abs(1./helixneg[4]);
2004 Double_t postrackradius = TMath::Abs(1./helixpos[4]);
a0b94e5c 2005
2eedd4ed 2006 Double_t pi = 3.14159265358979323846;
a0b94e5c 2007
2eedd4ed 2008 Double_t convpos[2];
2009 GetConvPosXY(ptrack,ntrack,b,convpos);
a0b94e5c 2010
2eedd4ed 2011 Double_t convposx = convpos[0];
2012 Double_t convposy = convpos[1];
a0b94e5c 2013
2eedd4ed 2014 Double_t helixcenterpos[2];
2015 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
a0b94e5c 2016
2eedd4ed 2017 Double_t helixcenterneg[2];
2018 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
a0b94e5c 2019
2eedd4ed 2020 Double_t xpos = helixcenterpos[0];
2021 Double_t ypos = helixcenterpos[1];
2022 Double_t xneg = helixcenterneg[0];
2023 Double_t yneg = helixcenterneg[1];
a0b94e5c 2024
2eedd4ed 2025 Double_t deltaXPos = convposx - xpos;
2026 Double_t deltaYPos = convposy - ypos;
a0b94e5c 2027
2eedd4ed 2028 Double_t deltaXNeg = convposx - xneg;
2029 Double_t deltaYNeg = convposy - yneg;
a0b94e5c 2030
2eedd4ed 2031 Double_t alphaPos = pi + TMath::ATan2(-deltaYPos,-deltaXPos);
2032 Double_t alphaNeg = pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
a0b94e5c 2033
2eedd4ed 2034 Double_t vertexXNeg = xneg + TMath::Abs(negtrackradius)*
2035 TMath::Cos(alphaNeg);
2036 Double_t vertexYNeg = yneg + TMath::Abs(negtrackradius)*
2037 TMath::Sin(alphaNeg);
a0b94e5c 2038
2eedd4ed 2039 Double_t vertexXPos = xpos + TMath::Abs(postrackradius)*
2040 TMath::Cos(alphaPos);
2041 Double_t vertexYPos = ypos + TMath::Abs(postrackradius)*
2042 TMath::Sin(alphaPos);
a0b94e5c 2043
2eedd4ed 2044 Double_t x0neg = helixneg[5];
2045 Double_t y0neg = helixneg[0];
a0b94e5c 2046
2eedd4ed 2047 Double_t x0pos = helixpos[5];
2048 Double_t y0pos = helixpos[0];
a0b94e5c 2049
2eedd4ed 2050 Double_t dNeg = TMath::Sqrt((vertexXNeg - x0neg)*(vertexXNeg - x0neg)
2051 +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
a0b94e5c 2052
2eedd4ed 2053 Double_t dPos = TMath::Sqrt((vertexXPos - x0pos)*(vertexXPos - x0pos)
2054 +(vertexYPos - y0pos)*(vertexYPos - y0pos));
a0b94e5c 2055
2eedd4ed 2056 Double_t rNeg = TMath::Sqrt(negtrackradius*negtrackradius -
2057 dNeg*dNeg/4.);
a0b94e5c 2058
2eedd4ed 2059 Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
2060 dPos*dPos/4.);
a0b94e5c 2061
2eedd4ed 2062 Double_t deltabetaNeg = 2*(pi + TMath::ATan2(-dNeg/2.,-rNeg));
2063 Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
a0b94e5c 2064
2eedd4ed 2065 Double_t deltaUNeg = negtrackradius*deltabetaNeg;
2066 Double_t deltaUPos = postrackradius*deltabetaPos;
a0b94e5c 2067
2eedd4ed 2068 Double_t zphaseNeg = ntrack->GetZ() + deltaUNeg * ntrack->GetTgl();
2069 Double_t zphasePos = ptrack->GetZ() + deltaUPos * ptrack->GetTgl();
a0b94e5c 2070
2eedd4ed 2071 Double_t convposz = (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
a0b94e5c 2072
2eedd4ed 2073 return convposz;
a0b94e5c 2074}
3b77b2d1 2075//
f14266e7 2076AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t /*event*/) const{
2eedd4ed 2077 /*
2078 if(fUseChargedTrackMultiplicityForBG == kTRUE){
2079 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
2080 }
2081 else{ // means we use #v0s as multiplicity
2082 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
2083 }
2084 */
2085 return NULL;
037dc2db 2086}
2087
2088Int_t AliV0Reader::CountESDTracks(){
2eedd4ed 2089 // see header file for documentation
2090 if(fNumberOfESDTracks == 0){ // count the good esd tracks
2091 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2092 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2093 if(!curTrack){
037dc2db 2094 continue;
2eedd4ed 2095 }
2096 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
037dc2db 2097 fNumberOfESDTracks++;
2eedd4ed 2098 }
2099 }
2100 }
037dc2db 2101
2eedd4ed 2102 return fNumberOfESDTracks;
037dc2db 2103}
2104
2105Bool_t AliV0Reader::CheckIfPi0IsMother(Int_t label){
2eedd4ed 2106 // see headerfile for documentation
2107 Bool_t iResult=kFALSE;
2108 // cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2109 if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 111){
2110 iResult=kTRUE;
2111 }
2112 return iResult;
5e55d806 2113}
48682642 2114
10e3319b 2115Bool_t AliV0Reader::CheckIfEtaIsMother(Int_t label){
2eedd4ed 2116 // see headerfile for documentation
2117 Bool_t iResult=kFALSE;
2118 // cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2119 if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 221){
2120 iResult=kTRUE;
2121 }
2122 return iResult;
10e3319b 2123}
2124
48682642 2125
3b77b2d1 2126
2127
2128Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, Double_t armenterosQtAlpha[2], Int_t kfProductionMethod){
2129
2130 AliKFParticle PosParticle = *positiveKFParticle;
2131 AliKFParticle NegParticle = *negativeKFParticle;
2132 AliKFParticle Gamma;
2133 if(kfProductionMethod < 3)
2134 Gamma.ConstructGamma(PosParticle, NegParticle);
2135 else if(kfProductionMethod == 3){
2136 Gamma += PosParticle;
2137 Gamma += NegParticle;
2138 }
2139
2140 Double_t VertexGamma[3] = {Gamma.GetX(), Gamma.GetY(), Gamma.GetZ()};
2141 PosParticle.TransportToPoint(VertexGamma);
2142 NegParticle.TransportToPoint(VertexGamma);
2143
2144 AliKFParticle::GetArmenterosPodolanski(PosParticle, NegParticle, armenterosQtAlpha);
2145
2146 return 1;
2147}
2148
2149
2150
2151
2152Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2eedd4ed 2153 //see header file for documentation
48682642 2154
2eedd4ed 2155 TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
2156 TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
2157 TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
48682642 2158
2eedd4ed 2159 Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2160 Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2161
2162 Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2163 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2164
2165
2166 Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2167
3b77b2d1 2168 armenterosQtAlpha[0]=qt;
2169 armenterosQtAlpha[1]=alfa;
2eedd4ed 2170
2171 return 1;
2172
2173}
2174
3b77b2d1 2175Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliESDv0* v0, Double_t armenterosQtAlpha[2])
2176{ //see header file for documentation
2177
2178 Double_t mn[3] = {0,0,0};
2179 Double_t mp[3] = {0,0,0};
2180 Double_t mm[3] = {0,0,0};
2181
2182
2183 Int_t pIndex = 0, nIndex = 0;
2184 pIndex = v0->GetPindex();
2185 nIndex = v0->GetNindex();
2186
2187 AliESDtrack* d[2];
2188 d[0] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(pIndex));
2189 d[1] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(nIndex));
2190
2191 Int_t sign[2];
2192 sign[0] = (int)d[0]->GetSign();
2193 sign[1] = (int)d[1]->GetSign();
2194
2195 Bool_t correct = kFALSE;
2196
2197
2198 if(-1 == sign[0] && 1 == sign[1]){
2199 correct = kFALSE;
2200 }
2201 else{
2202 correct = kTRUE;
2203 }
2204
2205
2206 if(correct){
2207 v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2208 v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2209 }
2210 else{
2211 v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2212 v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2213 }
2214 v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
2215
2216 TVector3 vecN(mn[0],mn[1],mn[2]);
2217 TVector3 vecP(mp[0],mp[1],mp[2]);
2218 TVector3 vecM(mm[0],mm[1],mm[2]);
2219
2220 Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
2221 Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
2222
2223 Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
2224 ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
2225 Double_t qt = vecP.Mag()*sin(thetaP);
2226
2227 armenterosQtAlpha[0]=qt;
2228 armenterosQtAlpha[1]=alfa;
2229
2230 return 1;
2231
2232}
2233
2234
2235
2236Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2eedd4ed 2237 //see header file for documentation
2238
2239 TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2240 TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2241 TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
48682642 2242
2eedd4ed 2243 Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2244 Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2245
2246 Float_t alfa;
2247 Float_t qt;
2248 if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2249 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2250 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2251 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2252 } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2253 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2254 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2255 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2256 } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2257 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2258 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2259 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2260 } else {
2261 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2262 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2263 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2264 }
2265
3b77b2d1 2266 armenterosQtAlpha[0]=qt;
2267 armenterosQtAlpha[1]=alfa;
2eedd4ed 2268 return 1;
2269
2270}
2271
3b77b2d1 2272Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const TParticle * gammaCandidate, Double_t armenterosQtAlpha[2] ){
2eedd4ed 2273 //see header file for documentation
48682642 2274
2eedd4ed 2275 TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2276 TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2277 TVector3 vecV0(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
2278
2279 Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2280 Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2281
2282 Float_t alfa;
2283 Float_t qt;
2284 if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2285 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2286 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2287 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2288 } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2289 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2290 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2291 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2292 } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2293 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2294 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2295 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2296 } else {
2297 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2298 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2299 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2300 }
2301
3b77b2d1 2302 armenterosQtAlpha[0]=qt;
2303 armenterosQtAlpha[1]=alfa;
2eedd4ed 2304 return 1;
48682642 2305
2306}
2307
2eedd4ed 2308
40051a3e 2309Int_t AliV0Reader::GetFirstTPCRow(Double_t radius){
2310
2311
2eedd4ed 2312 Int_t firstTPCRow=0;
2313 Double_t radiusI = 84.8;
2314 Double_t radiusO = 134.6;
2315 Double_t radiusOB = 198.;
2316 Double_t rSizeI = 0.75;
2317 Double_t rSizeO = 1.;
2318 Double_t rSizeOB = 1.5;
2319 Int_t nClsI=63;
2320 Int_t nClsIO=127;
2321
2322 if(radius <= radiusI){
2323 return firstTPCRow;
2324 }
2325 if(radius>radiusI && radius<=radiusO){
2326 firstTPCRow = (Int_t)((radius-radiusI)/rSizeI);
2327 }
2328 if(radius>radiusO && radius<=radiusOB){
2329 firstTPCRow = (Int_t)(nClsI+(radius-radiusO)/rSizeO);
2330 }
2331
2332 if(radius>radiusOB){
2333 firstTPCRow =(Int_t)(nClsIO+(radius-radiusOB)/rSizeOB);
2334 }
2335
2336
2337 return firstTPCRow;
40051a3e 2338}
2ea0dd4a 2339void AliV0Reader::SmearKFParticle(AliKFParticle * kfParticle)
2340{
2eedd4ed 2341 Double_t facPBrem = 1.;
2342 Double_t facPSig = 0.;
2343
2344 Double_t phi=0.;
2345 Double_t theta=0.;
2346 Double_t P=0.;
2347
2348 P=kfParticle->GetP();
2349 phi=kfParticle->GetPhi();
2350 if( kfParticle->GetP()!=0){
2351 theta=acos( kfParticle->Pz()/ kfParticle->GetP());
2352 }
2353
2354 if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
2355 facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
2356 }
2357
2358 if( fPBremSmearing != 1.){
2359 if(fBrem!=NULL){
2360 facPBrem = fBrem->GetRandom();
2361 }
2362 }
2363
2364 kfParticle->Px() = facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) ;
2365 kfParticle->Py() = facPBrem* (1+facPSig)* P*sin(theta)*sin(phi) ;
2366 kfParticle->Pz() = facPBrem* (1+facPSig)* P*cos(theta) ;
2367 kfParticle->E() = kfParticle->GetP();
2ea0dd4a 2368}
3b77b2d1 2369
2370///________________________________________________________________________
2371const AliExternalTrackParam *AliV0Reader::GetExternalTrackParam(AliESDv0 *v0,Int_t charge){
2372
2373 if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
2374
2375 Int_t label;
2376 if(charge>0)label=0;
2377 else label=1;
2378 // Check for sign flip
2379
2380 if(v0){
2381 if(!v0->GetParamN()||!v0->GetParamP())return 0x0;
2382 if(!GetTrack(v0->GetNindex())||!fESDEvent->GetTrack(v0->GetPindex()))return 0x0;
2383 if((GetTrack(v0->GetPindex()))->Charge()==charge){
2384// fCurrentTrackLabels[label]=v0->GetPindex();
2385 return v0->GetParamP();}
2386 if((GetTrack(v0->GetNindex()))->Charge()==charge){
2387// fCurrentTrackLabels[label]=v0->GetNindex();
2388 return v0->GetParamN();}
2389 }
2390 return 0x0;
2391}
2392
2393///________________________________________________________________________
2394AliVTrack *AliV0Reader::GetTrack(Int_t label){
2395 if(fESDEvent){
2396 return (AliESDtrack*)fESDEvent->GetTrack(label);
2397 }
2398// if(fAODEvent)return (AliAODTrack*)GetAODTrack(label);
2399 return 0x0;
2400}
2401
2402Double_t AliV0Reader::GetV0CosineOfPointingAngle(Double_t V0PointX, Double_t V0PointY, Double_t V0PointZ){
2403 // calculates the pointing angle of the recalculated V0
2404
2405 Double_t momV0[3]; //momentum of the V0
2406 fCurrentV0->GetPxPyPz(momV0[0],momV0[1],momV0[2]);
2407
2408 Double_t PosV0[3]; //Recalculated V0 Position vector
2409
2410 PosV0[0] = V0PointX - fESDEvent->GetPrimaryVertex()->GetX();
2411 PosV0[1] = V0PointY - fESDEvent->GetPrimaryVertex()->GetY();
2412 PosV0[2] = V0PointZ - fESDEvent->GetPrimaryVertex()->GetZ();
2413
2414 Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
2415 Double_t PosV02 = PosV0[0]*PosV0[0] + PosV0[1]*PosV0[1] + PosV0[2]*PosV0[2];
2416
2417 Double_t cosinePointingAngle = (PosV0[0]*momV0[0] + PosV0[1]*momV0[1] + PosV0[2]*momV0[2] ) / TMath::Sqrt(momV02 * PosV02);
2418
2419 return cosinePointingAngle;
2420}
2421
2422
2423Double_t AliV0Reader::GetPsiPair(AliESDv0* v0)
2424{
2425 //
2426 // Angle between daughter momentum plane and plane
2427 //
2428 Float_t magField = fESDEvent->GetMagneticField();
2429
2430 Double_t xyz[3] = {0.,0.,0.};
2431 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
2432
2433 Double_t mn[3] = {0,0,0};
2434 Double_t mp[3] = {0,0,0};
2435
2436 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
2437 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
2438
2439 Double_t deltat = 1.;
2440 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
2441 Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
2442
2443 Double_t momPosProp[3] = {0,0,0};
2444 Double_t momNegProp[3] = {0,0,0};
2445
2446 AliExternalTrackParam nt = *(v0->GetParamN());
2447 AliExternalTrackParam pt = *(v0->GetParamP());
2448
2449 Double_t psiPair = 4.;
2450 if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2451
2452 if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2453
2454 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
2455 nt.GetPxPyPz(momNegProp);
2456
2457 Double_t pEle =
2458 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
2459
2460 Double_t pPos =
2461 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
2462
2463 Double_t scalarproduct =
2464 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
2465
2466 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
2467
2468 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
2469
2470 return psiPair;
2471}
2472
7033a2a4 2473Int_t AliV0Reader::GetNumberOfHijingPlusPythiaPrimeries(Int_t excludeHeaderType){
3b77b2d1 2474
2475 // Calculate NPrimaries for LHC11a10b_*
2476
2477 Int_t nproduced = 0;
2478 AliGenCocktailEventHeader *cHeader = dynamic_cast<AliGenCocktailEventHeader*>(fMCEvent->GenEventHeader());
2479 if(cHeader){
2480 TList *genHeaders = cHeader->GetHeaders();
2481 AliGenEventHeader* gh = 0;
2482 for(Int_t i = 0; i<genHeaders->GetEntries();i++){
2483 gh = (AliGenEventHeader*)genHeaders->At(i);
2484 TString GeneratorName = gh->GetName();
7033a2a4 2485
3b77b2d1 2486 if(GeneratorName.CompareTo("Hijing") == 0){
2487 nproduced = nproduced + gh->NProduced();
7033a2a4 2488 // cout<<i<<" "<<GeneratorName<<" "<<gh->NProduced()<<endl;
3b77b2d1 2489 }
7033a2a4 2490 else if(GeneratorName.CompareTo("Pythia") == 0 && excludeHeaderType == 1){
3b77b2d1 2491 nproduced = nproduced + gh->NProduced();
7033a2a4 2492 // cout<<i<<" "<<GeneratorName<<" "<<gh->NProduced()<<endl;
3b77b2d1 2493 }
2494 }
2495 }
2496 if(!cHeader){
2497 nproduced = fMCStack->GetNprimary();
2498 }
2499
2500 // cout<<fMCStack->GetNprimary()-nproduced<<endl;
2501
2502 return nproduced;
2503}
2504
2505
2506Bool_t AliV0Reader::IsParticleFromBGEvent(Int_t index){
2507
2508 Bool_t particleFromBG = kFALSE;
2509
2510 if(index == -1) return kFALSE;
2511 if(index > fNumberOfPrimerisFromHijingAndPythia && index < fMCStack->GetNprimary()){
7033a2a4 2512 // cout<<fMCEvent->IsFromBGEvent(index)<<endl;
3b77b2d1 2513 // cout<<index<<" "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2514 return kTRUE;
2515 }
2516 // else cout<<"Passt Noch "<<index<<" "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2517 // cout<<fMCEvent->IsFromBGEvent(index)<<endl;
2518 TParticle *BGParticle = fMCStack->Particle(index);
7033a2a4 2519 if(BGParticle->GetMother(0) > -1) return kFALSE;
3b77b2d1 2520 Int_t indexMother = fMCStack->Particle(index)->GetMother(0);
2521 particleFromBG = IsParticleFromBGEvent(indexMother);
2522
2523 return kFALSE;
2524}