]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGHF/hfe/AliHFEsecVtx.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEsecVtx.cxx
... / ...
CommitLineData
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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// Secondary vertexing construction Class
17// Construct secondary vertex from Beauty hadron with electron and
18// hadrons, then apply selection criteria
19//
20// Authors:
21// MinJung Kweon <minjung@physi.uni-heidelberg.de>
22//
23
24#include <TH2F.h>
25#include <TIterator.h>
26#include <TParticle.h>
27
28#include <AliESDVertex.h>
29#include <AliESDEvent.h>
30#include <AliAODEvent.h>
31#include <AliVTrack.h>
32#include <AliESDtrack.h>
33#include <AliAODTrack.h>
34#include "AliHFEsecVtx.h"
35#include <AliKFParticle.h>
36#include <AliKFVertex.h>
37#include <AliLog.h>
38#include <AliMCEvent.h>
39#include <AliAODMCParticle.h>
40#include "AliHFEpairs.h"
41#include "AliHFEsecVtxs.h"
42#include "AliHFEtrackFilter.h"
43#include "AliHFEmcQA.h"
44#include "AliHFEtools.h"
45
46ClassImp(AliHFEsecVtx)
47
48//_______________________________________________________________________________________________
49AliHFEsecVtx::AliHFEsecVtx():
50 fFilter(0x0)
51 ,fESD1(0x0)
52 ,fAOD1(0x0)
53 ,fMCEvent(0x0)
54 ,fMCQA(0x0)
55 ,fUseMCPID(kFALSE)
56 ,fkSourceLabel()
57 ,fNparents(0)
58 ,fParentSelect()
59 ,fPtRng()
60 ,fDcaCut()
61 ,fNoOfHFEpairs(0)
62 ,fNoOfHFEsecvtxs(0)
63 ,fArethereSecVtx(0)
64 ,fHFEpairs(0x0)
65 ,fHFEsecvtxs(0x0)
66 ,fMCArray(0x0)
67 ,fPVx(-999)
68 ,fPVy(-999)
69 ,fPVx2(999)
70 ,fPVy2(-999)
71 ,fCosPhi(-1)
72 ,fSignedLxy(-1)
73 ,fSignedLxy2(-1)
74 ,fKFchi2(-1)
75 ,fInvmass(-1)
76 ,fInvmassSigma(-1)
77 ,fKFip(0)
78 ,fKFip2(0)
79 ,fNsectrk2prim(0)
80 ,fVtxchi2Tightcut(3.)
81 ,fVtxchi2Loosecut(5.)
82 ,fPairQA(0x0)
83 ,fSecvtxQA(0x0)
84 ,fSecVtxList(0x0)
85{
86 //
87 // Default constructor
88 //
89
90 Init();
91}
92
93//_______________________________________________________________________________________________
94AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
95 TObject(p)
96 ,fFilter(0x0)
97 ,fESD1(0x0)
98 ,fAOD1(0x0)
99 ,fMCEvent(0x0)
100 ,fMCQA(0x0)
101 ,fUseMCPID(p.fUseMCPID)
102 ,fkSourceLabel()
103 ,fNparents(p.fNparents)
104 ,fParentSelect()
105 ,fPtRng()
106 ,fDcaCut()
107 ,fNoOfHFEpairs(p.fNoOfHFEpairs)
108 ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
109 ,fArethereSecVtx(p.fArethereSecVtx)
110 ,fHFEpairs(0x0)
111 ,fHFEsecvtxs(0x0)
112 ,fMCArray(0x0)
113 ,fPVx(p.fPVx)
114 ,fPVy(p.fPVy)
115 ,fPVx2(p.fPVx2)
116 ,fPVy2(p.fPVy2)
117 ,fCosPhi(p.fCosPhi)
118 ,fSignedLxy(p.fSignedLxy)
119 ,fSignedLxy2(p.fSignedLxy2)
120 ,fKFchi2(p.fKFchi2)
121 ,fInvmass(p.fInvmass)
122 ,fInvmassSigma(p.fInvmassSigma)
123 ,fKFip(p.fKFip)
124 ,fKFip2(p.fKFip2)
125 ,fNsectrk2prim(p.fNsectrk2prim)
126 ,fVtxchi2Tightcut(p.fVtxchi2Tightcut)
127 ,fVtxchi2Loosecut(p.fVtxchi2Loosecut)
128 ,fPairQA(0x0)
129 ,fSecvtxQA(0x0)
130 ,fSecVtxList(0x0)
131{
132 //
133 // Copy constructor
134 //
135 fFilter = new AliHFEtrackFilter(*p.fFilter);
136}
137
138//_______________________________________________________________________________________________
139AliHFEsecVtx&
140AliHFEsecVtx::operator=(const AliHFEsecVtx &)
141{
142 //
143 // Assignment operator
144 //
145
146 AliInfo("Not yet implemented.");
147 return *this;
148}
149
150//_______________________________________________________________________________________________
151AliHFEsecVtx::~AliHFEsecVtx()
152{
153 //
154 // Destructor
155 //
156
157 //cout << "Analysis Done." << endl;
158 delete fFilter;
159}
160
161//__________________________________________
162void AliHFEsecVtx::Init()
163{
164 //
165 // set pdg code and index
166 //
167
168 fNparents = 7;
169
170 fParentSelect[0][0] = 411;
171 fParentSelect[0][1] = 421;
172 fParentSelect[0][2] = 431;
173 fParentSelect[0][3] = 4122;
174 fParentSelect[0][4] = 4132;
175 fParentSelect[0][5] = 4232;
176 fParentSelect[0][6] = 4332;
177
178 fParentSelect[1][0] = 511;
179 fParentSelect[1][1] = 521;
180 fParentSelect[1][2] = 531;
181 fParentSelect[1][3] = 5122;
182 fParentSelect[1][4] = 5132;
183 fParentSelect[1][5] = 5232;
184 fParentSelect[1][6] = 5332;
185
186 // momentum ranges to apply pt dependent cuts
187 fPtRng[0] = 1.0;
188 fPtRng[1] = 2.0;
189 fPtRng[2] = 2.5;
190 fPtRng[3] = 3.0;
191 fPtRng[4] = 5.0;
192 fPtRng[5] = 12.0;
193 fPtRng[6] = 20.0;
194
195 // momentum dependent DCA cut values (preliminary)
196 fDcaCut[0] = 0.04;
197 fDcaCut[1] = 0.03;
198 fDcaCut[2] = 0.02;
199 fDcaCut[3] = 0.015;
200 fDcaCut[4] = 0.01;
201 fDcaCut[5] = 0.005;
202
203 fFilter = new AliHFEtrackFilter("SecVtxFilter");
204 fFilter->MakeCutStepRecKineITSTPC();
205 fFilter->MakeCutStepPrimary();
206}
207
208Bool_t AliHFEsecVtx::Process(AliVTrack *signalTrack){
209 //
210 // Run Process
211 //
212 //if(signalTrack->Pt() < 1.0) return;
213
214
215 AliESDtrack *track = dynamic_cast<AliESDtrack *>(signalTrack);
216
217 if(!track){
218 AliDebug(1, "no esd track pointer, return\n");
219 return kFALSE;
220 }
221
222 Bool_t bTagged = kFALSE;
223
224 FillHistos(0,track); // wo any cuts
225
226 InitHFEpairs();
227 InitHFEsecvtxs();
228 AliESDtrack *htrack = 0x0;
229 fFilter->Flush();
230 fFilter->FilterTracks(fESD1);
231 TObjArray *candidates = fFilter->GetFilteredTracks();
232 TIterator *trackIter = candidates->MakeIterator();
233 while((htrack = dynamic_cast<AliESDtrack *>(trackIter->Next()))){
234 if(track->GetID() == htrack->GetID()) continue; // since it is for tagging single electron, don't need additional condition
235 if (htrack->Pt()<1.0) continue;
236 PairAnalysis(track, htrack, htrack->GetID()); // e-h pairing
237 }
238 delete trackIter;
239 for(int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
240 //if(HasMCData()){
241 AliHFEpairs *pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
242 //if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.)) // apply various cuts
243 if(pair->GetKFChi2()>5.) // only apply vertex chi2 cut for the moment
244 HFEpairs()->RemoveAt(ip);
245 //}
246 }
247 HFEpairs()->Compress();
248 if(HFEpairs()->GetEntriesFast()) FillHistos(1,track); // after paired
249 if(HFEpairs()->GetEntriesFast()) RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
250 if(HFEsecvtxs()->GetEntriesFast()) FillHistos(2,track); // after secvtxing
251 for(int ip=0; ip<HFEsecvtxs()->GetEntriesFast(); ip++){
252 AliHFEsecVtxs *secvtx=0x0;
253 secvtx = (AliHFEsecVtxs*) (HFEsecvtxs()->UncheckedAt(ip));
254 // here you apply cuts, then if it doesn't pass the cut, remove it from the HFEsecvtxs()
255 if(!(secvtx->GetInvmass()>2.0 && secvtx->GetInvmass()<5.2) || !(secvtx->GetSignedLxy2()>0.08 && secvtx->GetSignedLxy2()<1.5) || !(secvtx->GetKFIP2()>-0.1 && secvtx->GetKFIP2()<0.1))
256 HFEsecvtxs()->RemoveAt(ip);
257 }
258
259 // fill histos for raw spectra
260 if(HFEsecvtxs()->GetEntriesFast()) {
261 FillHistos(3,track); //after secvtx cut
262 bTagged = kTRUE;
263 }
264
265 DeleteHFEpairs();
266 DeleteHFEsecvtxs();
267
268 return bTagged;
269}
270
271//_______________________________________________________________________________________________
272void AliHFEsecVtx::CreateHistograms(TList * const qaList)
273{
274 //
275 // create histograms
276 //
277
278 if(!qaList) return;
279
280 fSecVtxList = qaList;
281 fSecVtxList->SetName("SecVtx");
282
283 MakeContainer();
284 MakeHistos(0);
285 MakeHistos(1);
286 MakeHistos(2);
287 MakeHistos(3);
288
289 /*
290 fkSourceLabel[kAll]="all";
291 fkSourceLabel[kDirectCharm]="directCharm";
292 fkSourceLabel[kDirectBeauty]="directBeauty";
293 fkSourceLabel[kBeautyCharm]="beauty2charm";
294 fkSourceLabel[kGamma]="gamma";
295 fkSourceLabel[kPi0]="pi0";
296 fkSourceLabel[kElse]="others";
297 fkSourceLabel[kBeautyGamma]="beauty22gamma";
298 fkSourceLabel[kBeautyPi0]="beauty22pi0";
299 fkSourceLabel[kBeautyElse]="beauty22others";
300
301 TString hname;
302 TString hnopt="secvtx_";
303 for (Int_t isource = 0; isource < 10; isource++ ){
304 }
305 */
306
307 //qaList->AddLast(fSecVtxList);
308}
309
310//_______________________________________________________________________________________________
311void AliHFEsecVtx::GetPrimaryCondition()
312{
313 //
314 // get primary characteristics and set
315 //
316
317 if (fESD1) {
318 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
319 if( primVtxCopy.GetNDF() <1 ) return;
320 fPVx = primVtxCopy.GetX();
321 fPVy = primVtxCopy.GetY();
322 }
323 else if(fAOD1) {
324 AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
325 if( primVtxCopy.GetNDF() <1 ) return;
326 fPVx = primVtxCopy.GetX();
327 fPVy = primVtxCopy.GetY();
328 }
329}
330
331//_______________________________________________________________________________________________
332void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
333{
334 //
335 // calculate e-h pair characteristics and tag pair
336 //
337
338 //later consider the below
339 Float_t dca1[2]={-999.,-999.}, dca2[2]={-999.,-999.};
340 Float_t cov1[3]={-999.,-999.,-999.}, cov2[3]={-999.,-999.,-999.};
341
342 Double_t dca1aod[2]={-999.,-999.}, dca2aod[2]={-999.,-999.};
343 Double_t cov1aod[3]={-999.,-999.,-999.}, cov2aod[3]={-999.,-999.,-999.};
344
345 if (IsAODanalysis()){
346 const AliAODVertex *primVtx = fAOD1->GetPrimaryVertex();
347 AliESDtrack esdTrk1(track1);
348 AliESDtrack esdTrk2(track2);
349 esdTrk1.PropagateToDCA(primVtx,0.,10000.,dca1aod,cov1aod);
350 esdTrk2.PropagateToDCA(primVtx,0.,10000.,dca2aod,cov2aod);
351 }
352 else {
353 ((AliESDtrack*)track1)->GetImpactParameters(dca1,cov1);
354 ((AliESDtrack*)track2)->GetImpactParameters(dca2,cov2);
355 }
356
357 // apply pt dependent dca cut on hadrons
358 /*for(int ibin=0; ibin<6; ibin++){
359 if((track2->Pt()>fPtRng[ibin] && track2->Pt()<fPtRng[ibin+1]) && TMath::Abs(dca2[0])<fDcaCut[ibin]) return;
360 }*/
361
362 // get KF particle input pid
363 Int_t pdg1 = GetPDG(track1);
364 Int_t pdg2 = GetPDG(track2);
365
366 if(pdg1==-1 || pdg2==-1) {
367 //printf("out if considered pid range \n");
368 return;
369 }
370
371 // create KF particle of pair
372 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
373 else AliKFParticle::SetField(fESD1->GetMagneticField());
374
375 AliKFParticle kfTrack[2];
376 kfTrack[0] = AliKFParticle(*track1, pdg1);
377 kfTrack[1] = AliKFParticle(*track2, pdg2);
378
379 AliKFParticle kfSecondary(kfTrack[0],kfTrack[1]);
380
381 //secondary vertex point from kf particle
382 Double_t kfx = kfSecondary.GetX();
383 Double_t kfy = kfSecondary.GetY();
384 //Double_t kfz = kfSecondary.GetZ();
385
386 //momentum at the decay point from kf particle
387 Double_t kfpx = kfSecondary.GetPx();
388 Double_t kfpy = kfSecondary.GetPy();
389 //Double_t kfpz = kfSecondary.GetPz();
390
391
392/* //directly use of ESD vertex
393 const AliESDVertex *pvertex = fESD1->GetPrimaryVertex();
394 Double_t xyzVtx[3];
395 pvertex->GetXYZ(xyzVtx);
396
397 Double_t dx = kfx-xyzVtx[0];
398 Double_t dy = kfy-xyzVtx[1];*/
399
400 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
401 if( primVtxCopy.GetNDF() <1 ) return;
402 fPVx = primVtxCopy.GetX();
403 fPVy = primVtxCopy.GetY();
404
405 // printf("esdx= %lf kfx= %lf esdy= %lf kfy= %lf\n",xyzVtx[0],fPVx,xyzVtx[1],fPVy);
406
407 Double_t dx = kfx-fPVx;
408 Double_t dy = kfy-fPVy;
409
410 // discriminating variables
411 // invariant mass of the KF particle
412 Double_t invmass = -1;
413 Double_t invmassSigma = -1;
414 kfSecondary.GetMass(invmass,invmassSigma);
415
416 // chi2 of the KF particle
417 Double_t kfchi2 = -1;
418 if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
419
420 // opening angle between two particles in XY plane
421 Double_t phi = kfTrack[0].GetAngleXY(kfTrack[1]);
422 Double_t cosphi = TMath::Cos(phi);
423
424 // DCA from primary to e-h KF particle (impact parameter of KF particle)
425 Double_t vtx[2]={fPVx, fPVy};
426 Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
427
428 // projection of kf vertex vector to the kf momentum direction
429 Double_t signedLxy=-999.;
430 if((dx*kfpx+dy*kfpy)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
431 if((dx*kfpx+dy*kfpy)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
432 //[the other way to think about]
433 //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
434 //if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
435
436 //recalculating primary vertex after removing secvtx tracks --------------------------
437 Int_t trkid[2];
438 trkid[0] = track1->GetID();
439 trkid[1] = track2->GetID();
440
441 RecalcPrimvtx(2, trkid, kfTrack);
442 Double_t dx2 = kfx-fPVx2;
443 Double_t dy2 = kfy-fPVy2;
444
445 // IP of sec particle recalculated based on recalculated primary vertex
446 Double_t vtx2[2]={fPVx2, fPVy2};
447 Double_t kfip2 = kfSecondary.GetDistanceFromVertexXY(vtx2);
448 // signed Lxy recalculated based on recalculated primary vertex
449 Double_t signedLxy2=-999.;
450 if((dx2*kfpx+dy2*kfpy)>0) signedLxy2 = TMath::Sqrt(dx2*dx2+dy2*dy2);
451 if((dx2*kfpx+dy2*kfpy)<0) signedLxy2 = -1*TMath::Sqrt(dx2*dx2+dy2*dy2);
452 //------------------------------------------------------------------------------------
453
454
455 Int_t paircode = -1;
456 if (HasMCData()) paircode = GetPairCode(track1,track2);
457
458 AliHFEpairs hfepair;
459 hfepair.SetTrkLabel(index2);
460 hfepair.SetInvmass(invmass);
461 hfepair.SetKFChi2(kfchi2);
462 hfepair.SetOpenangle(phi);
463 hfepair.SetCosOpenangle(cosphi);
464 hfepair.SetSignedLxy(signedLxy);
465 hfepair.SetSignedLxy2(signedLxy2);
466 hfepair.SetKFIP(kfip);
467 hfepair.SetKFIP2(kfip2);
468 hfepair.SetPairCode(paircode);
469 AddHFEpairToArray(&hfepair);
470 fNoOfHFEpairs++;
471
472 // fill into container for later QA
473 Double_t dataE[6];
474 dataE[0]=invmass;
475 dataE[1]=kfchi2;
476 dataE[2]=phi;
477 dataE[3]=signedLxy2;
478 dataE[4]=kfip2;
479 dataE[5]=paircode;
480 //if(cov1[0]>0) dataE[6]=Double_t(dca1[0]/cov1[0]);
481 //if(cov2[0]>0) dataE[7]=Double_t(dca2[0]/cov2[0]);
482 //dataE[6]=track1->Pt();
483 //dataE[7]=track2->Pt();
484 //dataE[6]=dca1[0]; //mjtmp
485 //dataE[7]=dca2[0]; //mjtmp
486 //dataE[8]=TMath::Abs(dca1[0]);
487 //dataE[9]=TMath::Abs(dca2[0]);
488
489 fPairQA->Fill(dataE);
490
491}
492
493//_______________________________________________________________________________________________
494void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
495{
496 //
497 // run secondary vertexing algorithm and do tagging
498 //
499
500 //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
501 FindSECVTXCandid(track);
502}
503
504//_______________________________________________________________________________________________
505void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
506{
507 //
508 // find secondary vertex candidate and store them into container
509 //
510
511 AliVTrack *htrack[20];
512 //Int_t htracklabel[20];
513 //Int_t paircode[20];
514 //Double_t vtxchi2[20];
515 //Double_t dataE[7]={-999.,-999.,-999.,-999.,-1.,0,0};
516
517 fVtxchi2Tightcut=3.; // tight cut for pair
518 fVtxchi2Loosecut=5.; // loose cut for secvtx
519
520 if (HFEpairs()->GetEntriesFast()>20){
521 AliDebug(3, "number of paired hadron is over maximum(20)");
522 return;
523 }
524
525 // get paired track objects
526 AliHFEpairs *pair=0x0;
527 if (IsAODanalysis()){
528 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
529 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
530 //htracklabel[ip] = pair->GetTrkLabel();
531 htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
532 //if(pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode[ip]=1;
533 //else paircode[ip]=0;
534 //vtxchi2[ip] = pair->GetKFChi2();
535 }
536 }
537 else{
538 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
539 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
540 //htracklabel[ip] = pair->GetTrkLabel();
541 htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
542 //if(pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode[ip]=1;
543 //else paircode[ip]=0;
544 //vtxchi2[ip] = pair->GetKFChi2();
545 }
546 }
547
548 Int_t nPairs = HFEpairs()->GetEntriesFast();
549
550
551 // 1 electron candidate + 1 track
552 if (nPairs == 1){
553 if (pair->GetKFChi2() < fVtxchi2Tightcut) { // you can also put single track cut -> here you apply very tight cut for the pair
554 Fill2TrkSECVTX(track, pair);
555 }
556 return;
557 }
558 //--------------------------------------------------------------
559
560 // 1 electron candidate + 2 tracks
561 if (nPairs == 2){
562 CalcSECVTXProperty(track, htrack[0], htrack[1]); // calculate secondary vertex property
563
564 if (fKFchi2 < fVtxchi2Loosecut) { // -> here you apply rather loose cut
565 Fill3TrkSECVTX(track, 0, 1);
566 }
567 else{ // if doesn't pass the sec vtx chi2 cut
568 for(int jp=0; jp<2; jp++){
569 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jp));
570 if (pair->GetKFChi2() < fVtxchi2Tightcut){
571 Fill2TrkSECVTX(track, pair);
572 }
573 }
574 }
575 return;
576 }
577 //--------------------------------------------------------------
578
579 // 1 electron candidate + 3 tracks
580 if (nPairs == 3){
581 CalcSECVTXProperty(track, htrack[0], htrack[1], htrack[2]); // calculate secondary vertex property
582
583 if (fKFchi2 < fVtxchi2Loosecut) {
584 Fill4TrkSECVTX(track, 0, 1, 2);
585 }
586 else {
587 fArethereSecVtx=0;
588 for (int i=0; i<nPairs-1; i++){
589 for (int j=i+1; j<nPairs; j++){
590 CalcSECVTXProperty(track, htrack[i], htrack[j]);
591 if (fKFchi2 < fVtxchi2Loosecut) {
592 fArethereSecVtx++;
593 Fill3TrkSECVTX(track, i, j);
594 }
595 }
596 }
597 if(!fArethereSecVtx){
598 for(int jp=0; jp<nPairs; jp++){
599 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jp));
600 if (pair->GetKFChi2() < fVtxchi2Tightcut){
601 Fill2TrkSECVTX(track, pair);
602 }
603 }
604 }
605 }
606 return;
607 }
608 //--------------------------------------------------------------
609
610 // 1 electron candidate + more than 3 tracks
611 if (nPairs > 3){
612 fArethereSecVtx=0;
613 for (int ih1=0; ih1<nPairs-2; ih1++){
614 for (int ih2=ih1+1; ih2<nPairs-1; ih2++){
615 for (int ih3=ih2+1; ih3<nPairs; ih3++){
616 CalcSECVTXProperty(track, htrack[ih1], htrack[ih2], htrack[ih3]); // calculate secondary vertex property
617 if (fKFchi2 < fVtxchi2Loosecut) {
618 fArethereSecVtx++;
619 Fill4TrkSECVTX(track, ih1, ih2, ih3);
620 }
621 }
622 }
623 }
624 if (!fArethereSecVtx){
625 fArethereSecVtx=0;
626 for (int i=0; i<nPairs-1; i++){
627 for (int j=i+1; j<nPairs; j++){
628 CalcSECVTXProperty(track, htrack[i], htrack[j]);
629 if (fKFchi2 < fVtxchi2Loosecut) {
630 fArethereSecVtx++;
631 Fill3TrkSECVTX(track, i, j);
632 }
633 }
634 }
635 }
636 if (!fArethereSecVtx){
637 for(int jp=0; jp<nPairs; jp++){
638 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jp));
639 if (pair->GetKFChi2() < fVtxchi2Tightcut){
640 Fill2TrkSECVTX(track, pair);
641 }
642 }
643 }
644 return;
645 }
646 //--------------------------------------------------------------
647
648}
649
650//_______________________________________________________________________________________________
651void AliHFEsecVtx::Fill4TrkSECVTX(AliVTrack* track, Int_t ipair, Int_t jpair, Int_t kpair)
652{
653 //
654 // fill 3 tracks' secondary vertex properties
655 //
656
657 Double_t dataE[9]={-999.,-999.,-999.,-999.,-1.,0,0,-999.,-999.};
658
659 Int_t paircode1 = 0, paircode2 = 0, paircode3 = 0;
660 Int_t htracklabel1 = 0, htracklabel2= 0;
661
662 if (HasMCData()){
663 AliHFEpairs *pair1=0x0;
664 AliHFEpairs *pair2=0x0;
665 AliHFEpairs *pair3=0x0;
666 pair1 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ipair));
667 pair2 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jpair));
668 pair3 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(kpair));
669
670 htracklabel1 = pair1->GetTrkLabel();
671 htracklabel2 = pair2->GetTrkLabel();
672
673 if (pair1->GetPairCode()==2 || pair1->GetPairCode()==3) paircode1=1;
674 else paircode1=0;
675 if (pair2->GetPairCode()==2 || pair2->GetPairCode()==3) paircode2=1;
676 else paircode2=0;
677 if (pair3->GetPairCode()==2 || pair3->GetPairCode()==3) paircode3=1;
678 else paircode3=0;
679 }
680
681 AliHFEsecVtxs hfesecvtx;
682 hfesecvtx.SetTrkLabel1(htracklabel1); // mj: not much meaningful for the moment
683 hfesecvtx.SetTrkLabel2(htracklabel2); // mj: not much meaningful for the moment
684 if(HasMCData()) hfesecvtx.SetMCCode(GetElectronSource(TMath::Abs(track->GetLabel())));
685 hfesecvtx.SetKFChi2(fKFchi2);
686 hfesecvtx.SetInvmass(fInvmass);
687 hfesecvtx.SetSignedLxy(fSignedLxy);
688 hfesecvtx.SetSignedLxy2(fSignedLxy2);
689 hfesecvtx.SetKFIP(fKFip);
690 hfesecvtx.SetKFIP2(fKFip2);
691 AddHFEsecvtxToArray(&hfesecvtx);
692 fNoOfHFEsecvtxs++;
693
694 dataE[0]=fInvmass;
695 dataE[1]=fKFchi2;
696 dataE[2]=fSignedLxy;
697 dataE[3]=fKFip;
698 if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
699 dataE[5]=4; //# of associated tracks
700 if(paircode1 & paircode2 & paircode3) dataE[6]=1;
701 else if(!paircode1 & !paircode2 & !paircode3) dataE[6]=0;
702 else dataE[6]=3;
703 dataE[7]=fSignedLxy2;
704 dataE[8]=track->Pt();
705 fSecvtxQA->Fill(dataE);
706}
707
708//_______________________________________________________________________________________________
709void AliHFEsecVtx::Fill3TrkSECVTX(AliVTrack* track, Int_t ipair, Int_t jpair)
710{
711 //
712 // fill 3 tracks' secondary vertex properties
713 //
714
715 Double_t dataE[9]={-999.,-999.,-999.,-999.,-1.,0,0,-999.,-999.};
716
717 Int_t paircode1 = 0, paircode2 = 0;
718 Int_t htracklabel1 = 0, htracklabel2 = 0;
719
720 if (HasMCData()){
721 AliHFEpairs *pair1=0x0;
722 AliHFEpairs *pair2=0x0;
723 pair1 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ipair));
724 pair2 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jpair));
725
726 htracklabel1 = pair1->GetTrkLabel();
727 htracklabel2 = pair2->GetTrkLabel();
728
729 if (pair1->GetPairCode()==2 || pair1->GetPairCode()==3) paircode1=1;
730 else paircode1=0;
731 if (pair2->GetPairCode()==2 || pair2->GetPairCode()==3) paircode2=1;
732 else paircode2=0;
733 }
734
735 // fill secondary vertex container
736 AliHFEsecVtxs hfesecvtx;
737 hfesecvtx.SetTrkLabel1(htracklabel1);
738 hfesecvtx.SetTrkLabel2(htracklabel2);
739 if(HasMCData()) hfesecvtx.SetMCCode(GetElectronSource(TMath::Abs(track->GetLabel())));
740 hfesecvtx.SetKFChi2(fKFchi2);
741 hfesecvtx.SetInvmass(fInvmass);
742 hfesecvtx.SetSignedLxy(fSignedLxy);
743 hfesecvtx.SetSignedLxy2(fSignedLxy2);
744 hfesecvtx.SetKFIP(fKFip);
745 hfesecvtx.SetKFIP2(fKFip2);
746 AddHFEsecvtxToArray(&hfesecvtx);
747 fNoOfHFEsecvtxs++;
748
749 // fill debugging THnSparse
750 dataE[0]=fInvmass;
751 dataE[1]=fKFchi2;
752 dataE[2]=fSignedLxy;
753 dataE[3]=fKFip;
754 if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
755 dataE[5]=3;
756 if(paircode1 & paircode2) dataE[6]=1;
757 else if(!paircode1 & !paircode2) dataE[6]=0;
758 else dataE[6]=3;
759 dataE[7]=fSignedLxy2;
760 dataE[8]=track->Pt();
761 fSecvtxQA->Fill(dataE);
762
763}
764
765//_______________________________________________________________________________________________
766void AliHFEsecVtx::Fill2TrkSECVTX(AliVTrack* track, const AliHFEpairs *pair)
767{
768 //
769 // fill 2 tracks' secondary vertex properties
770 //
771
772 Double_t dataE[9]={-999.,-999.,-999.,-999.,-1.,0,0,-999.,-999.};
773
774 Int_t paircode;
775 if (pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode=1;
776 else paircode=0;
777
778 // fill secondary vertex container
779 AliHFEsecVtxs hfesecvtx;
780 hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
781 hfesecvtx.SetTrkLabel2(-999);
782 if(HasMCData()) hfesecvtx.SetMCCode(GetElectronSource(TMath::Abs(track->GetLabel())));
783 hfesecvtx.SetInvmass(pair->GetInvmass());
784 hfesecvtx.SetKFChi2(pair->GetKFChi2());
785 hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
786 hfesecvtx.SetSignedLxy2(pair->GetSignedLxy2());
787 hfesecvtx.SetKFIP(pair->GetKFIP());
788 hfesecvtx.SetKFIP2(pair->GetKFIP2());
789 AddHFEsecvtxToArray(&hfesecvtx);
790 fNoOfHFEsecvtxs++;
791
792 // fill debugging THnSparse
793 dataE[0]=pair->GetInvmass();
794 dataE[1]=pair->GetKFChi2();
795 dataE[2]=pair->GetSignedLxy();
796 dataE[3]=pair->GetKFIP();
797 if (HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
798 dataE[5]=2; //# of associated tracks
799 dataE[6]=paircode;
800 dataE[7]=pair->GetSignedLxy2();
801 dataE[8]=track->Pt();
802 fSecvtxQA->Fill(dataE);
803
804}
805
806//_______________________________________________________________________________________________
807void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
808{
809 //
810 // calculate secondary vertex properties
811 //
812
813 // get KF particle input pid
814 Int_t pdg1 = GetPDG(track1);
815 Int_t pdg2 = GetPDG(track2);
816 Int_t pdg3 = GetPDG(track3);
817
818 if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
819 //printf("out if considered pid range \n");
820 return;
821 }
822
823 // create KF particle of pair
824 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
825 else AliKFParticle::SetField(fESD1->GetMagneticField());
826 AliKFParticle kfTrack[3];
827 kfTrack[0] = AliKFParticle(*track1, pdg1);
828 kfTrack[1] = AliKFParticle(*track2, pdg2);
829 kfTrack[2] = AliKFParticle(*track3, pdg3);
830
831 AliKFParticle kfSecondary(kfTrack[0],kfTrack[1],kfTrack[2]);
832 //AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
833
834 //secondary vertex point from kf particle
835 Double_t kfx = kfSecondary.GetX();
836 Double_t kfy = kfSecondary.GetY();
837 //Double_t kfz = kfSecondary.GetZ();
838
839 //momentum at the decay point from kf particle
840 Double_t kfpx = kfSecondary.GetPx();
841 Double_t kfpy = kfSecondary.GetPy();
842 //Double_t kfpz = kfSecondary.GetPz();
843
844 Double_t dx = kfx-fPVx;
845 Double_t dy = kfy-fPVy;
846
847 // discriminating variables ----------------------------------------------------------
848
849 if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
850
851 // invariant mass of the KF particle
852 kfSecondary.GetMass(fInvmass,fInvmassSigma);
853
854 // DCA from primary to e-h KF particle (impact parameter of KF particle)
855 Double_t vtx[2]={fPVx, fPVy};
856 fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
857
858 if((dx*kfpx+dy*kfpy)>0) fSignedLxy= TMath::Sqrt(dx*dx+dy*dy);
859 if((dx*kfpx+dy*kfpy)<0) fSignedLxy= -1*TMath::Sqrt(dx*dx+dy*dy);
860 //[the other way to think about] - projection of kf vertex vector to the kf momentum direction
861 //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
862 //if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
863
864
865 //recalculating primary vertex after removing secvtx tracks --------------------------
866 Int_t trkid[3];
867 trkid[0] = track1->GetID();
868 trkid[1] = track2->GetID();
869 trkid[2] = track3->GetID();
870
871 RecalcPrimvtx(3, trkid, kfTrack);
872 Double_t dx2 = kfx-fPVx2;
873 Double_t dy2 = kfy-fPVy2;
874
875 // IP of sec particle recalculated based on recalculated primary vertex
876 Double_t vtx2[2]={fPVx2, fPVy2};
877 fKFip2 = kfSecondary.GetDistanceFromVertexXY(vtx2);
878 // signed Lxy recalculated based on recalculated primary vertex
879 if((dx2*kfpx+dy2*kfpy)>0) fSignedLxy2= TMath::Sqrt(dx2*dx2+dy2*dy2);
880 if((dx2*kfpx+dy2*kfpy)<0) fSignedLxy2= -1*TMath::Sqrt(dx2*dx2+dy2*dy2);
881 //------------------------------------------------------------------------------------
882
883}
884
885//_______________________________________________________________________________________________
886void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3, AliVTrack* track4)
887{
888 //
889 // calculate secondary vertex properties
890 //
891
892 // get KF particle input pid
893 Int_t pdg1 = GetPDG(track1);
894 Int_t pdg2 = GetPDG(track2);
895 Int_t pdg3 = GetPDG(track3);
896 Int_t pdg4 = GetPDG(track4);
897
898 if(pdg1==-1 || pdg2==-1 || pdg3==-1 || pdg4==-1) {
899 //printf("out if considered pid range \n");
900 return;
901 }
902
903 // create KF particle of pair
904 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
905 else AliKFParticle::SetField(fESD1->GetMagneticField());
906
907 AliKFParticle kfTrack[4];
908 kfTrack[0] = AliKFParticle(*track1, pdg1);
909 kfTrack[1] = AliKFParticle(*track2, pdg2);
910 kfTrack[2] = AliKFParticle(*track3, pdg3);
911 kfTrack[3] = AliKFParticle(*track4, pdg4);
912
913 AliKFParticle kfSecondary(kfTrack[0],kfTrack[1],kfTrack[2],kfTrack[3]);
914
915 //secondary vertex point from kf particle
916 Double_t kfx = kfSecondary.GetX();
917 Double_t kfy = kfSecondary.GetY();
918 //Double_t kfz = kfSecondary.GetZ();
919
920 //momentum at the decay point from kf particle
921 Double_t kfpx = kfSecondary.GetPx();
922 Double_t kfpy = kfSecondary.GetPy();
923 //Double_t kfpz = kfSecondary.GetPz();
924
925 Double_t dx = kfx-fPVx;
926 Double_t dy = kfy-fPVy;
927
928 // discriminating variables ----------------------------------------------------------
929
930 if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
931
932 // invariant mass of the KF particle
933 kfSecondary.GetMass(fInvmass,fInvmassSigma);
934
935 // DCA from primary to e-h KF particle (impact parameter of KF particle)
936 Double_t vtx[2]={fPVx, fPVy};
937 fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
938
939 if((dx*kfpx+dy*kfpy)>0) fSignedLxy= TMath::Sqrt(dx*dx+dy*dy);
940 if((dx*kfpx+dy*kfpy)<0) fSignedLxy= -1*TMath::Sqrt(dx*dx+dy*dy);
941 //[the other way to think about] - projection of kf vertex vector to the kf momentum direction
942 //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
943 //if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
944
945 //recalculating primary vertex after removing secvtx tracks --------------------------
946 Int_t trkid[4];
947 trkid[0] = track1->GetID();
948 trkid[1] = track2->GetID();
949 trkid[2] = track3->GetID();
950 trkid[3] = track4->GetID();
951
952 RecalcPrimvtx(4, trkid, kfTrack);
953 Double_t dx2 = kfx-fPVx2;
954 Double_t dy2 = kfy-fPVy2;
955
956 // IP of sec particle recalculated based on recalculated primary vertex
957 Double_t vtx2[2]={fPVx2, fPVy2};
958 fKFip2 = kfSecondary.GetDistanceFromVertexXY(vtx2);
959 // signed Lxy recalculated based on recalculated primary vertex
960 if((dx2*kfpx+dy2*kfpy)>0) fSignedLxy2= TMath::Sqrt(dx2*dx2+dy2*dy2);
961 if((dx2*kfpx+dy2*kfpy)<0) fSignedLxy2= -1*TMath::Sqrt(dx2*dx2+dy2*dy2);
962 //------------------------------------------------------------------------------------
963
964}
965
966//_______________________________________________________________________________________________
967void AliHFEsecVtx::RecalcPrimvtx(Int_t nkftrk, const Int_t * const trkid, const AliKFParticle * const kftrk){
968 //
969 // reccalculate primary vertex after removing considering track in the calculation
970 //
971
972 const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
973
974 AliKFVertex kfESDprimary;
975 Int_t n = primvtx->GetNIndices();
976 fNsectrk2prim = 0;
977 fPVx2 = -999.;
978 fPVy2 = -999.;
979
980 if (n>0 && primvtx->GetStatus()){
981 kfESDprimary = AliKFVertex(*primvtx);
982 UShort_t *priIndex = primvtx->GetIndices();
983 for(Int_t j=0; j<nkftrk; j++){
984 for (Int_t i=0;i<n;i++){
985 Int_t idx = Int_t(priIndex[i]);
986 if (idx == trkid[j]){
987 kfESDprimary -= kftrk[j];
988 fNsectrk2prim++;
989 }
990 }
991 }
992 }
993
994 fPVx2 = kfESDprimary.GetX();
995 fPVy2 = kfESDprimary.GetY();
996
997}
998
999//_______________________________________________________________________________________________
1000Int_t AliHFEsecVtx::GetMCPID(const AliESDtrack *track)
1001{
1002 //
1003 // return mc pid
1004 //
1005
1006 AliMCParticle *mctrack = NULL;
1007 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return 0;
1008 TParticle *mcpart = mctrack->Particle();
1009
1010 if ( !mcpart ) return 0;
1011 Int_t pdgCode = mcpart->GetPdgCode();
1012
1013 return pdgCode;
1014}
1015
1016//_______________________________________________________________________________________________
1017Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
1018{
1019 //
1020 // return pdg code of the origin(source) of the pair
1021 //
1022 //
1023 // ---*---*---*-----ancester A----- track1
1024 // |____*______
1025 // |______ track2
1026 // => if they originated from same ancester,
1027 // then return "the absolute value of pdg code of ancester A"
1028 //
1029 // ---*---*---B hadron-----ancester A----- track1
1030 // |____*______
1031 // |______ track2
1032 // => if they originated from same ancester, and this ancester originally comes from B hadrons
1033 // then return -1*"the absolute value of pdg code of ancester A"
1034 //
1035 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
1036 //
1037
1038 if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
1039
1040 AliMCParticle *mctrack = NULL;
1041 AliMCParticle *mctrack1 = NULL;
1042 AliMCParticle *mctrack2 = NULL;
1043 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(trk1->GetLabel()))))) return 0;
1044 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(trk2->GetLabel()))))) return 0;
1045 TParticle *part1 = mctrack1->Particle();
1046 TParticle *part2 = mctrack2->Particle();
1047
1048 TParticle* part2cp = part2;
1049 if (!(part1) || !(part2)) return 0;
1050
1051 Int_t srcpdg = 0;
1052
1053 //if the two tracks' mother's label is same, get the mother info
1054 //in case of charm, check if it originated from beauty
1055 for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
1056 Int_t label1 = part1->GetFirstMother();
1057 if (label1 < 0) return 0;
1058
1059 for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
1060 Int_t label2 = part2->GetFirstMother();
1061 if (label2 < 0) break;
1062
1063 if (label1 == label2){ //check if two tracks are originated from same mother
1064 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label2))))) return 0;
1065 TParticle* commonmom = mctrack->Particle();
1066
1067 srcpdg = TMath::Abs(commonmom->GetPdgCode());
1068
1069 //check ancester to see if it is originally from beauty
1070 for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
1071 Int_t ancesterlabel = commonmom->GetFirstMother();
1072 if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
1073
1074 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(ancesterlabel))))) return 0;
1075 TParticle* commonancester = mctrack->Particle();
1076
1077 Int_t ancesterpdg = TMath::Abs(commonancester->GetPdgCode());
1078
1079 for (Int_t l=0; l<fNparents; l++){
1080 if (TMath::Abs(ancesterpdg)==fParentSelect[1][l]){
1081 srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
1082 return srcpdg;
1083 }
1084 }
1085 commonmom = commonancester;
1086 }
1087 }
1088 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label2))))) return 0;
1089 part2 = mctrack->Particle(); //if their mother is different, go to earlier generation of 2nd particle
1090
1091 if (!(part2)) break;
1092 }
1093 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label1))))) return 0;
1094 part1 = mctrack->Particle(); //if their mother is different, go to earlier generation of 1st particle
1095 part2 = part2cp;
1096 if (!(part1)) return 0;
1097 }
1098
1099 return srcpdg;
1100}
1101
1102//_______________________________________________________________________________________________
1103Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
1104{
1105
1106 //
1107 // return pdg code of the origin(source) of the pair
1108 //
1109 //
1110 // ---*---*---*-----ancester A----- track1
1111 // |____*______
1112 // |______ track2
1113 // => if they originated from same ancester,
1114 // then return "the absolute value of pdg code of ancester A"
1115 //
1116 // ---*---*---B hadron-----ancester A----- track1
1117 // |____*______
1118 // |______ track2
1119 // => if they originated from same ancester, and this ancester originally comes from B hadrons
1120 // then return -1*"the absolute value of pdg code of ancester A"
1121 //
1122 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
1123 //
1124
1125 if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
1126 AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
1127 AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
1128 AliAODMCParticle *part2cp = part2;
1129 if (!(part1) || !(part2)) return 0;
1130
1131 Int_t srcpdg = 0;
1132
1133 //if the two tracks' mother's label is same, get the mother info
1134 //in case of charm, check if it originated from beauty
1135 for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
1136 Int_t label1 = part1->GetMother();
1137 if (label1 < 0) return 0;
1138
1139 for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
1140 Int_t label2 = part2->GetMother();
1141 if (label2 < 0) return 0;
1142
1143 if (label1 == label2){ //check if two tracks are originated from same mother
1144 AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
1145 srcpdg = TMath::Abs(commonmom->GetPdgCode());
1146
1147 //check ancester to see if it is originally from beauty
1148 for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
1149 Int_t ancesterlabel = commonmom->GetMother();
1150 if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
1151
1152 AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
1153 Int_t ancesterpdg = TMath::Abs(commonancester->GetPdgCode());
1154
1155 for (Int_t l=0; l<fNparents; l++){
1156 if (TMath::Abs(ancesterpdg)==fParentSelect[1][l]){
1157 srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
1158 return srcpdg;
1159 }
1160 }
1161 commonmom = commonancester;
1162 }
1163 }
1164 part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
1165 if (!(part2)) break;
1166 }
1167 part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
1168 part2 = part2cp;
1169 if (!(part1)) return 0;
1170 }
1171
1172 return srcpdg;
1173}
1174
1175//_______________________________________________________________________________________________
1176Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
1177{
1178 //
1179 // return pair code which is predefinded as:
1180 // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
1181 // kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
1182 //
1183
1184 Int_t srcpdg = -1;
1185 Int_t srccode = kElse;
1186
1187 if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
1188 else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
1189
1190 if (srcpdg < 0) srccode = kBeautyElse;
1191 for (Int_t i=0; i<fNparents; i++){
1192 if (TMath::Abs(srcpdg)==fParentSelect[0][i]) {
1193 if (srcpdg>0) srccode = kDirectCharm;
1194 if (srcpdg<0) srccode = kBeautyCharm;
1195 }
1196 if (TMath::Abs(srcpdg)==fParentSelect[1][i]) {
1197 if (srcpdg>0) srccode = kDirectBeauty;
1198 if (srcpdg<0) return kElse;
1199 }
1200 }
1201 if (srcpdg == 22) srccode = kGamma;
1202 if (srcpdg == -22) srccode = kBeautyGamma;
1203 if (srcpdg == 111) srccode = kPi0;
1204 if (srcpdg == -111) srccode = kBeautyPi0;
1205
1206 return srccode;
1207}
1208
1209//_______________________________________________________________________________________________
1210Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
1211{
1212 //
1213 // return decay electron's origin
1214 //
1215
1216 if (iTrack < 0) {
1217 AliDebug(1, "Stack label is negative, return\n");
1218 return -1;
1219 }
1220
1221 AliMCParticle *mctrack = NULL;
1222 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iTrack))))) return -1;
1223 TParticle *mcpart = mctrack->Particle();
1224
1225 if(!mcpart){
1226 AliDebug(1, "no mc particle, return\n");
1227 return -1;
1228 }
1229
1230 if ( TMath::Abs(mcpart->GetPdgCode()) != 11 ) return kMisID;
1231
1232// if ( TMath::Abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
1233
1234 Int_t iLabel = mcpart->GetFirstMother();
1235 if (iLabel<0){
1236 AliDebug(1, "Stack label is negative, return\n");
1237 return -1;
1238 }
1239
1240 Int_t origin = -1;
1241 Bool_t isFinalOpenCharm = kFALSE;
1242
1243 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1244 TParticle *partMother = mctrack->Particle();
1245
1246 Int_t maPdgcode = partMother->GetPdgCode();
1247
1248 // if the mother is charmed hadron
1249 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1250
1251 for (Int_t i=0; i<fNparents; i++){
1252 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1253 isFinalOpenCharm = kTRUE;
1254 }
1255 }
1256 if (!isFinalOpenCharm) return -1;
1257
1258 // iterate until you find B hadron as a mother or become top ancester
1259 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1260
1261 Int_t jLabel = partMother->GetFirstMother();
1262 if (jLabel == -1){
1263 origin = kDirectCharm;
1264 return origin;
1265 }
1266 if (jLabel < 0){ // safety protection
1267 AliDebug(1, "Stack label is negative, return\n");
1268 return -1;
1269 }
1270
1271 // if there is an ancester
1272 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1273 TParticle *grandMa = mctrack->Particle();
1274
1275 Int_t grandMaPDG = grandMa->GetPdgCode();
1276
1277 for (Int_t j=0; j<fNparents; j++){
1278 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1279 origin = kBeautyCharm;
1280 return origin;
1281 }
1282 }
1283
1284 partMother = grandMa;
1285 } // end of iteration
1286 } // end of if
1287 else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
1288 for (Int_t i=0; i<fNparents; i++){
1289 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1290 origin = kDirectBeauty;
1291 return origin;
1292 }
1293 }
1294 } // end of if
1295
1296 //============ gamma ================
1297 else if ( TMath::Abs(maPdgcode) == 22 ) {
1298 origin = kGamma;
1299
1300 // iterate until you find B hadron as a mother or become top ancester
1301 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1302
1303 Int_t jLabel = partMother->GetFirstMother();
1304 if (jLabel == -1){
1305 origin = kGamma;
1306 return origin;
1307 }
1308 if (jLabel < 0){ // safety protection
1309 AliDebug(1, "Stack label is negative, return\n");
1310 return -1;
1311 }
1312
1313 // if there is an ancester
1314 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1315 TParticle *grandMa = mctrack->Particle();
1316
1317 Int_t grandMaPDG = grandMa->GetPdgCode();
1318
1319 for (Int_t j=0; j<fNparents; j++){
1320 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1321 origin = kBeautyGamma;
1322 return origin;
1323 }
1324 }
1325
1326 partMother = grandMa;
1327 } // end of iteration
1328
1329 return origin;
1330 } // end of if
1331
1332 //============ pi0 ================
1333 else if ( TMath::Abs(maPdgcode) == 111 ) {
1334 origin = kPi0;
1335
1336 // iterate until you find B hadron as a mother or become top ancester
1337 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1338
1339 Int_t jLabel = partMother->GetFirstMother();
1340 if (jLabel == -1){
1341 origin = kPi0;
1342 return origin;
1343 }
1344 if (jLabel < 0){ // safety protection
1345 AliDebug(1, "Stack label is negative, return\n");
1346 return -1;
1347 }
1348
1349 // if there is an ancester
1350 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1351 TParticle *grandMa = mctrack->Particle();
1352
1353 Int_t grandMaPDG = grandMa->GetPdgCode();
1354
1355 for (Int_t j=0; j<fNparents; j++){
1356 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1357 origin = kBeautyPi0;
1358 return origin;
1359 }
1360 }
1361
1362 partMother = grandMa;
1363 } // end of iteration
1364
1365 return origin;
1366 } // end of if
1367
1368 else {
1369 origin = kElse;
1370 // iterate until you find B hadron as a mother or become top ancester
1371 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1372
1373 Int_t jLabel = partMother->GetFirstMother();
1374 if (jLabel == -1){
1375 origin = kElse;
1376 return origin;
1377 }
1378 if (jLabel < 0){ // safety protection
1379 AliDebug(1, "Stack label is negative, return\n");
1380 return -1;
1381 }
1382
1383 // if there is an ancester
1384 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1385 TParticle *grandMa = mctrack->Particle();
1386
1387 Int_t grandMaPDG = grandMa->GetPdgCode();
1388
1389 for (Int_t j=0; j<fNparents; j++){
1390 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1391 origin = kBeautyElse;
1392 return origin;
1393 }
1394 }
1395
1396 partMother = grandMa;
1397 } // end of iteration
1398 }
1399
1400 return origin;
1401}
1402
1403//_______________________________________________________________________________________________
1404Int_t AliHFEsecVtx::GetPDG(const AliVTrack *track)
1405{
1406 //
1407 // get KF particle input pdg for mass hypothesis
1408 //
1409
1410 Int_t pdgCode=-1;
1411
1412 if (fUseMCPID && HasMCData()){
1413 pdgCode = GetMCPDG(track);
1414 }
1415 else if(fESD1){
1416 Int_t pid=0;
1417 Double_t prob;
1418 GetESDPID((AliESDtrack*)track, pid, prob);
1419 switch(pid){
1420 case 0: pdgCode = 11; break;
1421 case 1: pdgCode = 13; break;
1422 case 2: pdgCode = 211; break;
1423 case 3: pdgCode = 321; break;
1424 case 4: pdgCode = 2212; break;
1425 default: pdgCode = -1;
1426 }
1427 }
1428 else if(fAOD1){
1429 Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
1430 switch(pid){
1431 case 0: pdgCode = 11; break;
1432 case 1: pdgCode = 13; break;
1433 case 2: pdgCode = 211; break;
1434 case 3: pdgCode = 321; break;
1435 case 4: pdgCode = 2212; break;
1436 default: pdgCode = -1;
1437 }
1438 }
1439
1440 return pdgCode;
1441}
1442
1443//_______________________________________________________________________________________________
1444Int_t AliHFEsecVtx::GetMCPDG(const AliVTrack *track)
1445{
1446 //
1447 // return mc pdg code
1448 //
1449
1450 Int_t label = TMath::Abs(track->GetLabel());
1451 Int_t pdgCode;
1452 AliMCParticle *mctrack = NULL;
1453
1454 if (IsAODanalysis()) {
1455 AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
1456 if ( !mcpart ) return 0;
1457 pdgCode = mcpart->GetPdgCode();
1458 }
1459 else {
1460 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label))))) return 0;
1461 TParticle *mcpart = mctrack->Particle();
1462
1463 if ( !mcpart ) return 0;
1464 pdgCode = mcpart->GetPdgCode();
1465 }
1466
1467 return pdgCode;
1468}
1469
1470//______________________________________________________________________________
1471void AliHFEsecVtx::GetESDPID(const AliESDtrack *track, Int_t &recpid, Double_t &recprob)
1472{
1473 //
1474 // calculate likehood for esd pid
1475 //
1476
1477 recpid = -1;
1478 recprob = -1;
1479
1480 Int_t ipid=-1;
1481 Double_t max=0.;
1482
1483 Double_t probability[5];
1484
1485 // get probability of the diffenrent particle types
1486 track->GetESDpid(probability);
1487
1488 // find most probable particle in ESD pid
1489 // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
1490 ipid = TMath::LocMax(5,probability);
1491 max = TMath::MaxElement(5,probability);
1492
1493 recpid = ipid;
1494 recprob = max;
1495}
1496
1497//_____________________________________________________________________________
1498void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
1499{
1500 //
1501 // Add a HFE pair to the array
1502 //
1503
1504 Int_t n = HFEpairs()->GetEntriesFast();
1505 if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
1506 new((*HFEpairs())[n]) AliHFEpairs(*pair);
1507}
1508
1509//_____________________________________________________________________________
1510TClonesArray *AliHFEsecVtx::HFEpairs()
1511{
1512 //
1513 // Returns the list of HFE pairs
1514 //
1515
1516 if (!fHFEpairs) {
1517 fHFEpairs = new TClonesArray("AliHFEpairs", 200);
1518 }
1519 return fHFEpairs;
1520}
1521
1522//_____________________________________________________________________________
1523void AliHFEsecVtx::DeleteHFEpairs()
1524{
1525 //
1526 // Delete the list of HFE pairs
1527 //
1528
1529 if (fHFEpairs) {
1530 fHFEpairs->Delete();
1531 //delete fHFEpairs;
1532 }
1533}
1534
1535//_____________________________________________________________________________
1536void AliHFEsecVtx::InitHFEpairs()
1537{
1538 //
1539 // Initialization should be done before make all possible pairs for a given electron candidate
1540 //
1541
1542 fNoOfHFEpairs = 0;
1543}
1544
1545//_____________________________________________________________________________
1546void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
1547{
1548 //
1549 // Add a HFE secondary vertex to the array
1550 //
1551
1552 Int_t n = HFEsecvtxs()->GetEntriesFast();
1553 if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
1554 new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
1555}
1556
1557//_____________________________________________________________________________
1558TClonesArray *AliHFEsecVtx::HFEsecvtxs()
1559{
1560 //
1561 // Returns the list of HFE secvtx
1562 //
1563
1564 if (!fHFEsecvtxs) {
1565 fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
1566 }
1567 return fHFEsecvtxs;
1568}
1569
1570//_____________________________________________________________________________
1571void AliHFEsecVtx::DeleteHFEsecvtxs()
1572{
1573 //
1574 // Delete the list of HFE pairs
1575 //
1576
1577 if (fHFEsecvtxs) {
1578 fHFEsecvtxs->Delete();
1579 //delete fHFEsecvtx;
1580 }
1581}
1582
1583//_____________________________________________________________________________
1584void AliHFEsecVtx::InitHFEsecvtxs()
1585{
1586 //
1587 // Initialization should be done
1588 //
1589
1590 fNoOfHFEsecvtxs = 0;
1591}
1592
1593//____________________________________________________________
1594void AliHFEsecVtx::MakeContainer(){
1595
1596 //
1597 // make container
1598 //
1599
1600 const Int_t nDimPair=6;
1601 Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 13};
1602 //Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 13, 60, 60, 2000, 2000};
1603 const Double_t kInvmassmin = 0., kInvmassmax = 20.;
1604 const Double_t kKFChi2min = 0, kKFChi2max= 50;
1605 const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
1606 const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
1607 const Double_t kKFIPmin = -10, kKFIPmax= 10;
1608 const Double_t kPairCodemin = -1, kPairCodemax= 12;
1609 //const Double_t kPtmin = 0, kPtmax= 30;
1610 //const Double_t kDCAsigmin = -5, kDCAsigmax= 5;
1611
1612 Double_t* binEdgesPair[nDimPair];
1613 for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1614 binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
1615
1616 for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
1617 for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
1618 for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
1619 for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
1620 for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
1621 for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
1622 //for(Int_t i=0; i<=nBinPair[6]; i++) binEdgesPair[6][i]=(Double_t)kPtmin + (kPtmax - kPtmin)/nBinPair[6]*(Double_t)i;
1623 //for(Int_t i=0; i<=nBinPair[7]; i++) binEdgesPair[7][i]=binEdgesPair[6][i];
1624 //for(Int_t i=0; i<=nBinPair[6]; i++) binEdgesPair[6][i]=(Double_t)kDCAsigmin + (kDCAsigmax - kDCAsigmin)/nBinPair[6]*(Double_t)i;
1625 //for(Int_t i=0; i<=nBinPair[7]; i++) binEdgesPair[7][i]=binEdgesPair[6][i];
1626
1627 fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code; dca1; dca2", nDimPair, nBinPair);
1628 //fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code; pt1; pt2; dca1; dca2", nDimPair, nBinPair);
1629 for(Int_t idim = 0; idim < nDimPair; idim++){
1630 fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1631 }
1632
1633 fSecVtxList->AddAt(fPairQA,0);
1634
1635 const Int_t nDimSecvtx=9;
1636 Double_t* binEdgesSecvtx[nDimSecvtx];
1637 Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 13, 10, 4, 2000, 500};
1638 const Double_t kNtrksmin = 0, kNtrksmax= 10;
1639 const Double_t kTrueBmin = 0, kTrueBmax= 4;
1640 const Double_t kPtmin = 0, kPtmax= 50;
1641 for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1642 binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
1643
1644 for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
1645 for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
1646 for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
1647 for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
1648 for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
1649 for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
1650 for(Int_t i=0; i<=nBinSecvtx[6]; i++) binEdgesSecvtx[6][i]=(Double_t)kTrueBmin + (kTrueBmax - kTrueBmin)/nBinSecvtx[6]*(Double_t)i;
1651 for(Int_t i=0; i<=nBinSecvtx[7]; i++) binEdgesSecvtx[7][i]=binEdgesPair[3][i];
1652 for(Int_t i=0; i<=nBinSecvtx[8]; i++) binEdgesSecvtx[8][i]=(Double_t)kPtmin + (kPtmax - kPtmin)/nBinSecvtx[8]*(Double_t)i;
1653
1654 fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
1655 for(Int_t idim = 0; idim < nDimSecvtx; idim++){
1656 fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
1657 }
1658
1659 fSecVtxList->AddAt(fSecvtxQA,1);
1660
1661 for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1662 delete[] binEdgesPair[ivar];
1663 for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1664 delete[] binEdgesSecvtx[ivar];
1665}
1666
1667//____________________________________________________________
1668void AliHFEsecVtx::MakeHistos(Int_t step){
1669
1670 //
1671 // make container
1672 //
1673
1674 TString hname=Form("step%d",step);
1675 step = step*7;
1676
1677 const Double_t kPtbound[2] = {0.1, 20.};
1678 Int_t iBin[1];
1679 iBin[0] = 44; // bins in pt
1680 Double_t* binEdges[1];
1681 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
1682
1683 fSecVtxList->AddAt(new TH1F(hname+"taggedElec", "pT of e", iBin[0],binEdges[0]), step);
1684 fSecVtxList->AddAt(new TH1F(hname+"charmElec", "pT of e", iBin[0],binEdges[0]), step+1);
1685 fSecVtxList->AddAt(new TH1F(hname+"beautyElec", "pT of e", iBin[0],binEdges[0]), step+2);
1686 fSecVtxList->AddAt(new TH1F(hname+"conversionElec", "pT of e", iBin[0],binEdges[0]), step+3);
1687 fSecVtxList->AddAt(new TH1F(hname+"ebgElec", "pT of e", iBin[0],binEdges[0]), step+4);
1688 fSecVtxList->AddAt(new TH1F(hname+"hcontaminElec", "pT of e", iBin[0],binEdges[0]), step+5);
1689 fSecVtxList->AddAt(new TH1F(hname+"elseElec", "pT of e", iBin[0],binEdges[0]), step+6);
1690 delete[] binEdges[0];
1691}
1692
1693//____________________________________________________________
1694void AliHFEsecVtx::FillHistos(Int_t step, const AliESDtrack *track){
1695
1696 //
1697 // make container
1698 //
1699
1700 step = step*7;
1701
1702 AliMCParticle *mctrack = NULL;
1703 TParticle* mcpart = NULL;
1704
1705 (static_cast<TH1F *>(fSecVtxList->At(step)))->Fill(track->Pt()); // electrons tagged
1706
1707 if(HasMCData() && fMCQA){
1708 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return;
1709 mcpart = mctrack->Particle();
1710
1711 Int_t esource=fMCQA->GetElecSource(mcpart,kTRUE);
1712 if(esource==1) {
1713 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+1)))) return;
1714 (static_cast<TH1F *>(fSecVtxList->At(step+1)))->Fill(mcpart->Pt()); //charm
1715 }
1716 else if(esource==2 || esource==3) {
1717 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+2)))) return;
1718 (static_cast<TH1F *>(fSecVtxList->At(step+2)))->Fill(mcpart->Pt()); //beauty
1719 }
1720 else if(esource>12 && esource<19) {
1721 //else if(esource==4) {
1722 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+3)))) return;
1723 (static_cast<TH1F *>(fSecVtxList->At(step+3)))->Fill(mcpart->Pt()); //conversion
1724 }
1725 else if(esource==7) {
1726 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+5)))) return;
1727 (static_cast<TH1F *>(fSecVtxList->At(step+5)))->Fill(mcpart->Pt()); //contamination
1728 }
1729 else if(!(esource<0)) {
1730 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+4)))) return;
1731 (static_cast<TH1F *>(fSecVtxList->At(step+4)))->Fill(mcpart->Pt()); //e backgrounds
1732 }
1733 else {
1734 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+6)))) return;
1735 (static_cast<TH1F *>(fSecVtxList->At(step+6)))->Fill(mcpart->Pt()); //something else?
1736 }
1737 }
1738
1739}