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