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