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