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