]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/GammaConv/AliV0Reader.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix, some mc histograms where filled outside mc flag
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliV0Reader.cxx
CommitLineData
a0b94e5c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
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. *
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"
29#include "AliESDtrack.h"
30#include "AliMCEvent.h"
31#include "AliKFVertex.h"
32
33#include "AliStack.h"
34#include "AliMCEventHandler.h"
9640a3d1 35#include "AliTPCpidESD.h"
a0b94e5c 36
37class iostream;
38class AliESDv0;
39class TFormula;
40
41using namespace std;
42
43ClassImp(AliV0Reader)
44
45
46
47AliV0Reader::AliV0Reader() :
48 TObject(),
49 fMCStack(NULL),
50 fMCTruth(NULL),
51 fMCEvent(NULL), // for CF
52 fChain(NULL),
53 fESDHandler(NULL),
54 fESDEvent(NULL),
4a6157dc 55 fCFManager(NULL),
9640a3d1 56 fTPCpid(NULL),
a0b94e5c 57 fHistograms(NULL),
58 fCurrentV0IndexNumber(0),
59 fCurrentV0(NULL),
60 fCurrentNegativeKFParticle(NULL),
61 fCurrentPositiveKFParticle(NULL),
62 fCurrentMotherKFCandidate(NULL),
63 fCurrentNegativeESDTrack(NULL),
64 fCurrentPositiveESDTrack(NULL),
65 fNegativeTrackLorentzVector(NULL),
66 fPositiveTrackLorentzVector(NULL),
67 fMotherCandidateLorentzVector(NULL),
68 fCurrentXValue(0),
69 fCurrentYValue(0),
70 fCurrentZValue(0),
71 fPositiveTrackPID(0),
72 fNegativeTrackPID(0),
73 fNegativeMCParticle(NULL),
74 fPositiveMCParticle(NULL),
75 fMotherMCParticle(NULL),
76 fMotherCandidateKFMass(0),
77 fMotherCandidateKFWidth(0),
78 fUseKFParticle(kTRUE),
79 fUseESDTrack(kFALSE),
80 fDoMC(kFALSE),
81 fMaxR(10000),// 100 meter(outside of ALICE)
82 fEtaCut(0.),
83 fPtCut(0.),
84 fMaxZ(0.),
85 fLineCutZRSlope(0.),
86 fLineCutZValue(0.),
87 fChi2CutConversion(0.),
88 fChi2CutMeson(0.),
89 fPIDProbabilityCutNegativeParticle(0),
90 fPIDProbabilityCutPositiveParticle(0),
9640a3d1 91 fDodEdxSigmaCut(kFALSE),
92 fPIDnSigmaAboveElectronLine(100),
93 fPIDnSigmaBelowElectronLine(-100),
94 fPIDnSigmaAbovePionLine(-100),
95 fPIDMinPnSigmaAbovePionLine(100),
a0b94e5c 96 fXVertexCut(0.),
97 fYVertexCut(0.),
98 fZVertexCut(0.),
99 fNSigmaMass(0.),
100 fUseImprovedVertex(kFALSE),
101 fUseOwnXYZCalculation(kFALSE),
1e7846f4 102 fDoCF(kFALSE),
77880bd8 103 fUseOnFlyV0Finder(kTRUE),
cb90a330 104 fUpdateV0AlreadyCalled(kFALSE),
a0b94e5c 105 fCurrentEventGoodV0s(),
106 fPreviousEventGoodV0s()
107{
9640a3d1 108 fTPCpid = new AliTPCpidESD;
a0b94e5c 109}
110
111
112AliV0Reader::AliV0Reader(const AliV0Reader & original) :
113 TObject(original),
114 fMCStack(original.fMCStack),
115 fMCTruth(original.fMCTruth),
116 fMCEvent(original.fMCEvent), // for CF
117 fChain(original.fChain),
118 fESDHandler(original.fESDHandler),
119 fESDEvent(original.fESDEvent),
4a6157dc 120 fCFManager(original.fCFManager),
9640a3d1 121 fTPCpid(original.fTPCpid),
a0b94e5c 122 fHistograms(original.fHistograms),
123 fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),
124 fCurrentV0(original.fCurrentV0),
125 fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),
126 fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),
127 fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),
128 fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),
129 fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),
130 fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),
131 fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),
132 fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),
133 fCurrentXValue(original.fCurrentXValue),
134 fCurrentYValue(original.fCurrentYValue),
135 fCurrentZValue(original.fCurrentZValue),
136 fPositiveTrackPID(original.fPositiveTrackPID),
137 fNegativeTrackPID(original.fNegativeTrackPID),
138 fNegativeMCParticle(original.fNegativeMCParticle),
139 fPositiveMCParticle(original.fPositiveMCParticle),
140 fMotherMCParticle(original.fMotherMCParticle),
141 fMotherCandidateKFMass(original.fMotherCandidateKFMass),
142 fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),
143 fUseKFParticle(kTRUE),
144 fUseESDTrack(kFALSE),
145 fDoMC(kFALSE),
146 fMaxR(original.fMaxR),
147 fEtaCut(original.fEtaCut),
148 fPtCut(original.fPtCut),
149 fMaxZ(original.fMaxZ),
150 fLineCutZRSlope(original.fLineCutZRSlope),
151 fLineCutZValue(original.fLineCutZValue),
152 fChi2CutConversion(original.fChi2CutConversion),
153 fChi2CutMeson(original.fChi2CutMeson),
154 fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),
155 fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),
9640a3d1 156 fDodEdxSigmaCut(original.fDodEdxSigmaCut),
157 fPIDnSigmaAboveElectronLine(original.fPIDnSigmaAboveElectronLine),
158 fPIDnSigmaBelowElectronLine(original.fPIDnSigmaBelowElectronLine),
159 fPIDnSigmaAbovePionLine(original.fPIDnSigmaAbovePionLine),
160 fPIDMinPnSigmaAbovePionLine(original.fPIDMinPnSigmaAbovePionLine),
a0b94e5c 161 fXVertexCut(original.fXVertexCut),
162 fYVertexCut(original.fYVertexCut),
163 fZVertexCut(original.fZVertexCut),
164 fNSigmaMass(original.fNSigmaMass),
165 fUseImprovedVertex(original.fUseImprovedVertex),
166 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
1e7846f4 167 fDoCF(original.fDoCF),
77880bd8 168 fUseOnFlyV0Finder(original.fUseOnFlyV0Finder),
cb90a330 169 fUpdateV0AlreadyCalled(original.fUpdateV0AlreadyCalled),
a0b94e5c 170 fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
171 fPreviousEventGoodV0s(original.fPreviousEventGoodV0s)
172{
173
174}
175
176
177AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)
178{
179 // assignment operator
180 return *this;
181}
9640a3d1 182AliV0Reader::~AliV0Reader()
183{
184 if(fTPCpid){
185 delete fTPCpid;
186 }
187}
a0b94e5c 188
189void AliV0Reader::Initialize(){
190 //see header file for documentation
cb90a330 191 fUpdateV0AlreadyCalled = kFALSE;
a0b94e5c 192 // Get the input handler from the manager
193 fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
194 if(fESDHandler == NULL){
195 //print warning here
196 }
197
198 // Get pointer to esd event from input handler
199 fESDEvent = fESDHandler->GetEvent();
200 if(fESDEvent == NULL){
201 //print warning here
202 }
203
204 //Get pointer to MCTruth
205 fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
206 if(fMCTruth == NULL){
207 //print warning here
61374d97 208 fDoMC = kFALSE;
a0b94e5c 209 }
210
211 //Get pointer to the mc stack
61374d97 212 if(fMCTruth){
213 fMCStack = fMCTruth->MCEvent()->Stack();
214 if(fMCStack == NULL){
215 //print warning here
216 }
a0b94e5c 217 }
218
a0b94e5c 219 // for CF
220 //Get pointer to the mc event
1e7846f4 221 if(fDoCF){
222 fMCEvent = fMCTruth->MCEvent();
223 if(fMCEvent == NULL){
224 //print warning here
225 fDoCF = kFALSE;
226 }
227 }
a0b94e5c 228
229 AliKFParticle::SetField(fESDEvent->GetMagneticField());
230
231}
232
233AliESDv0* AliV0Reader::GetV0(Int_t index){
234 //see header file for documentation
235 fCurrentV0 = fESDEvent->GetV0(index);
236 UpdateV0Information();
237 return fCurrentV0;
238}
239
240Bool_t AliV0Reader::CheckForPrimaryVertex(){
241 return fESDEvent->GetPrimaryVertex()->GetNContributors()>0;
242}
243
77880bd8 244
245Bool_t AliV0Reader::CheckV0FinderStatus(Int_t index){
246
247 if(fUseOnFlyV0Finder){
77880bd8 248 if(!GetV0(index)->GetOnFlyStatus()){
249 return kFALSE;
250 }
251 }
252 if(!fUseOnFlyV0Finder){
77880bd8 253 if(!GetV0(index)->GetOnFlyStatus()){
254 return kFALSE;
255 }
256 }
77880bd8 257 return kTRUE;
258}
259
260
261
a0b94e5c 262Bool_t AliV0Reader::NextV0(){
263 //see header file for documentation
264
265 Bool_t iResult=kFALSE;
266 while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){
267 fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);
268
269 // moved it up here so that the correction framework can access pt and eta information
270 if(UpdateV0Information() == kFALSE){
271 fCurrentV0IndexNumber++;
272 continue;
273 }
cb90a330 274
a0b94e5c 275 Double_t containerInput[3];
1e7846f4 276 if(fDoCF){
277 containerInput[0] = GetMotherCandidatePt();
278 containerInput[1] = GetMotherCandidateEta();
279 containerInput[2] = GetMotherCandidateMass();
280 }
cb90a330 281
282
a0b94e5c 283 //checks if on the fly mode is set
284 if ( !fCurrentV0->GetOnFlyStatus() ){
285 if(fHistograms != NULL){
286 fHistograms->FillHistogram("ESD_CutGetOnFly_InvMass",GetMotherCandidateMass());
287 }
288 fCurrentV0IndexNumber++;
289 continue;
290 }
1e7846f4 291 if(fDoCF){
292 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly); // for CF
293 }
cb90a330 294
a0b94e5c 295 //checks if we have a prim vertex
296 if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
297 if(fHistograms != NULL){
298 fHistograms->FillHistogram("ESD_CutNContributors_InvMass",GetMotherCandidateMass());
299 }
300 fCurrentV0IndexNumber++;
301 continue;
302 }
1e7846f4 303 if(fDoCF){
304 fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors); // for CF
305 }
a0b94e5c 306
307 //Check the pid probability
308 if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
309 if(fHistograms != NULL){
310 fHistograms->FillHistogram("ESD_CutPIDProb_InvMass",GetMotherCandidateMass());
311 }
312 fCurrentV0IndexNumber++;
313 continue;
314 }
1e7846f4 315 if(fDoCF){
316 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID); // for CF
317 }
a0b94e5c 318
a0b94e5c 319 if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
320 if(fHistograms != NULL){
321 fHistograms->FillHistogram("ESD_CutR_InvMass",GetMotherCandidateMass());
322 }
323 fCurrentV0IndexNumber++;
324 continue;
1e7846f4 325 }
326 if(fDoCF){
327 fCFManager->GetParticleContainer()->Fill(containerInput,kStepR); // for CF
328 }
a0b94e5c 329
330
331 if((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue > GetXYRadius() ){ // cuts out regions where we do not reconstruct
332 if(fHistograms != NULL){
333 fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
334 }
335 fCurrentV0IndexNumber++;
336 continue;
1e7846f4 337 }
338 if(fDoCF){
339 fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine); // for CF
340 }
a0b94e5c 341
342 if(TMath::Abs(fCurrentZValue) > fMaxZ ){ // cuts out regions where we do not reconstruct
343 if(fHistograms != NULL){
344 fHistograms->FillHistogram("ESD_CutZ_InvMass",GetMotherCandidateMass());
345 }
346 fCurrentV0IndexNumber++;
347 continue;
1e7846f4 348 }
349 if(fDoCF){
350 fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ); // for CF
351 }
a0b94e5c 352
353 /* Moved further up so corr framework can work
354 if(UpdateV0Information() == kFALSE){
355 fCurrentV0IndexNumber++;
356 continue;
357 }
358 */
9640a3d1 359
a0b94e5c 360
361 if(fUseKFParticle){
362 if(fCurrentMotherKFCandidate->GetNDF()<=0){
363 if(fHistograms != NULL){
364 fHistograms->FillHistogram("ESD_CutNDF_InvMass",GetMotherCandidateMass());
365 }
366 fCurrentV0IndexNumber++;
367 continue;
368 }
1e7846f4 369 if(fDoCF){
370 fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF); // for CF
371 }
a0b94e5c 372
373 Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
374 if(chi2V0 > fChi2CutConversion || chi2V0 <=0){
375 if(fHistograms != NULL){
376 fHistograms->FillHistogram("ESD_CutChi2_InvMass",GetMotherCandidateMass());
377 }
378 fCurrentV0IndexNumber++;
379 continue;
380 }
1e7846f4 381 if(fDoCF){
382 fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2); // for CF
383 }
a0b94e5c 384
385 if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){
386 if(fHistograms != NULL){
387 fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
388 }
389 fCurrentV0IndexNumber++;
390 continue;
391 }
1e7846f4 392 if(fDoCF){
393 fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta); // for CF
394 }
a0b94e5c 395
396 if(fMotherCandidateLorentzVector->Pt()<fPtCut){
397 if(fHistograms != NULL){
398 fHistograms->FillHistogram("ESD_CutPt_InvMass",GetMotherCandidateMass());
399 }
400 fCurrentV0IndexNumber++;
401 continue;
402 }
1e7846f4 403 if(fDoCF){
404 fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt); // for CF
405 }
a0b94e5c 406
407 }
408 else if(fUseESDTrack){
409 //TODO
410 }
9640a3d1 411
412 if(fHistograms != NULL){
413 fHistograms->FillHistogram("ESD_GoodV0s_InvMass",GetMotherCandidateMass());
414 }
415
a0b94e5c 416 fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
417
418 iResult=kTRUE;//means we have a v0 who survived all the cuts applied
419
420 fCurrentV0IndexNumber++;
421
422 break;
423 }
424 return iResult;
425}
426
427Bool_t AliV0Reader::UpdateV0Information(){
428 //see header file for documentation
429
430 Bool_t iResult=kTRUE; // for taking out not refitted, kinks and like sign tracks
431
432 Bool_t switchTracks = kFALSE;
433
434 fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
435 fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
436
437 if(fCurrentNegativeESDTrack->GetSign() == fCurrentPositiveESDTrack->GetSign()){ // avoid like sign
438 iResult=kFALSE;
cb90a330 439 if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
a0b94e5c 440 fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
441 }
442 }
443
444 if(fCurrentPositiveESDTrack->GetSign() == -1 && fCurrentNegativeESDTrack->GetSign() == 1){ // switch wrong signed tracks
445 fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
446 fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
447 switchTracks = kTRUE;
448 }
449
450 if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ||
451 !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
452 // if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) ||
453 // !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
a0b94e5c 454 iResult=kFALSE;
cb90a330 455 if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
a0b94e5c 456 fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
457 }
458 }
459
460 if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
461 fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
462
463 iResult=kFALSE;
cb90a330 464 if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
a0b94e5c 465 fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
466 }
467 }
9640a3d1 468
469 if(fDodEdxSigmaCut == kTRUE){
470
471 if( fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
472 fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ||
473 fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
474 fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
475 iResult=kFALSE;
cb90a330 476 if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
9640a3d1 477 fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
478 }
479 }
480 if( fCurrentPositiveESDTrack->P()>fPIDMinPnSigmaAbovePionLine){
481 if(fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
482 fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
483 fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
484 iResult=kFALSE;
cb90a330 485 if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
9640a3d1 486 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
487 }
488 }
489 }
490
491 if( fCurrentNegativeESDTrack->P()>fPIDMinPnSigmaAbovePionLine){
492 if(fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
493 fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
494 fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
495 iResult=kFALSE;
cb90a330 496 if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
9640a3d1 497 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
498 }
499 }
500 }
501 }
a0b94e5c 502
503 if(fCurrentNegativeKFParticle != NULL){
504 delete fCurrentNegativeKFParticle;
505 }
506 if(switchTracks == kFALSE){
507 fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);
508 }
509 else{
510 fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fNegativeTrackPID);
511 }
512
513 if(fCurrentPositiveKFParticle != NULL){
514 delete fCurrentPositiveKFParticle;
515 }
516 if(switchTracks == kFALSE){
517 fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);
518 }
519 else{
520 fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fPositiveTrackPID);
521 }
522
523 if(fCurrentMotherKFCandidate != NULL){
524 delete fCurrentMotherKFCandidate;
525 }
526 fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
527
528
529 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
530 fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
531 }
532
533 if(fUseImprovedVertex == kTRUE){
534 AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
535 primaryVertexImproved+=*fCurrentMotherKFCandidate;
536 fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
537 }
538
539 fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);
540
541
542 if(fNegativeTrackLorentzVector != NULL){
543 delete fNegativeTrackLorentzVector;
544 }
545 if(fUseKFParticle){
546 fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());
547 }
548 else if(fUseESDTrack){
549 fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());
550 }
551
552 if(fPositiveTrackLorentzVector != NULL){
553 delete fPositiveTrackLorentzVector;
554 }
555 if(fUseKFParticle){
556 fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());
557 }
558 else if(fUseESDTrack){
559 fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());
560 }
561
562 if(fMotherCandidateLorentzVector != NULL){
563 delete fMotherCandidateLorentzVector;
564 }
565 if(fUseKFParticle){
566 fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
567 }
568 else if(fUseESDTrack){
569 fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
570 }
571
572 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
573 fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.);
574 }
575
576
577 if(fDoMC == kTRUE){
578 fMotherMCParticle= NULL;
579 fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
580 fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
581 if(fPositiveMCParticle->GetMother(0)>-1){
582 fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));
583 }
584 }
585
a0b94e5c 586
587 // for CF
588 Double_t containerInput[3];
1e7846f4 589 if(fDoCF){
590 containerInput[0] = GetMotherCandidatePt();
591 containerInput[1] = GetMotherCandidateEta();
592 containerInput[2] = GetMotherCandidateMass();
593
594 fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign); // for CF
595 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit); // for CF
596 fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks); // for CF
597 }
cb90a330 598
599
600 if(fUseOwnXYZCalculation == kFALSE){
601 fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
602 }
603 else{
604 Double_t convpos[2];
605 convpos[0]=0;
606 convpos[1]=0;
607 GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
608 fCurrentXValue = convpos[0];
609 fCurrentYValue = convpos[1];
610 fCurrentZValue = GetConvPosZ(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField());
611 }
612
613
614 fUpdateV0AlreadyCalled = kTRUE;
615
a0b94e5c 616 return iResult;
617}
618
619
620
621Bool_t AliV0Reader::HasSameMCMother(){
622 //see header file for documentation
623
624 Bool_t iResult = kFALSE;
625 if(fDoMC == kTRUE){
626 if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){
627 if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))
628 if(fMotherMCParticle){
629 iResult = kTRUE;
630 }
631 }
632 }
633 return iResult;
634}
635
636Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){
637 //see header file for documentation
638
639 Bool_t iResult=kFALSE;
640
641 Double_t *posProbArray = new Double_t[10];
642 Double_t *negProbArray = new Double_t[10];
643 AliESDtrack* negTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
644 AliESDtrack* posTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
645
646 negTrack->GetTPCpid(negProbArray);
647 posTrack->GetTPCpid(posProbArray);
648
4a6157dc 649 // if(negProbArray != NULL && posProbArray != NULL){ // this is not allowed anymore for some reason(RC19)
650 if(negProbArray && posProbArray){
a0b94e5c 651 if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){
652 iResult=kTRUE;
653 }
654 }
655 delete [] posProbArray;
656 delete [] negProbArray;
657 return iResult;
658}
659
660void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){
4a6157dc 661 // see header file for documentation
662
a0b94e5c 663 Double_t *posProbArray = new Double_t[10];
664 Double_t *negProbArray = new Double_t[10];
665 AliESDtrack* negTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
666 AliESDtrack* posTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
667
668 negTrack->GetTPCpid(negProbArray);
669 posTrack->GetTPCpid(posProbArray);
670
4a6157dc 671 // if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
672 if(negProbArray && posProbArray){
a0b94e5c 673 negPIDProb = negProbArray[GetSpeciesIndex(-1)];
674 posPIDProb = posProbArray[GetSpeciesIndex(1)];
675 }
676 delete [] posProbArray;
677 delete [] negProbArray;
678}
679
680void AliV0Reader::UpdateEventByEventData(){
681 //see header file for documentation
682
683 if(fCurrentEventGoodV0s.size() >0 ){
684 // fPreviousEventGoodV0s.clear();
685 // fPreviousEventGoodV0s = fCurrentEventGoodV0s;
686 if(fPreviousEventGoodV0s.size()>19){
687 for(UInt_t nCurrent=0;nCurrent<fCurrentEventGoodV0s.size();nCurrent++){
688 fPreviousEventGoodV0s.erase(fPreviousEventGoodV0s.begin());
689 fPreviousEventGoodV0s.push_back(fCurrentEventGoodV0s.at(nCurrent));
690 }
691 }
692 else{
693 for(UInt_t nCurrent=0;nCurrent<fCurrentEventGoodV0s.size();nCurrent++){
694 if(fPreviousEventGoodV0s.size()<20){
695 fPreviousEventGoodV0s.push_back(fCurrentEventGoodV0s.at(nCurrent));
696 }
697 else{
698 fPreviousEventGoodV0s.erase(fPreviousEventGoodV0s.begin());
699 fPreviousEventGoodV0s.push_back(fCurrentEventGoodV0s.at(nCurrent));
700 }
701 }
702 }
703 }
704 fCurrentEventGoodV0s.clear();
705
706 fCurrentV0IndexNumber=0;
707}
708
709
710Double_t AliV0Reader::GetNegativeTrackPhi() const{
711 //see header file for documentation
712
713 Double_t offset=0;
714 if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
715 offset = -2*TMath::Pi();
716 }
717 return fNegativeTrackLorentzVector->Phi()+offset;
718}
719
720Double_t AliV0Reader::GetPositiveTrackPhi() const{
721 //see header file for documentation
722
723 Double_t offset=0;
724 if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
725 offset = -2*TMath::Pi();
726 }
727 return fPositiveTrackLorentzVector->Phi()+offset;
728}
729
730Double_t AliV0Reader::GetMotherCandidatePhi() const{
731 //see header file for documentation
732
733 Double_t offset=0;
734 if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
735 offset = -2*TMath::Pi();
736 }
737 return fMotherCandidateLorentzVector->Phi()+offset;
738}
739
740
741Double_t AliV0Reader::GetMotherCandidateRapidity() const{
742 //see header file for documentation
743
744 Double_t rapidity=0;
745 if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
746 else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
747 return rapidity;
748
749}
750
751
752
753
754
755Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
756 //see header file for documentation
757
758 Int_t iResult = 10; // Unknown particle
759
760 if(chargeOfTrack==-1){ //negative track
761 switch(abs(fNegativeTrackPID)){
762 case 11: //electron
763 iResult = 0;
764 break;
765 case 13: //muon
766 iResult = 1;
767 break;
768 case 211: //pion
769 iResult = 2;
770 break;
771 case 321: //kaon
772 iResult = 3;
773 break;
774 case 2212: //proton
775 iResult = 4;
776 break;
777 case 22: //photon
778 iResult = 5;
779 break;
780 case 111: //pi0
781 iResult = 6;
782 break;
783 case 2112: //neutron
784 iResult = 7;
785 break;
786 case 311: //K0
787 iResult = 8;
788 break;
789
790 //Put in here for kSPECIES::kEleCon ????
791 }
792 }
793 else if(chargeOfTrack==1){ //positive track
794 switch(abs(fPositiveTrackPID)){
795 case 11: //electron
796 iResult = 0;
797 break;
798 case 13: //muon
799 iResult = 1;
800 break;
801 case 211: //pion
802 iResult = 2;
803 break;
804 case 321: //kaon
805 iResult = 3;
806 break;
807 case 2212: //proton
808 iResult = 4;
809 break;
810 case 22: //photon
811 iResult = 5;
812 break;
813 case 111: //pi0
814 iResult = 6;
815 break;
816 case 2112: //neutron
817 iResult = 7;
818 break;
819 case 311: //K0
820 iResult = 8;
821 break;
822
823 //Put in here for kSPECIES::kEleCon ????
824 }
825 }
826 else{
827 //Wrong parameter.. Print warning
828 }
829 return iResult;
830}
831
4a6157dc 832Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
a0b94e5c 833 // see header file for documentation
834
835 Double_t pi = 3.14159265358979323846;
836
837 Double_t helix[6];
838 track->GetHelixParameters(helix,b);
839
840 Double_t xpos = helix[5];
841 Double_t ypos = helix[0];
842 Double_t radius = TMath::Abs(1./helix[4]);
843 Double_t phi = helix[2];
844
845 if(phi < 0){
846 phi = phi + 2*pi;
847 }
848
849 phi -= pi/2.;
850 Double_t xpoint = radius * TMath::Cos(phi);
851 Double_t ypoint = radius * TMath::Sin(phi);
852
853 if(charge > 0){
854 xpoint = - xpoint;
855 ypoint = - ypoint;
856 }
857
858 if(charge < 0){
859 xpoint = xpoint;
860 ypoint = ypoint;
861 }
862 center[0] = xpos + xpoint;
863 center[1] = ypos + ypoint;
864
865 return 1;
866}
867
4a6157dc 868Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
a0b94e5c 869 //see header file for documentation
870
871 Double_t helixcenterpos[2];
872 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
873
874 Double_t helixcenterneg[2];
875 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
876
877 Double_t poshelix[6];
878 ptrack->GetHelixParameters(poshelix,b);
879 Double_t posradius = TMath::Abs(1./poshelix[4]);
880
881 Double_t neghelix[6];
882 ntrack->GetHelixParameters(neghelix,b);
883 Double_t negradius = TMath::Abs(1./neghelix[4]);
884
885 Double_t xpos = helixcenterpos[0];
886 Double_t ypos = helixcenterpos[1];
887 Double_t xneg = helixcenterneg[0];
888 Double_t yneg = helixcenterneg[1];
889
890 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
891 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
892
893 return 1;
894}
895
896
897
4a6157dc 898Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
a0b94e5c 899 //see header file for documentation
900
901 Double_t helixpos[6];
902 ptrack->GetHelixParameters(helixpos,b);
903
904 Double_t helixneg[6];
905 ntrack->GetHelixParameters(helixneg,b);
906
907 Double_t negtrackradius = TMath::Abs(1./helixneg[4]);
908 Double_t postrackradius = TMath::Abs(1./helixpos[4]);
909
910 Double_t pi = 3.14159265358979323846;
911
912 Double_t convpos[2];
913 GetConvPosXY(ptrack,ntrack,b,convpos);
914
915 Double_t convposx = convpos[0];
916 Double_t convposy = convpos[1];
917
918 Double_t helixcenterpos[2];
919 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
920
921 Double_t helixcenterneg[2];
922 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
923
924 Double_t xpos = helixcenterpos[0];
925 Double_t ypos = helixcenterpos[1];
926 Double_t xneg = helixcenterneg[0];
927 Double_t yneg = helixcenterneg[1];
928
4a6157dc 929 Double_t deltaXPos = convposx - xpos;
930 Double_t deltaYPos = convposy - ypos;
a0b94e5c 931
4a6157dc 932 Double_t deltaXNeg = convposx - xneg;
933 Double_t deltaYNeg = convposy - yneg;
a0b94e5c 934
4a6157dc 935 Double_t alphaPos = pi + TMath::ATan2(-deltaYPos,-deltaXPos);
936 Double_t alphaNeg = pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
a0b94e5c 937
4a6157dc 938 Double_t vertexXNeg = xneg + TMath::Abs(negtrackradius)*
939 TMath::Cos(alphaNeg);
940 Double_t vertexYNeg = yneg + TMath::Abs(negtrackradius)*
941 TMath::Sin(alphaNeg);
a0b94e5c 942
4a6157dc 943 Double_t vertexXPos = xpos + TMath::Abs(postrackradius)*
944 TMath::Cos(alphaPos);
945 Double_t vertexYPos = ypos + TMath::Abs(postrackradius)*
946 TMath::Sin(alphaPos);
a0b94e5c 947
948 Double_t x0neg = helixneg[5];
949 Double_t y0neg = helixneg[0];
950
951 Double_t x0pos = helixpos[5];
952 Double_t y0pos = helixpos[0];
953
4a6157dc 954 Double_t dNeg = TMath::Sqrt((vertexXNeg - x0neg)*(vertexXNeg - x0neg)
955 +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
a0b94e5c 956
4a6157dc 957 Double_t dPos = TMath::Sqrt((vertexXPos - x0pos)*(vertexXPos - x0pos)
958 +(vertexYPos - y0pos)*(vertexYPos - y0pos));
a0b94e5c 959
4a6157dc 960 Double_t rNeg = TMath::Sqrt(negtrackradius*negtrackradius -
961 dNeg*dNeg/4.);
a0b94e5c 962
4a6157dc 963 Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
964 dPos*dPos/4.);
a0b94e5c 965
4a6157dc 966 Double_t deltabetaNeg = 2*(pi + TMath::ATan2(-dNeg/2.,-rNeg));
967 Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
a0b94e5c 968
4a6157dc 969 Double_t deltaUNeg = negtrackradius*deltabetaNeg;
970 Double_t deltaUPos = postrackradius*deltabetaPos;
a0b94e5c 971
4a6157dc 972 Double_t zphaseNeg = ntrack->GetZ() + deltaUNeg * ntrack->GetTgl();
973 Double_t zphasePos = ptrack->GetZ() + deltaUPos * ptrack->GetTgl();
a0b94e5c 974
975 Double_t convposz =
4a6157dc 976 (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
a0b94e5c 977
978 return convposz;
979}