]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEsecVtx.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update of the HFE package
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEsecVtx.cxx
CommitLineData
9bcfd1ab 1 /**************************************************************************
259c3296 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 **************************************************************************/
50685501 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//
259c3296 23
259c3296 24#include <TH2F.h>
70da6c5a 25#include <TIterator.h>
259c3296 26#include <TParticle.h>
27
faee3b18 28#include <AliESDVertex.h>
259c3296 29#include <AliESDEvent.h>
9bcfd1ab 30#include <AliAODEvent.h>
31#include <AliVTrack.h>
259c3296 32#include <AliESDtrack.h>
9bcfd1ab 33#include <AliAODTrack.h>
259c3296 34#include "AliHFEsecVtx.h"
35#include <AliKFParticle.h>
36#include <AliKFVertex.h>
37#include <AliLog.h>
faee3b18 38#include <AliMCEvent.h>
9bcfd1ab 39#include <AliAODMCParticle.h>
40#include "AliHFEpairs.h"
41#include "AliHFEsecVtxs.h"
70da6c5a 42#include "AliHFEtrackFilter.h"
3a72645a 43#include "AliHFEmcQA.h"
44#include "AliHFEtools.h"
259c3296 45
46ClassImp(AliHFEsecVtx)
47
48//_______________________________________________________________________________________________
49AliHFEsecVtx::AliHFEsecVtx():
70da6c5a 50 fFilter(0x0)
51 ,fESD1(0x0)
9bcfd1ab 52 ,fAOD1(0x0)
faee3b18 53 ,fMCEvent(0x0)
3a72645a 54 ,fMCQA(0x0)
9bcfd1ab 55 ,fUseMCPID(kFALSE)
56 ,fkSourceLabel()
57 ,fNparents(0)
58 ,fParentSelect()
70da6c5a 59 ,fPtRng()
60 ,fDcaCut()
9bcfd1ab 61 ,fNoOfHFEpairs(0)
62 ,fNoOfHFEsecvtxs(0)
faee3b18 63 ,fArethereSecVtx(0)
9bcfd1ab 64 ,fHFEpairs(0x0)
65 ,fHFEsecvtxs(0x0)
66 ,fMCArray(0x0)
faee3b18 67 ,fPVx(-999)
68 ,fPVy(-999)
69 ,fPVx2(999)
70 ,fPVy2(-999)
9bcfd1ab 71 ,fCosPhi(-1)
72 ,fSignedLxy(-1)
faee3b18 73 ,fSignedLxy2(-1)
9bcfd1ab 74 ,fKFchi2(-1)
75 ,fInvmass(-1)
76 ,fInvmassSigma(-1)
77 ,fKFip(0)
faee3b18 78 ,fKFip2(0)
79 ,fNsectrk2prim(0)
80 ,fVtxchi2Tightcut(3.)
81 ,fVtxchi2Loosecut(5.)
9bcfd1ab 82 ,fPairQA(0x0)
83 ,fSecvtxQA(0x0)
84 ,fSecVtxList(0x0)
259c3296 85{
9bcfd1ab 86 //
87 // Default constructor
88 //
259c3296 89
9bcfd1ab 90 Init();
259c3296 91}
92
93//_______________________________________________________________________________________________
94AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
9bcfd1ab 95 TObject(p)
70da6c5a 96 ,fFilter(0x0)
9bcfd1ab 97 ,fESD1(0x0)
98 ,fAOD1(0x0)
faee3b18 99 ,fMCEvent(0x0)
3a72645a 100 ,fMCQA(0x0)
9bcfd1ab 101 ,fUseMCPID(p.fUseMCPID)
102 ,fkSourceLabel()
103 ,fNparents(p.fNparents)
104 ,fParentSelect()
70da6c5a 105 ,fPtRng()
106 ,fDcaCut()
9bcfd1ab 107 ,fNoOfHFEpairs(p.fNoOfHFEpairs)
108 ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
faee3b18 109 ,fArethereSecVtx(p.fArethereSecVtx)
9bcfd1ab 110 ,fHFEpairs(0x0)
111 ,fHFEsecvtxs(0x0)
112 ,fMCArray(0x0)
113 ,fPVx(p.fPVx)
114 ,fPVy(p.fPVy)
faee3b18 115 ,fPVx2(p.fPVx2)
116 ,fPVy2(p.fPVy2)
9bcfd1ab 117 ,fCosPhi(p.fCosPhi)
118 ,fSignedLxy(p.fSignedLxy)
faee3b18 119 ,fSignedLxy2(p.fSignedLxy2)
9bcfd1ab 120 ,fKFchi2(p.fKFchi2)
121 ,fInvmass(p.fInvmass)
122 ,fInvmassSigma(p.fInvmassSigma)
123 ,fKFip(p.fKFip)
faee3b18 124 ,fKFip2(p.fKFip2)
125 ,fNsectrk2prim(p.fNsectrk2prim)
126 ,fVtxchi2Tightcut(p.fVtxchi2Tightcut)
127 ,fVtxchi2Loosecut(p.fVtxchi2Loosecut)
9bcfd1ab 128 ,fPairQA(0x0)
129 ,fSecvtxQA(0x0)
130 ,fSecVtxList(0x0)
259c3296 131{
9bcfd1ab 132 //
133 // Copy constructor
134 //
70da6c5a 135 fFilter = new AliHFEtrackFilter(*p.fFilter);
259c3296 136}
137
138//_______________________________________________________________________________________________
139AliHFEsecVtx&
140AliHFEsecVtx::operator=(const AliHFEsecVtx &)
141{
9bcfd1ab 142 //
259c3296 143 // Assignment operator
9bcfd1ab 144 //
259c3296 145
146 AliInfo("Not yet implemented.");
147 return *this;
148}
149
150//_______________________________________________________________________________________________
151AliHFEsecVtx::~AliHFEsecVtx()
152{
9bcfd1ab 153 //
154 // Destructor
155 //
259c3296 156
9bcfd1ab 157 //cout << "Analysis Done." << endl;
70da6c5a 158 delete fFilter;
259c3296 159}
160
161//__________________________________________
162void AliHFEsecVtx::Init()
163{
9bcfd1ab 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;
70da6c5a 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();
259c3296 206}
207
c2690925 208Bool_t AliHFEsecVtx::Process(AliVTrack *signalTrack){
faee3b18 209 //
210 // Run Process
211 //
3a72645a 212 //if(signalTrack->Pt() < 1.0) return;
213
214
70da6c5a 215 AliESDtrack *track = dynamic_cast<AliESDtrack *>(signalTrack);
3a72645a 216
bf892a6a 217 if(!track){
218 AliDebug(1, "no esd track pointer, return\n");
c2690925 219 return kFALSE;
bf892a6a 220 }
221
c2690925 222 Bool_t bTagged = kFALSE;
223
3a72645a 224 FillHistos(0,track); // wo any cuts
225
70da6c5a 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;
3a72645a 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 }
70da6c5a 247 HFEpairs()->Compress();
3a72645a 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
70da6c5a 251 for(int ip=0; ip<HFEsecvtxs()->GetEntriesFast(); ip++){
252 AliHFEsecVtxs *secvtx=0x0;
253 secvtx = (AliHFEsecVtxs*) (HFEsecvtxs()->UncheckedAt(ip));
3a72645a 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);
70da6c5a 257 }
3a72645a 258
259 // fill histos for raw spectra
c2690925 260 if(HFEsecvtxs()->GetEntriesFast()) {
261 FillHistos(3,track); //after secvtx cut
262 bTagged = kTRUE;
263 }
3a72645a 264
70da6c5a 265 DeleteHFEpairs();
266 DeleteHFEsecvtxs();
c2690925 267
268 return bTagged;
70da6c5a 269}
270
9bcfd1ab 271//_______________________________________________________________________________________________
e3fc062d 272void AliHFEsecVtx::CreateHistograms(TList * const qaList)
9bcfd1ab 273{
274 //
275 // create histograms
276 //
e3fc062d 277
278 if(!qaList) return;
9bcfd1ab 279
e3fc062d 280 fSecVtxList = qaList;
9bcfd1ab 281 fSecVtxList->SetName("SecVtx");
282
283 MakeContainer();
3a72645a 284 MakeHistos(0);
285 MakeHistos(1);
286 MakeHistos(2);
287 MakeHistos(3);
288
9bcfd1ab 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
e3fc062d 307 //qaList->AddLast(fSecVtxList);
259c3296 308}
309
9bcfd1ab 310//_______________________________________________________________________________________________
311void AliHFEsecVtx::GetPrimaryCondition()
259c3296 312{
9bcfd1ab 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();
faee3b18 321 fPVy = primVtxCopy.GetY();
9bcfd1ab 322 }
323 else if(fAOD1) {
324 AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
325 if( primVtxCopy.GetNDF() <1 ) return;
326 fPVx = primVtxCopy.GetX();
faee3b18 327 fPVy = primVtxCopy.GetY();
9bcfd1ab 328 }
259c3296 329}
330
331//_______________________________________________________________________________________________
9bcfd1ab 332void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
333{
334 //
335 // calculate e-h pair characteristics and tag pair
336 //
259c3296 337
70da6c5a 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
faee3b18 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
70da6c5a 345 if (IsAODanalysis()){
faee3b18 346 const AliAODVertex *primVtx = fAOD1->GetPrimaryVertex();
70da6c5a 347 AliESDtrack esdTrk1(track1);
348 AliESDtrack esdTrk2(track2);
faee3b18 349 esdTrk1.PropagateToDCA(primVtx,0.,10000.,dca1aod,cov1aod);
350 esdTrk2.PropagateToDCA(primVtx,0.,10000.,dca2aod,cov2aod);
70da6c5a 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
faee3b18 358 /*for(int ibin=0; ibin<6; ibin++){
70da6c5a 359 if((track2->Pt()>fPtRng[ibin] && track2->Pt()<fPtRng[ibin+1]) && TMath::Abs(dca2[0])<fDcaCut[ibin]) return;
faee3b18 360 }*/
70da6c5a 361
9bcfd1ab 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());
9bcfd1ab 374
faee3b18 375 AliKFParticle kfTrack[2];
376 kfTrack[0] = AliKFParticle(*track1, pdg1);
377 kfTrack[1] = AliKFParticle(*track2, pdg2);
378
379 AliKFParticle kfSecondary(kfTrack[0],kfTrack[1]);
9bcfd1ab 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
faee3b18 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
9bcfd1ab 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
faee3b18 421 Double_t phi = kfTrack[0].GetAngleXY(kfTrack[1]);
9bcfd1ab 422 Double_t cosphi = TMath::Cos(phi);
423
faee3b18 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
9bcfd1ab 428 // projection of kf vertex vector to the kf momentum direction
429 Double_t signedLxy=-999.;
70da6c5a 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);
9bcfd1ab 432 //[the other way to think about]
70da6c5a 433 //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
434 //if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
9bcfd1ab 435
faee3b18 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
9bcfd1ab 454
9bcfd1ab 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);
faee3b18 465 hfepair.SetSignedLxy2(signedLxy2);
9bcfd1ab 466 hfepair.SetKFIP(kfip);
faee3b18 467 hfepair.SetKFIP2(kfip2);
9bcfd1ab 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;
faee3b18 477 dataE[3]=signedLxy2;
478 dataE[4]=kfip2;
9bcfd1ab 479 dataE[5]=paircode;
70da6c5a 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]);
faee3b18 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
9bcfd1ab 489 fPairQA->Fill(dataE);
70da6c5a 490
259c3296 491}
492
493//_______________________________________________________________________________________________
9bcfd1ab 494void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
259c3296 495{
9bcfd1ab 496 //
497 // run secondary vertexing algorithm and do tagging
498 //
259c3296 499
9bcfd1ab 500 //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
501 FindSECVTXCandid(track);
502}
259c3296 503
9bcfd1ab 504//_______________________________________________________________________________________________
505void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
506{
507 //
508 // find secondary vertex candidate and store them into container
509 //
510
511 AliVTrack *htrack[20];
faee3b18 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
9bcfd1ab 520 if (HFEpairs()->GetEntriesFast()>20){
521 AliDebug(3, "number of paired hadron is over maximum(20)");
522 return;
523 }
259c3296 524
9bcfd1ab 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));
faee3b18 530 //htracklabel[ip] = pair->GetTrkLabel();
9bcfd1ab 531 htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
faee3b18 532 //if(pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode[ip]=1;
533 //else paircode[ip]=0;
534 //vtxchi2[ip] = pair->GetKFChi2();
9bcfd1ab 535 }
536 }
537 else{
538 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
539 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
faee3b18 540 //htracklabel[ip] = pair->GetTrkLabel();
9bcfd1ab 541 htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
faee3b18 542 //if(pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode[ip]=1;
543 //else paircode[ip]=0;
544 //vtxchi2[ip] = pair->GetKFChi2();
9bcfd1ab 545 }
546 }
faee3b18 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);
9bcfd1ab 555 }
556 return;
557 }
faee3b18 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
9bcfd1ab 582
faee3b18 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 }
259c3296 621 }
faee3b18 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;
9bcfd1ab 679 }
faee3b18 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//_______________________________________________________________________________________________
c2690925 766void AliHFEsecVtx::Fill2TrkSECVTX(AliVTrack* track, const AliHFEpairs *pair)
faee3b18 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
259c3296 804}
805
806//_______________________________________________________________________________________________
9bcfd1ab 807void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
259c3296 808{
9bcfd1ab 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());
faee3b18 826 AliKFParticle kfTrack[3];
827 kfTrack[0] = AliKFParticle(*track1, pdg1);
828 kfTrack[1] = AliKFParticle(*track2, pdg2);
829 kfTrack[2] = AliKFParticle(*track3, pdg3);
9bcfd1ab 830
faee3b18 831 AliKFParticle kfSecondary(kfTrack[0],kfTrack[1],kfTrack[2]);
832 //AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
9bcfd1ab 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();
259c3296 838
9bcfd1ab 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();
259c3296 843
9bcfd1ab 844 Double_t dx = kfx-fPVx;
845 Double_t dy = kfy-fPVy;
dbe3abbe 846
9bcfd1ab 847 // discriminating variables ----------------------------------------------------------
dbe3abbe 848
9bcfd1ab 849 if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
259c3296 850
9bcfd1ab 851 // invariant mass of the KF particle
852 kfSecondary.GetMass(fInvmass,fInvmassSigma);
259c3296 853
9bcfd1ab 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);
259c3296 857
70da6c5a 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);
faee3b18 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//_______________________________________________________________________________________________
3a72645a 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 //
faee3b18 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
259c3296 997}
998
999//_______________________________________________________________________________________________
3a72645a 1000Int_t AliHFEsecVtx::GetMCPID(const AliESDtrack *track)
259c3296 1001{
9bcfd1ab 1002 //
1003 // return mc pid
1004 //
259c3296 1005
faee3b18 1006 AliMCParticle *mctrack = NULL;
1007 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return 0;
1008 TParticle *mcpart = mctrack->Particle();
1009
9bcfd1ab 1010 if ( !mcpart ) return 0;
1011 Int_t pdgCode = mcpart->GetPdgCode();
259c3296 1012
9bcfd1ab 1013 return pdgCode;
259c3296 1014}
1015
1016//_______________________________________________________________________________________________
9bcfd1ab 1017Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
259c3296 1018{
9bcfd1ab 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;
faee3b18 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
9bcfd1ab 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
faee3b18 1064 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label2))))) return 0;
1065 TParticle* commonmom = mctrack->Particle();
1066
9bcfd1ab 1067 srcpdg = 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
faee3b18 1074 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(ancesterlabel))))) return 0;
1075 TParticle* commonancester = mctrack->Particle();
1076
9bcfd1ab 1077 Int_t ancesterpdg = abs(commonancester->GetPdgCode());
1078
1079 for (Int_t l=0; l<fNparents; l++){
1080 if (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 }
faee3b18 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
9bcfd1ab 1091 if (!(part2)) break;
1092 }
faee3b18 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
9bcfd1ab 1095 part2 = part2cp;
1096 if (!(part1)) return 0;
1097 }
1098
1099 return srcpdg;
259c3296 1100}
1101
1102//_______________________________________________________________________________________________
9bcfd1ab 1103Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
259c3296 1104{
1105
9bcfd1ab 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 = 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 = abs(commonancester->GetPdgCode());
1154
1155 for (Int_t l=0; l<fNparents; l++){
1156 if (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;
259c3296 1173}
1174
1175//_______________________________________________________________________________________________
9bcfd1ab 1176Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
259c3296 1177{
9bcfd1ab 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 (abs(srcpdg)==fParentSelect[0][i]) {
1193 if (srcpdg>0) srccode = kDirectCharm;
1194 if (srcpdg<0) srccode = kBeautyCharm;
1195 }
1196 if (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;
259c3296 1205
9bcfd1ab 1206 return srccode;
259c3296 1207}
1208
1209//_______________________________________________________________________________________________
9bcfd1ab 1210Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
259c3296 1211{
9bcfd1ab 1212 //
1213 // return decay electron's origin
1214 //
259c3296 1215
9bcfd1ab 1216 if (iTrack < 0) {
1217 AliDebug(1, "Stack label is negative, return\n");
1218 return -1;
1219 }
259c3296 1220
faee3b18 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 ( abs(mcpart->GetPdgCode()) != 11 ) return kMisID;
41d0c3b3 1231
9bcfd1ab 1232// if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
259c3296 1233
9bcfd1ab 1234 Int_t iLabel = mcpart->GetFirstMother();
1235 if (iLabel<0){
1236 AliDebug(1, "Stack label is negative, return\n");
1237 return -1;
1238 }
259c3296 1239
9bcfd1ab 1240 Int_t origin = -1;
1241 Bool_t isFinalOpenCharm = kFALSE;
259c3296 1242
faee3b18 1243 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1244 TParticle *partMother = mctrack->Particle();
1245
9bcfd1ab 1246 Int_t maPdgcode = partMother->GetPdgCode();
259c3296 1247
9bcfd1ab 1248 // if the mother is charmed hadron
1249 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 1250
9bcfd1ab 1251 for (Int_t i=0; i<fNparents; i++){
1252 if (abs(maPdgcode)==fParentSelect[0][i]){
1253 isFinalOpenCharm = kTRUE;
1254 }
1255 }
1256 if (!isFinalOpenCharm) return -1;
259c3296 1257
9bcfd1ab 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
259c3296 1260
9bcfd1ab 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
faee3b18 1272 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1273 TParticle *grandMa = mctrack->Particle();
1274
9bcfd1ab 1275 Int_t grandMaPDG = grandMa->GetPdgCode();
1276
1277 for (Int_t j=0; j<fNparents; j++){
1278 if (abs(grandMaPDG)==fParentSelect[1][j]){
1279 origin = kBeautyCharm;
1280 return origin;
dbe3abbe 1281 }
9bcfd1ab 1282 }
1283
1284 partMother = grandMa;
1285 } // end of iteration
1286 } // end of if
1287 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1288 for (Int_t i=0; i<fNparents; i++){
1289 if (abs(maPdgcode)==fParentSelect[1][i]){
1290 origin = kDirectBeauty;
1291 return origin;
1292 }
1293 }
1294 } // end of if
dbe3abbe 1295
9bcfd1ab 1296 //============ gamma ================
1297 else if ( abs(maPdgcode) == 22 ) {
1298 origin = kGamma;
259c3296 1299
9bcfd1ab 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
259c3296 1302
9bcfd1ab 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
faee3b18 1314 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1315 TParticle *grandMa = mctrack->Particle();
1316
9bcfd1ab 1317 Int_t grandMaPDG = grandMa->GetPdgCode();
1318
1319 for (Int_t j=0; j<fNparents; j++){
1320 if (abs(grandMaPDG)==fParentSelect[1][j]){
1321 origin = kBeautyGamma;
1322 return origin;
1323 }
1324 }
259c3296 1325
9bcfd1ab 1326 partMother = grandMa;
1327 } // end of iteration
259c3296 1328
9bcfd1ab 1329 return origin;
1330 } // end of if
259c3296 1331
9bcfd1ab 1332 //============ pi0 ================
1333 else if ( abs(maPdgcode) == 111 ) {
1334 origin = kPi0;
259c3296 1335
9bcfd1ab 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
259c3296 1338
9bcfd1ab 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
faee3b18 1350 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1351 TParticle *grandMa = mctrack->Particle();
1352
9bcfd1ab 1353 Int_t grandMaPDG = grandMa->GetPdgCode();
1354
1355 for (Int_t j=0; j<fNparents; j++){
1356 if (abs(grandMaPDG)==fParentSelect[1][j]){
1357 origin = kBeautyPi0;
1358 return origin;
dbe3abbe 1359 }
9bcfd1ab 1360 }
dbe3abbe 1361
9bcfd1ab 1362 partMother = grandMa;
1363 } // end of iteration
259c3296 1364
9bcfd1ab 1365 return origin;
1366 } // end of if
259c3296 1367
9bcfd1ab 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
259c3296 1372
9bcfd1ab 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
faee3b18 1384 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1385 TParticle *grandMa = mctrack->Particle();
1386
9bcfd1ab 1387 Int_t grandMaPDG = grandMa->GetPdgCode();
1388
1389 for (Int_t j=0; j<fNparents; j++){
1390 if (abs(grandMaPDG)==fParentSelect[1][j]){
1391 origin = kBeautyElse;
1392 return origin;
259c3296 1393 }
9bcfd1ab 1394 }
259c3296 1395
9bcfd1ab 1396 partMother = grandMa;
1397 } // end of iteration
1398 }
dbe3abbe 1399
9bcfd1ab 1400 return origin;
259c3296 1401}
1402
1403//_______________________________________________________________________________________________
e156c3bb 1404Int_t AliHFEsecVtx::GetPDG(const AliVTrack *track)
259c3296 1405{
9bcfd1ab 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 }
dbe3abbe 1439
9bcfd1ab 1440 return pdgCode;
259c3296 1441}
1442
1443//_______________________________________________________________________________________________
3a72645a 1444Int_t AliHFEsecVtx::GetMCPDG(const AliVTrack *track)
259c3296 1445{
9bcfd1ab 1446 //
1447 // return mc pdg code
1448 //
1449
1450 Int_t label = TMath::Abs(track->GetLabel());
1451 Int_t pdgCode;
faee3b18 1452 AliMCParticle *mctrack = NULL;
9bcfd1ab 1453
1454 if (IsAODanalysis()) {
1455 AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
1456 if ( !mcpart ) return 0;
1457 pdgCode = mcpart->GetPdgCode();
1458 }
1459 else {
faee3b18 1460 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label))))) return 0;
1461 TParticle *mcpart = mctrack->Particle();
1462
9bcfd1ab 1463 if ( !mcpart ) return 0;
1464 pdgCode = mcpart->GetPdgCode();
1465 }
1466
1467 return pdgCode;
259c3296 1468}
1469
9bcfd1ab 1470//______________________________________________________________________________
c2690925 1471void AliHFEsecVtx::GetESDPID(const AliESDtrack *track, Int_t &recpid, Double_t &recprob)
259c3296 1472{
9bcfd1ab 1473 //
1474 // calculate likehood for esd pid
1475 //
259c3296 1476
9bcfd1ab 1477 recpid = -1;
1478 recprob = -1;
259c3296 1479
9bcfd1ab 1480 Int_t ipid=-1;
1481 Double_t max=0.;
259c3296 1482
9bcfd1ab 1483 Double_t probability[5];
259c3296 1484
9bcfd1ab 1485 // get probability of the diffenrent particle types
1486 track->GetESDpid(probability);
259c3296 1487
9bcfd1ab 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);
259c3296 1492
9bcfd1ab 1493 recpid = ipid;
1494 recprob = max;
259c3296 1495}
1496
9bcfd1ab 1497//_____________________________________________________________________________
1498void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
259c3296 1499{
9bcfd1ab 1500 //
1501 // Add a HFE pair to the array
1502 //
259c3296 1503
9bcfd1ab 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);
259c3296 1507}
1508
9bcfd1ab 1509//_____________________________________________________________________________
1510TClonesArray *AliHFEsecVtx::HFEpairs()
259c3296 1511{
9bcfd1ab 1512 //
1513 // Returns the list of HFE pairs
1514 //
1515
1516 if (!fHFEpairs) {
1517 fHFEpairs = new TClonesArray("AliHFEpairs", 200);
1518 }
1519 return fHFEpairs;
259c3296 1520}
1521
9bcfd1ab 1522//_____________________________________________________________________________
1523void AliHFEsecVtx::DeleteHFEpairs()
259c3296 1524{
9bcfd1ab 1525 //
1526 // Delete the list of HFE pairs
1527 //
1528
1529 if (fHFEpairs) {
1530 fHFEpairs->Delete();
1531 //delete fHFEpairs;
1532 }
259c3296 1533}
1534
9bcfd1ab 1535//_____________________________________________________________________________
1536void AliHFEsecVtx::InitHFEpairs()
dbe3abbe 1537{
9bcfd1ab 1538 //
1539 // Initialization should be done before make all possible pairs for a given electron candidate
1540 //
dbe3abbe 1541
9bcfd1ab 1542 fNoOfHFEpairs = 0;
dbe3abbe 1543}
1544
9bcfd1ab 1545//_____________________________________________________________________________
1546void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
259c3296 1547{
9bcfd1ab 1548 //
1549 // Add a HFE secondary vertex to the array
1550 //
259c3296 1551
9bcfd1ab 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}
259c3296 1556
9bcfd1ab 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}
259c3296 1569
9bcfd1ab 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}
259c3296 1582
9bcfd1ab 1583//_____________________________________________________________________________
1584void AliHFEsecVtx::InitHFEsecvtxs()
1585{
1586 //
1587 // Initialization should be done
1588 //
259c3296 1589
9bcfd1ab 1590 fNoOfHFEsecvtxs = 0;
259c3296 1591}
1592
9bcfd1ab 1593//____________________________________________________________
1594void AliHFEsecVtx::MakeContainer(){
1595
1596 //
1597 // make container
1598 //
1599
1600 const Int_t nDimPair=6;
faee3b18 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};
9bcfd1ab 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;
faee3b18 1608 const Double_t kPairCodemin = -1, kPairCodemax= 12;
70da6c5a 1609 //const Double_t kPtmin = 0, kPtmax= 30;
faee3b18 1610 //const Double_t kDCAsigmin = -5, kDCAsigmax= 5;
9bcfd1ab 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;
faee3b18 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];
9bcfd1ab 1626
faee3b18 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);
9bcfd1ab 1629 for(Int_t idim = 0; idim < nDimPair; idim++){
1630 fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1631 }
1632
1633 fSecVtxList->AddAt(fPairQA,0);
1634
faee3b18 1635 const Int_t nDimSecvtx=9;
9bcfd1ab 1636 Double_t* binEdgesSecvtx[nDimSecvtx];
faee3b18 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;
9bcfd1ab 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;
faee3b18 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;
9bcfd1ab 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);
faee3b18 1660
70da6c5a 1661 for(Int_t ivar = 0; ivar < nDimPair; ivar++)
e156c3bb 1662 delete[] binEdgesPair[ivar];
70da6c5a 1663 for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
e156c3bb 1664 delete[] binEdgesSecvtx[ivar];
259c3296 1665}
3a72645a 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);
e156c3bb 1690 delete[] binEdges[0];
3a72645a 1691}
1692
1693//____________________________________________________________
1694void AliHFEsecVtx::FillHistos(Int_t step, const AliESDtrack *track){
1695
1696 //
1697 // make container
1698 //
bf892a6a 1699
3a72645a 1700 step = step*7;
1701
1702 AliMCParticle *mctrack = NULL;
1703 TParticle* mcpart = NULL;
1704
ccc37cdc 1705 (static_cast<TH1F *>(fSecVtxList->At(step)))->Fill(track->Pt()); // electrons tagged
3a72645a 1706
e3ae862b 1707 if(HasMCData() && fMCQA){
3a72645a 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);
1712 if(esource==1) {
e97c2edf 1713 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+1)))) return;
ccc37cdc 1714 (static_cast<TH1F *>(fSecVtxList->At(step+1)))->Fill(mcpart->Pt()); //charm
3a72645a 1715 }
1716 else if(esource==2 || esource==3) {
e97c2edf 1717 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+2)))) return;
ccc37cdc 1718 (static_cast<TH1F *>(fSecVtxList->At(step+2)))->Fill(mcpart->Pt()); //beauty
3a72645a 1719 }
c2690925 1720 else if(esource>12 && esource<19) {
1721 //else if(esource==4) {
e97c2edf 1722 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+3)))) return;
ccc37cdc 1723 (static_cast<TH1F *>(fSecVtxList->At(step+3)))->Fill(mcpart->Pt()); //conversion
3a72645a 1724 }
1725 else if(esource==7) {
e97c2edf 1726 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+5)))) return;
ccc37cdc 1727 (static_cast<TH1F *>(fSecVtxList->At(step+5)))->Fill(mcpart->Pt()); //contamination
3a72645a 1728 }
1729 else if(!(esource<0)) {
e97c2edf 1730 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+4)))) return;
ccc37cdc 1731 (static_cast<TH1F *>(fSecVtxList->At(step+4)))->Fill(mcpart->Pt()); //e backgrounds
3a72645a 1732 }
1733 else {
e97c2edf 1734 //if(!(dynamic_cast<TH1F *>(fSecVtxList->At(step+6)))) return;
ccc37cdc 1735 (static_cast<TH1F *>(fSecVtxList->At(step+6)))->Fill(mcpart->Pt()); //something else?
3a72645a 1736 }
1737 }
1738
1739}