]>
Commit | Line | Data |
---|---|---|
70da6c5a | 1 | /************************************************************************* |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | // | |
16 | // Class for electrons from beauty study | |
17 | // Counting electrons from beauty | |
18 | // by DCA cuts, background subtraction | |
19 | // | |
20 | // Authors: | |
21 | // Hongyan Yang <hongyan@physi.uni-heidelberg.de> | |
faee3b18 | 22 | // Carlo Bombonati <Carlo.Bombonati@cern.ch> |
70da6c5a | 23 | // |
24 | ||
25 | #include "TMath.h" | |
26 | #include "TList.h" | |
27 | #include "AliLog.h" | |
28 | ||
29 | #include <TParticle.h> | |
30 | #include <TDatabasePDG.h> | |
70da6c5a | 31 | #include "THnSparse.h" |
faee3b18 | 32 | |
70da6c5a | 33 | #include "AliMCEvent.h" |
faee3b18 | 34 | #include "AliMCVertex.h" |
70da6c5a | 35 | #include "AliMCParticle.h" |
faee3b18 | 36 | #include "AliStack.h" |
37 | ||
38 | #include "AliESDEvent.h" | |
70da6c5a | 39 | #include "AliESDtrack.h" |
40 | ||
faee3b18 | 41 | #include "AliKFParticle.h" |
42 | #include "AliKFVertex.h" | |
43 | ||
44 | #include "AliVertexerTracks.h" | |
45 | ||
70da6c5a | 46 | |
47 | #include "AliHFEdisplacedElectrons.h" | |
48 | ||
49 | ClassImp(AliHFEdisplacedElectrons) | |
50 | ||
51 | //__________________________________________________________ | |
52 | const Float_t AliHFEdisplacedElectrons::fgkDcaMinPtIntv[13] = { | |
53 | // | |
54 | // define DCA low limits for single electrons cut in each Pt bin | |
55 | // preliminary numbers | |
56 | // these numbers should be replaced by the best numbers determined by | |
57 | // cut efficiency study: signal/background (not used right now) | |
58 | // | |
59 | 450, // 0.0 - 0.5 | |
60 | 450, // 0.5 - 1.0 | |
61 | 450, // 1.0 - 1.5 | |
62 | 500, // 1.5 - 2.0 | |
63 | 400, // 2.0 - 2.5 | |
64 | 300, // 2.5 - 3.0 | |
65 | 300, // 3.0 - 4.0 | |
66 | 300, // 4.0 - 5.0 | |
67 | 200, // 5.0 - 7.0 | |
68 | 150, // 7.0 - 9.0 | |
69 | 150, // 9.0 - 12.0 | |
70 | 100, //12.0 - 16.0 | |
71 | 50}; //16.0 - 20.0 | |
72 | //__________________________________________________________ | |
73 | const Float_t AliHFEdisplacedElectrons::fgkPtIntv[14] = { | |
74 | // | |
75 | // define pT bins for spectra of single electrons | |
76 | // | |
77 | 0.0,0.5,1.0,1.5,2.0,2.5,3.0,4.0,5.0,7.0,9.0,12.0,16.0,20.0}; | |
78 | ||
79 | //__________________________________________________________ | |
80 | const Char_t *AliHFEdisplacedElectrons::fgkKineVar[3] = { | |
81 | "y", "phi", "pt"}; | |
82 | ||
83 | //__________________________________________________________ | |
84 | const Char_t *AliHFEdisplacedElectrons::fgkKineVarTitle[3] ={ | |
85 | "rapidity;y;dN/dy", "azimuthal;#phi;dN/d#phi", "transverse momentum;p_{T};dN/dp_{T}", | |
86 | }; | |
87 | ||
88 | //__________________________________________________________ | |
89 | AliHFEdisplacedElectrons::AliHFEdisplacedElectrons(): | |
faee3b18 | 90 | fDeDebugLevel(0) |
91 | , fNclustersITS(0) | |
92 | , fMinNprimVtxContributor(1) | |
70da6c5a | 93 | , fTHnSparseDcaMcEleInfo(NULL) |
faee3b18 | 94 | , fTHnSparseDcaEsdEleInfo(NULL) |
70da6c5a | 95 | , fTHnSparseDcaDataEleInfo(NULL) |
faee3b18 | 96 | , fDeOutputList(0x0) |
70da6c5a | 97 | { |
98 | // | |
99 | // default constructor | |
100 | // | |
101 | ||
102 | } | |
103 | ||
104 | //__________________________________________________________ | |
91c7e1ec | 105 | AliHFEdisplacedElectrons::AliHFEdisplacedElectrons(const AliHFEdisplacedElectrons &ref): |
106 | TObject(ref) | |
faee3b18 | 107 | , fDeDebugLevel(ref.fDeDebugLevel) |
108 | , fNclustersITS(ref.fNclustersITS) | |
109 | , fMinNprimVtxContributor(ref.fMinNprimVtxContributor) | |
91c7e1ec | 110 | , fTHnSparseDcaMcEleInfo(ref.fTHnSparseDcaMcEleInfo) |
faee3b18 | 111 | , fTHnSparseDcaEsdEleInfo(ref.fTHnSparseDcaEsdEleInfo) |
91c7e1ec | 112 | , fTHnSparseDcaDataEleInfo(ref.fTHnSparseDcaDataEleInfo) |
faee3b18 | 113 | , fDeOutputList(ref.fDeOutputList) |
70da6c5a | 114 | |
115 | { | |
116 | // | |
117 | // copy constructor | |
118 | // | |
119 | } | |
120 | ||
121 | ||
122 | //__________________________________________________________ | |
91c7e1ec | 123 | AliHFEdisplacedElectrons&AliHFEdisplacedElectrons::operator=(const AliHFEdisplacedElectrons &ref) |
70da6c5a | 124 | { |
125 | // | |
126 | // Assignment operator | |
127 | // | |
91c7e1ec | 128 | |
129 | if(this == &ref) return *this; | |
130 | AliHFEdisplacedElectrons::operator=(ref); | |
131 | ||
faee3b18 | 132 | fDeDebugLevel = ref.fDeDebugLevel; |
133 | fNclustersITS = ref.fNclustersITS; | |
134 | fMinNprimVtxContributor = ref.fMinNprimVtxContributor; | |
91c7e1ec | 135 | fTHnSparseDcaMcEleInfo = ref.fTHnSparseDcaMcEleInfo; |
faee3b18 | 136 | fTHnSparseDcaEsdEleInfo = ref.fTHnSparseDcaEsdEleInfo; |
91c7e1ec | 137 | fTHnSparseDcaDataEleInfo = ref.fTHnSparseDcaDataEleInfo; |
faee3b18 | 138 | fDeOutputList = ref.fDeOutputList; |
91c7e1ec | 139 | |
70da6c5a | 140 | return *this; |
141 | } | |
142 | ||
143 | //__________________________________________________________ | |
144 | AliHFEdisplacedElectrons::~AliHFEdisplacedElectrons() | |
145 | { | |
146 | // | |
147 | // default constructor | |
148 | // | |
149 | ||
70da6c5a | 150 | if(fTHnSparseDcaMcEleInfo) |
151 | delete fTHnSparseDcaMcEleInfo; | |
faee3b18 | 152 | if(fTHnSparseDcaEsdEleInfo) |
153 | delete fTHnSparseDcaEsdEleInfo; | |
70da6c5a | 154 | if(fTHnSparseDcaDataEleInfo) |
155 | delete fTHnSparseDcaDataEleInfo; | |
156 | ||
faee3b18 | 157 | if(fDeOutputList){ |
158 | fDeOutputList->Clear(); | |
159 | delete fDeOutputList; | |
70da6c5a | 160 | |
161 | } | |
162 | ||
163 | Printf("analysis done\n"); | |
164 | } | |
165 | ||
166 | //__________________________________________________________ | |
167 | void AliHFEdisplacedElectrons::InitAnalysis(){ | |
168 | // | |
169 | // init analysis (no intialization yet) | |
170 | // | |
171 | ||
172 | ||
173 | Printf("initialize analysis\n"); | |
174 | ||
175 | } | |
176 | ||
177 | ||
178 | //__________________________________________________________ | |
179 | void AliHFEdisplacedElectrons::CreateOutputs(TList* const displacedList){ | |
180 | ||
181 | // | |
faee3b18 | 182 | // create output fDeOutputList |
70da6c5a | 183 | // |
184 | ||
185 | // THnSparseF | |
faee3b18 | 186 | // 8+? interested electron sources: others, photon conv, direct photon, pi0, eta, b, b->c, c + pion or missid, and missid pion |
187 | // 41 possible DCA cuts XY: | |
70da6c5a | 188 | // 13 pT bins |
189 | // 10 rapidity bins | |
190 | // 10 azimuthal angle phi bins | |
191 | ||
192 | ||
193 | if(!displacedList) return; | |
194 | ||
faee3b18 | 195 | fDeOutputList = displacedList; |
196 | fDeOutputList -> SetName("displacedElectrons"); | |
70da6c5a | 197 | |
70da6c5a | 198 | // electron source |
faee3b18 | 199 | Int_t nBinsEleSource = 10; |
70da6c5a | 200 | Double_t minEleSource = -1.5; |
69ac0e6f | 201 | Double_t maxEleSource = 8.5; |
70da6c5a | 202 | Double_t *binLimEleSource = new Double_t[nBinsEleSource+1]; |
203 | for(Int_t i=0; i<=nBinsEleSource; i++) | |
204 | binLimEleSource[i] = minEleSource + i*(maxEleSource-minEleSource)/nBinsEleSource; | |
205 | ||
faee3b18 | 206 | // dca bins: XY |
207 | Int_t nBinsDcaXY = kNDcaMin-1; // 41 bins | |
208 | Double_t minDcaXY = -20.5; | |
209 | Double_t maxDcaXY = 20.5; | |
210 | Double_t dcaXYBinWidth = (maxDcaXY-minDcaXY)/nBinsDcaXY; | |
211 | Double_t *binLimDcaXY = new Double_t[nBinsDcaXY+1]; | |
212 | for(Int_t i=0; i<=nBinsDcaXY; i++) | |
213 | binLimDcaXY[i] = minDcaXY + i*dcaXYBinWidth; | |
214 | ||
215 | ||
216 | // dca bins: Z | |
217 | Int_t nBinsDcaZ = kNDcaMin/2; // 21 bins | |
218 | Double_t minDcaZ = -10.5; | |
219 | Double_t maxDcaZ = 10.5; | |
220 | Double_t dcaZBinWidth = (maxDcaZ-minDcaZ)/nBinsDcaZ; | |
221 | Double_t *binLimDcaZ = new Double_t[nBinsDcaZ+1]; | |
222 | for(Int_t i=0; i<=nBinsDcaZ; i++) | |
223 | binLimDcaZ[i] = minDcaZ + i*dcaZBinWidth; | |
70da6c5a | 224 | |
70da6c5a | 225 | // pt bins |
226 | Int_t nBinsPt = kNPtIntv-1; | |
227 | Double_t *binLimPt = new Double_t[nBinsPt+1]; | |
228 | for(Int_t i=0; i<=nBinsPt; i++) | |
229 | binLimPt[i] = fgkPtIntv[i]; // variable bins | |
230 | ||
231 | // rapidity bins | |
232 | Int_t nBinsRap = 10; | |
233 | Double_t minRap = -1.0; | |
234 | Double_t maxRap = 1.0; | |
235 | Double_t *binLimRap = new Double_t[nBinsRap+1]; | |
236 | for(Int_t i=0; i<=nBinsRap; i++) | |
237 | binLimRap[i] = minRap + i*(maxRap-minRap)/nBinsRap; | |
238 | ||
239 | // azumuthal phi angle | |
240 | Int_t nBinsPhi = 10; | |
241 | Double_t minPhi = 0; | |
242 | Double_t maxPhi = 2*TMath::Pi(); | |
243 | Double_t *binLimPhi = new Double_t[nBinsPhi+1]; | |
244 | for(Int_t i=0; i<=nBinsPhi; i++) | |
245 | binLimPhi[i] = minPhi + i*(maxPhi-minPhi)/nBinsPhi; | |
246 | ||
69ac0e6f | 247 | // production radii |
248 | Int_t nBinsR = 30; | |
249 | Double_t minR = 0; | |
250 | Double_t maxR = 15; | |
251 | Double_t *binLimR = new Double_t[nBinsR+1]; | |
252 | for(Int_t i=0; i<=nBinsR; i++) | |
253 | binLimR[i] = minR + i*(maxR-minR)/nBinsR; | |
254 | ||
255 | // status | |
256 | Int_t nBinsStat = 11; | |
257 | Double_t minStat = 0; | |
258 | Double_t maxStat = 11; | |
259 | Double_t *binLimStat = new Double_t[nBinsStat+1]; | |
260 | for(Int_t i=0; i<=nBinsStat; i++) | |
261 | binLimStat[i] = minStat + i*(maxStat-minStat)/nBinsStat; | |
70da6c5a | 262 | |
faee3b18 | 263 | fTHnSparseDcaMcEleInfo = 0x0; |
264 | fTHnSparseDcaEsdEleInfo = 0x0; | |
265 | fTHnSparseDcaDataEleInfo = 0x0; | |
70da6c5a | 266 | |
faee3b18 | 267 | // for MC only: MC electron ID |
69ac0e6f | 268 | const Int_t nVarMc = 8; |
269 | Int_t iBinMc[nVarMc] = {nBinsEleSource, nBinsDcaXY, nBinsDcaZ, nBinsPt, nBinsRap, nBinsPhi, nBinsR, nBinsStat}; | |
270 | ||
faee3b18 | 271 | fTHnSparseDcaMcEleInfo = new THnSparseF("dcaMcElectronInfo", |
69ac0e6f | 272 | "MC electrons;electron source ID;mc dcaXY [50 #mum];mc dcaZ [100 #mum];mc pT [GeV/c];mc y [rapidity];mc #phi [rad];mc R [cm];Status Code", |
faee3b18 | 273 | nVarMc, iBinMc); |
69ac0e6f | 274 | |
faee3b18 | 275 | fTHnSparseDcaMcEleInfo->SetBinEdges(0, binLimEleSource); // electron source |
276 | fTHnSparseDcaMcEleInfo->SetBinEdges(1, binLimDcaXY); // dca xy cut | |
277 | fTHnSparseDcaMcEleInfo->SetBinEdges(2, binLimDcaZ); // dca z cut | |
278 | fTHnSparseDcaMcEleInfo->SetBinEdges(3, binLimPt); // pt | |
279 | fTHnSparseDcaMcEleInfo->SetBinEdges(4, binLimRap); // rapidity | |
280 | fTHnSparseDcaMcEleInfo->SetBinEdges(5, binLimPhi); // phi | |
69ac0e6f | 281 | fTHnSparseDcaMcEleInfo->SetBinEdges(6, binLimR); |
282 | fTHnSparseDcaMcEleInfo->SetBinEdges(7, binLimStat); | |
283 | fTHnSparseDcaMcEleInfo->Sumw2(); | |
284 | ||
faee3b18 | 285 | // for ESD with MC: HFE pid and MC pid |
69ac0e6f | 286 | const Int_t nVarEsd = 9; |
287 | Int_t iBin[nVarEsd] = {nBinsEleSource,nBinsDcaXY, nBinsDcaZ, nBinsDcaXY, nBinsPt, nBinsRap, nBinsPhi, nBinsR, nBinsStat}; | |
faee3b18 | 288 | |
289 | fTHnSparseDcaEsdEleInfo | |
290 | = new THnSparseF("dcaEsdElectronInfo", | |
69ac0e6f | 291 | "ESD electrons;electron source ID;esd dcaXY [50 #mum];esd dcaZ [100 #mum];esd dcaXY KF [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];R [cm];Status Code", |
faee3b18 | 292 | nVarEsd, iBin); |
293 | ||
294 | fTHnSparseDcaEsdEleInfo->SetBinEdges(0, binLimEleSource); // electron source | |
295 | fTHnSparseDcaEsdEleInfo->SetBinEdges(1, binLimDcaXY); // dca xy without current track | |
296 | fTHnSparseDcaEsdEleInfo->SetBinEdges(2, binLimDcaZ); // dca z without current track | |
297 | fTHnSparseDcaEsdEleInfo->SetBinEdges(3, binLimDcaXY); // dca xy kf without current track | |
298 | fTHnSparseDcaEsdEleInfo->SetBinEdges(4, binLimPt); // pt esd | |
299 | fTHnSparseDcaEsdEleInfo->SetBinEdges(5, binLimRap); // rapidity esd | |
300 | fTHnSparseDcaEsdEleInfo->SetBinEdges(6, binLimPhi); // phi esd | |
69ac0e6f | 301 | fTHnSparseDcaEsdEleInfo->SetBinEdges(7, binLimR); |
302 | fTHnSparseDcaEsdEleInfo->SetBinEdges(8, binLimStat); | |
faee3b18 | 303 | fTHnSparseDcaEsdEleInfo->Sumw2(); |
304 | ||
69ac0e6f | 305 | |
faee3b18 | 306 | // for ESD data: HFE pid |
307 | const Int_t nVarData = 6; | |
308 | Int_t iBinData[nVarData] = {nBinsDcaXY, nBinsDcaZ, nBinsDcaXY, nBinsPt, nBinsRap, nBinsPhi}; | |
309 | ||
310 | fTHnSparseDcaDataEleInfo | |
311 | = new THnSparseF("dcaDataElectronInfo", | |
312 | "Data electrons;dcaXY [50 #mum];dcaZ [100 #mum];dcaXYkf [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];", | |
313 | nVarData, iBinData); | |
314 | fTHnSparseDcaDataEleInfo->SetBinEdges(0, binLimDcaXY); // dca xy cut w/o | |
315 | fTHnSparseDcaDataEleInfo->SetBinEdges(1, binLimDcaZ); // dca z cut w/o | |
316 | fTHnSparseDcaDataEleInfo->SetBinEdges(2, binLimDcaXY); // dca xy kf cut | |
317 | fTHnSparseDcaDataEleInfo->SetBinEdges(3, binLimPt); // pt | |
318 | fTHnSparseDcaDataEleInfo->SetBinEdges(4, binLimRap); // rapidity | |
319 | fTHnSparseDcaDataEleInfo->SetBinEdges(5, binLimPhi); // phi | |
320 | fTHnSparseDcaDataEleInfo->Sumw2(); | |
70da6c5a | 321 | |
faee3b18 | 322 | fDeOutputList -> AddAt(fTHnSparseDcaMcEleInfo, kMcElectron); |
323 | fDeOutputList -> AddAt(fTHnSparseDcaEsdEleInfo, kEsdElectron); | |
324 | fDeOutputList -> AddAt(fTHnSparseDcaDataEleInfo, kDataElectron); | |
325 | ||
326 | AliInfo("THnSparse histograms are created\n"); | |
327 | fDeOutputList->Print(); | |
328 | ||
329 | } | |
330 | ||
331 | ||
332 | //__________________________________________________________ | |
c2690925 | 333 | void AliHFEdisplacedElectrons::FillMcOutput(const AliESDEvent *const fESD, AliMCEvent* const fMC, const AliMCParticle* const mctrack) |
faee3b18 | 334 | { |
70da6c5a | 335 | |
faee3b18 | 336 | // fill output |
337 | //0. after mc event cut | |
338 | //1. after checking stack, mcpart etc are valid | |
339 | //2. after PID | |
340 | //3. after event and track selection | |
341 | ||
342 | AliStack *stack = fMC->Stack(); | |
343 | TParticle *part = mctrack->Particle(); | |
344 | ||
345 | // obtain impact parameters in xy and z | |
346 | AliMCVertex *mcPrimVtx = (AliMCVertex *)fMC->GetPrimaryVertex(); | |
347 | Double_t mcPrimV[3]; | |
348 | mcPrimV[0] = mcPrimVtx->GetX(); | |
349 | mcPrimV[1] = mcPrimVtx->GetY(); | |
350 | mcPrimV[2] = mcPrimVtx->GetZ(); | |
351 | ||
69ac0e6f | 352 | Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]); |
faee3b18 | 353 | |
354 | Float_t vx = part->Vx(); // in cm | |
355 | Float_t vy = part->Vy(); // in cm | |
356 | Float_t vz = part->Vz(); // in cm | |
357 | ||
358 | Float_t vxy = TMath::Sqrt(vx*vx+vy*vy); | |
70da6c5a | 359 | |
faee3b18 | 360 | Float_t mcpx = part->Px(); |
361 | Float_t mcpy = part->Py(); | |
362 | Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy); | |
363 | Float_t mceta = part->Eta(); | |
364 | Float_t mcphi = part->Phi(); | |
70da6c5a | 365 | |
69ac0e6f | 366 | // begin new stuff |
367 | Double_t mcR = part->R(); | |
368 | Int_t mcStatus = part->GetStatusCode(); | |
369 | // end new staff | |
370 | ||
faee3b18 | 371 | Int_t pdg = part->GetPdgCode(); |
70da6c5a | 372 | |
faee3b18 | 373 | Int_t charge = 1; |
374 | if(pdg==kPDGelectron || pdg==-kPDGpion) charge = -1; | |
70da6c5a | 375 | |
faee3b18 | 376 | // calculate mcDca ------------------------------------------------------------------ |
377 | const Float_t conv[2] = {1.783/1.6, 2.99792458}; | |
378 | Float_t magneticField = 0; // initialized as 5kG | |
379 | magneticField = fESD->GetMagneticField(); // in kG | |
380 | Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter | |
70da6c5a | 381 | |
faee3b18 | 382 | Float_t radius; |
383 | radius = TMath::Abs(radiusMc); | |
384 | Float_t nx = mcpx/mcpt; | |
385 | Float_t ny = mcpy/mcpt; | |
386 | Double_t dxy = vxy - mcVtxXY; // in cm | |
387 | Double_t dvx = vx - mcPrimV[0]; // in cm | |
388 | Double_t dvy = vy - mcPrimV[1]; // in cm | |
389 | Float_t mcDcaXYm = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters | |
390 | Double_t mcDca[2] = {mcDcaXYm*100, vz}; // in cm | |
391 | Double_t mcDcaXY = mcDca[0]*1.0e4; // conv dca in cm to dca in micron | |
392 | Double_t mcDcaZ = mcDca[1]*1.0e4; | |
393 | ||
69ac0e6f | 394 | const Int_t nvarMC=8; |
faee3b18 | 395 | Double_t var[nvarMC]; |
396 | var[0] = -1; | |
397 | ||
398 | if(TMath::Abs(pdg)==kPDGelectron) { | |
70da6c5a | 399 | |
faee3b18 | 400 | Int_t eleLabel = mctrack->GetLabel(); |
401 | Int_t sourcePdg = ElectronFromSource(stack, eleLabel); | |
402 | ||
403 | if(sourcePdg==kPDGgamma){ // check direct photon or not | |
404 | if(ElePhotonDirect(stack, eleLabel)!=-1) | |
405 | var[0] = kEleDirectPhotonConv; | |
406 | else | |
407 | var[0] = kElePhotonConv; | |
408 | if(fDeDebugLevel>=10) | |
409 | printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]); | |
410 | } // photon or direct photon -> e | |
411 | ||
412 | if(sourcePdg==kPDGpi0){ | |
413 | var[0] = kElePi0; | |
414 | if(fDeDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]); | |
415 | } | |
416 | ||
417 | if(sourcePdg==kPDGeta){ | |
418 | var[0] = kEleEta; | |
419 | if(fDeDebugLevel>=10) | |
420 | printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]); | |
421 | } // from eta -> e | |
422 | ||
423 | if(TMath::Abs(sourcePdg%10000)/100==kPDGbeauty || // for special intermediate meson states: like 100553 | |
424 | TMath::Abs(sourcePdg)/1000==kPDGbeauty || | |
425 | TMath::Abs(sourcePdg)/100==kPDGbeauty || | |
426 | TMath::Abs(sourcePdg)==kPDGbeauty){ | |
427 | var[0]=kEleB; | |
428 | if(fDeDebugLevel>=10) | |
429 | printf("beauty======electron %d is from %d with id %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]); | |
430 | } // direct beauty -> e | |
431 | if(TMath::Abs(sourcePdg%10000)/100==kPDGcharm || // for special intermediate meson states: like 1004** | |
432 | TMath::Abs(sourcePdg)/1000==kPDGcharm || | |
433 | TMath::Abs(sourcePdg)/100==kPDGcharm || | |
434 | TMath::Abs(sourcePdg)==kPDGcharm){ | |
435 | // two cases: | |
436 | // electron from b->c->e | |
437 | // electron from c->e | |
438 | if(ElectronFromCharm(stack, eleLabel)!=-1){ | |
439 | var[0] = ElectronFromCharm(stack, eleLabel); | |
440 | if(fDeDebugLevel>=10) | |
441 | printf("charm----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]); | |
442 | } | |
443 | } // charm electrons: b->c->e or c->e | |
444 | ||
445 | if(fDeDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]); | |
446 | } // electron source endif | |
447 | ||
448 | else | |
449 | if(TMath::Abs(pdg)==kPDGpion) | |
450 | var[0] = kPion; | |
451 | ||
452 | if(TMath::Abs(mcDcaXY)<1000) var[1] = Int_t(mcDcaXY)/50; // larger than 1mm should go to the last bin | |
453 | else | |
454 | var[1] = ((mcDcaXY>0)?1:-1)*20; | |
455 | ||
456 | if(TMath::Abs(mcDcaZ)<1000) var[2] = Int_t(mcDcaZ)/100; | |
457 | else | |
458 | var[2] = ((mcDcaZ>0)?1:-1)*10; | |
459 | ||
69ac0e6f | 460 | var[3] = mcpt; // pt |
461 | var[4] = mceta; // eta | |
462 | var[5] = mcphi; // phi | |
463 | var[6] = mcR; // production radius | |
464 | var[7] = mcStatus; // internal status | |
465 | ||
faee3b18 | 466 | |
e97c2edf | 467 | (static_cast<THnSparseF *>(fDeOutputList->At(kMcElectron)))->Fill(var); |
70da6c5a | 468 | } |
469 | ||
470 | ||
faee3b18 | 471 | |
70da6c5a | 472 | //__________________________________________________________ |
c2690925 | 473 | void AliHFEdisplacedElectrons::FillEsdOutput(const AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack, AliStack * const stack) |
70da6c5a | 474 | { |
faee3b18 | 475 | // after esd event selection, esd track cuts, and hfe pid |
476 | // this is the case for ESD tracks, with MC information | |
70da6c5a | 477 | |
70da6c5a | 478 | AliESDtrack *track = esdTrack; |
69ac0e6f | 479 | Double_t pt = track->Pt(); |
70da6c5a | 480 | Double_t eta = track->Eta(); |
481 | Double_t phi = track->Phi(); | |
482 | ||
faee3b18 | 483 | Int_t eleLabel = track->GetLabel(); |
484 | if(eleLabel<0 || eleLabel>stack->GetNtrack()) return; | |
70da6c5a | 485 | |
faee3b18 | 486 | TParticle *particle = stack->Particle(eleLabel); |
70da6c5a | 487 | if(!particle) return; |
70da6c5a | 488 | |
69ac0e6f | 489 | // begin new stuff |
490 | Double_t mcR = particle->R(); | |
491 | Int_t mcStatus = particle->GetStatusCode(); | |
492 | // end new staff | |
493 | ||
faee3b18 | 494 | // obtain impact parameters in xy and z |
70da6c5a | 495 | Float_t magneticField = 5; // initialized as 5kG |
496 | magneticField = fESDEvent->GetMagneticField(); // in kG | |
faee3b18 | 497 | Double_t beampiperadius=3.; |
498 | const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex(); | |
70da6c5a | 499 | |
faee3b18 | 500 | |
69ac0e6f | 501 | const Int_t nvarESD = 9; |
faee3b18 | 502 | Double_t var[nvarESD]; |
503 | var[0] = -1; | |
70da6c5a | 504 | |
faee3b18 | 505 | Int_t sourcePdg = -1; |
70da6c5a | 506 | |
69ac0e6f | 507 | if(TMath::Abs(particle->GetPdgCode())==kPDGelectron){ |
508 | ||
faee3b18 | 509 | sourcePdg = ElectronFromSource(stack, eleLabel); |
70da6c5a | 510 | |
faee3b18 | 511 | if(sourcePdg==kPDGgamma){ // check direct photon or not |
70da6c5a | 512 | if(ElePhotonDirect(stack, eleLabel)!=-1) |
513 | var[0] = kEleDirectPhotonConv; | |
514 | else | |
515 | var[0] = kElePhotonConv; | |
faee3b18 | 516 | if(fDeDebugLevel>=10) |
70da6c5a | 517 | printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]); |
69ac0e6f | 518 | } else if(sourcePdg==kPDGpi0){ // check dalitz |
70da6c5a | 519 | var[0] = kElePi0; |
faee3b18 | 520 | if(fDeDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]); |
69ac0e6f | 521 | } else if(sourcePdg==kPDGeta){ // check eta |
70da6c5a | 522 | var[0] = kEleEta; |
faee3b18 | 523 | if(fDeDebugLevel>=10) |
70da6c5a | 524 | printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]); |
69ac0e6f | 525 | } else if(IsB(sourcePdg)) var[0] = kEleB; // check beauty |
526 | else if(IsC(sourcePdg)) var[0] = CheckCharm(stack,eleLabel); // check charm | |
70da6c5a | 527 | |
69ac0e6f | 528 | } else if(TMath::Abs(particle->GetPdgCode())==kPDGpion) var[0] = kEleMissIDpion; |
529 | else var[0] = kEleMissID; | |
530 | // ---- PID END ---- | |
531 | ||
532 | ||
faee3b18 | 533 | // excluding current track |
534 | // ---- beginning --- method from Andrea D 28.05.2010 | |
535 | AliVertexerTracks *vertexer = new AliVertexerTracks(magneticField); | |
536 | vertexer->SetITSMode(); | |
537 | vertexer->SetMinClusters(fNclustersITS); | |
538 | Int_t skipped[2]; | |
539 | skipped[0] = (Int_t)track->GetID(); | |
540 | vertexer->SetSkipTracks(1,skipped); | |
541 | AliESDVertex *vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESDEvent); | |
542 | delete vertexer; vertexer = NULL; | |
543 | if(vtxESDSkip->GetNContributors()<fMinNprimVtxContributor) return; | |
544 | // -- ending --- method from Andrea D 28.05.2010 | |
545 | ||
546 | Double_t dz[2]; // error of dca in cm | |
547 | Double_t covardz[3]; | |
548 | if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection | |
549 | ||
550 | Double_t dcaXY = dz[0]*1.0e4; // conv dca in cm to dca in micron | |
551 | Double_t dcaZ = dz[1]*1.0e4; | |
552 | ||
553 | if(fDeDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]); | |
554 | ||
555 | if(TMath::Abs(dcaXY)<1000) var[1] = Int_t(dcaXY)/50; // larger than 1mm should go to the last bin | |
556 | else | |
557 | var[1] = ((dcaXY>0)?1:-1)*20; | |
558 | ||
559 | if(TMath::Abs(dcaZ)<1000) var[2] = Int_t(dcaZ)/100; | |
560 | else | |
561 | var[2] = ((dcaZ>0)?1:-1)*10; | |
562 | ||
563 | // calculate dca using AliKFParticle class------------------------------------------------------------------ | |
564 | Float_t kfDcaXY = 0; | |
565 | Int_t trkID = track->GetID(); | |
566 | AliKFParticle::SetField(magneticField); | |
567 | AliKFParticle kfParticle(*track, particle->GetPdgCode()); | |
568 | // prepare kfprimary vertex | |
569 | AliKFVertex kfESDprimary; | |
570 | // Reconstruct Primary Vertex (with ESD tracks) | |
571 | Int_t n=primVtx->GetNIndices(); | |
572 | if (n>0 && primVtx->GetStatus()){ | |
573 | kfESDprimary = AliKFVertex(*primVtx); | |
574 | UShort_t *priIndex = primVtx->GetIndices(); | |
575 | for (Int_t i=0;i<n;i++){ | |
576 | Int_t idx = Int_t(priIndex[i]); | |
577 | if (idx == trkID){ | |
578 | kfESDprimary -= kfParticle; | |
579 | kfDcaXY = kfParticle.GetDistanceFromVertexXY(kfESDprimary)*1e4; | |
580 | } // remove current track from this calculation | |
581 | } // loop over all primary vertex contributors | |
70da6c5a | 582 | } |
faee3b18 | 583 | // end of KF dca |
584 | ||
585 | if(TMath::Abs(kfDcaXY)<1000) | |
586 | var[3] = Int_t(kfDcaXY)/50; | |
587 | else | |
588 | var[3] = (kfDcaXY>0?1:-1)*20; | |
589 | ||
590 | var[4] = pt; // pt | |
591 | var[5] = eta; // eta | |
592 | var[6] = phi; // phi | |
69ac0e6f | 593 | var[7] = mcR; |
594 | var[8] = mcStatus; | |
faee3b18 | 595 | |
e97c2edf | 596 | (static_cast<THnSparseF *>(fDeOutputList->At(kEsdElectron)))->Fill(var); |
70da6c5a | 597 | |
598 | } | |
599 | ||
70da6c5a | 600 | //__________________________________________________________ |
c2690925 | 601 | void AliHFEdisplacedElectrons::FillDataOutput(const AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack) |
70da6c5a | 602 | { |
faee3b18 | 603 | |
604 | // this is pure data, without MC information at all | |
605 | // fill output: with HFE pid selection of electrons after all track quality cuts | |
70da6c5a | 606 | |
70da6c5a | 607 | if(!esdTrack) return; |
608 | AliESDtrack *track = esdTrack; | |
faee3b18 | 609 | |
70da6c5a | 610 | Double_t pt = track->Pt(); |
611 | Double_t eta = track->Eta(); | |
612 | Double_t phi = track->Phi(); | |
613 | ||
614 | // obtain impact parameters in xy and y | |
91c7e1ec | 615 | const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex(); |
70da6c5a | 616 | |
617 | Float_t magneticField = 5; // initialized as 5kG | |
91c7e1ec | 618 | magneticField = fESDEvent->GetMagneticField(); // in kG |
faee3b18 | 619 | Double_t beampiperadius=3.; |
620 | ||
621 | ||
622 | const Int_t nvarData = 6; | |
623 | Double_t varData[nvarData]; | |
624 | ||
625 | // | |
626 | // excluding current track | |
627 | // | |
628 | ||
629 | //------ beginning --- method from Andrea D 28.05.2010 | |
630 | AliVertexerTracks *vertexer = new AliVertexerTracks(fESDEvent->GetMagneticField()); | |
631 | vertexer->SetITSMode(); | |
632 | vertexer->SetMinClusters(fNclustersITS); | |
633 | Int_t skipped[2]; | |
634 | skipped[0] = (Int_t)track->GetID(); | |
635 | vertexer->SetSkipTracks(1,skipped); | |
636 | AliESDVertex *vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESDEvent); | |
637 | delete vertexer; vertexer = NULL; | |
638 | if(vtxESDSkip->GetNContributors()<fMinNprimVtxContributor) return; | |
639 | //------ ending --- method from Andrea D 28.05.2010 | |
640 | ||
70da6c5a | 641 | Double_t dz[2]; // error of dca in cm |
642 | Double_t covardz[3]; | |
faee3b18 | 643 | |
644 | if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection | |
70da6c5a | 645 | |
faee3b18 | 646 | Double_t dcaXY = dz[0]*1.0e4; // conv dca in cm to dca in micron |
647 | Double_t dcaZ = dz[1]*1.0e4; | |
648 | ||
649 | if(TMath::Abs(dcaXY)<1000) varData[0] = Int_t(dcaXY)/50; // larger than 1mm should go to the last bin | |
650 | else | |
651 | varData[0] = (dcaXY>0?1:-1)*20; | |
652 | if(TMath::Abs(dcaZ)<1000) varData[1] = Int_t(dcaZ)/100; | |
653 | else | |
654 | varData[1] = (dcaZ>0?1:-1)*10; | |
655 | ||
656 | ||
657 | // calculate dca using AliKFParticle class------------------------------------------------------------------ | |
658 | Float_t kfDcaXY = 0; | |
659 | Int_t trkID = track->GetID(); | |
660 | AliKFParticle::SetField(magneticField); | |
661 | AliKFParticle kfParticle(*track, -11*track->Charge()); | |
662 | // prepare kfprimary vertex | |
663 | AliKFVertex kfESDprimary; | |
664 | // Reconstruct Primary Vertex (with ESD tracks) | |
665 | Int_t n=primVtx->GetNIndices(); | |
666 | if (n>0 && primVtx->GetStatus()){ | |
667 | kfESDprimary = AliKFVertex(*primVtx); | |
668 | UShort_t *priIndex = primVtx->GetIndices(); | |
669 | for (Int_t i=0;i<n;i++){ | |
670 | Int_t idx = Int_t(priIndex[i]); | |
671 | if (idx == trkID){ | |
672 | kfESDprimary -= kfParticle; | |
673 | kfDcaXY = kfParticle.GetDistanceFromVertexXY(kfESDprimary)*1e4; | |
674 | } // remove current track from this calculation | |
675 | } // loop over all primary vertex contributors | |
676 | } | |
677 | if(TMath::Abs(kfDcaXY)<1000) | |
678 | varData[2] = Int_t(kfDcaXY)/50; | |
70da6c5a | 679 | else |
faee3b18 | 680 | varData[2] = ((kfDcaXY)>0?1:-1)*20; |
681 | // end of KF dca | |
70da6c5a | 682 | |
faee3b18 | 683 | varData[3] = pt; // pt |
684 | varData[4] = eta; //eta | |
685 | varData[5] = phi; // phi | |
70da6c5a | 686 | |
e97c2edf | 687 | (static_cast<THnSparseF *>(fDeOutputList->At(kDataElectron)))->Fill(varData); |
70da6c5a | 688 | |
689 | } | |
690 | ||
691 | //__________________________________________________________ | |
692 | Int_t AliHFEdisplacedElectrons::ElectronFromSource(AliStack * const stack, Int_t label) const | |
693 | { | |
694 | ||
695 | // | |
696 | // return electron source label via electron label | |
697 | // | |
698 | ||
699 | // kEleSource is supposed to be either photon conv, direct photon conv, pi0, eta, beauty --> 0, 1, 2, 3, 4 | |
700 | ||
701 | Int_t eleLabel = label; | |
702 | ||
703 | if(eleLabel<0 || eleLabel>stack->GetNtrack()) return -1; | |
704 | ||
705 | TParticle *particle = stack->Particle(eleLabel); | |
706 | if(!particle) return -1; | |
707 | ||
708 | Int_t motherLabel = particle->GetFirstMother(); | |
709 | if(motherLabel<0 || motherLabel>stack->GetNtrack()) return -1; | |
710 | TParticle *motherPart = stack->Particle(motherLabel); | |
711 | if(!motherPart) return -1; | |
712 | ||
713 | Int_t pdgCode = TMath::Abs(motherPart->GetPdgCode()); | |
714 | ||
715 | if(pdgCode==kPDGelectron) { | |
faee3b18 | 716 | if(fDeDebugLevel>=10) printf("particle label: %d...(motherLabel=%d : motherPdg=%d) grandmother's pdg code was returned...%d \n", |
70da6c5a | 717 | label, motherLabel, pdgCode, ElectronFromSource(stack, motherLabel)); |
718 | return ElectronFromSource(stack, motherLabel); | |
719 | } | |
720 | ||
721 | return pdgCode; | |
722 | ||
723 | } | |
724 | ||
725 | //__________________________________________________________ | |
726 | Int_t AliHFEdisplacedElectrons::ElectronFromCharm(AliStack * const stack, Int_t label) const | |
727 | { | |
728 | // | |
729 | // separate electron: kEleC from c->eX, kEleBC from b->c->eX | |
730 | // | |
731 | ||
732 | Int_t motherLabel = label; | |
733 | TParticle *motherParticle = stack->Particle(motherLabel); // mother part | |
734 | if(!motherParticle) return -1; | |
735 | Int_t gMotherLabel = motherParticle->GetFirstMother(); // grand mother | |
736 | if(gMotherLabel<0 || gMotherLabel>stack->GetNtrack()) return -1; | |
737 | ||
738 | TParticle *gMotherPart = stack->Particle(gMotherLabel); | |
739 | if(!gMotherPart) return -1; | |
740 | ||
741 | Int_t pdgCode = gMotherPart->GetPdgCode(); | |
742 | if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty || // for special intermediate meson states: like 100553 | |
743 | TMath::Abs(pdgCode)/1000==kPDGbeauty || TMath::Abs(pdgCode)/100==kPDGbeauty || TMath::Abs(pdgCode)==kPDGbeauty) { | |
faee3b18 | 744 | if(fDeDebugLevel>=10) printf("this electron label %d is from mother %d, and finally from %d\n", label, ElectronFromSource(stack, label), pdgCode); |
70da6c5a | 745 | return kEleBC; |
746 | } // for sure it is from BC | |
747 | ||
748 | else | |
749 | if(TMath::Abs(pdgCode%10000)/100==kPDGcharm || // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4) | |
750 | TMath::Abs(pdgCode)/1000==kPDGcharm || | |
751 | TMath::Abs(pdgCode)/100==kPDGcharm || | |
752 | TMath::Abs(pdgCode)==kPDGcharm){ | |
753 | ||
69ac0e6f | 754 | if(CharmFromBeauty(stack, gMotherLabel)!=-1) return kEleBC; |
755 | else return kEleC; | |
70da6c5a | 756 | |
757 | } | |
758 | ||
759 | else | |
760 | return -1; | |
761 | ||
762 | } | |
763 | ||
764 | //__________________________________________________________ | |
765 | Int_t AliHFEdisplacedElectrons::CharmFromBeauty(AliStack * const stack, Int_t hfLabel) const | |
766 | { | |
767 | ||
768 | // | |
769 | // check if charm meson/hadron is from beauty decay | |
770 | // -1, not from beauty | |
771 | // return the label of the mother of this charm hadron | |
772 | // | |
773 | ||
774 | Int_t charmLabel = hfLabel; | |
775 | ||
776 | // AliStack *stack = fMC->Stack(); | |
777 | if(charmLabel<0 || charmLabel>stack->GetNtrack()) return -1; | |
778 | TParticle *particle = stack->Particle(charmLabel); | |
779 | if(!particle) return -1; | |
780 | ||
781 | Int_t motherCharmLabel = particle->GetFirstMother(); | |
782 | if(motherCharmLabel<0 || motherCharmLabel>stack->GetNtrack()) return -1; | |
783 | ||
784 | TParticle *motherCharmPart = stack->Particle(motherCharmLabel); | |
785 | if(!motherCharmPart) return -1; | |
786 | ||
787 | Int_t pdgCode = motherCharmPart->GetPdgCode(); | |
788 | ||
789 | if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty || // for special intermediate meson states: like 100553 | |
790 | TMath::Abs(pdgCode)/1000==kPDGbeauty || | |
791 | TMath::Abs(pdgCode)/100==kPDGbeauty || | |
792 | TMath::Abs(pdgCode)==kPDGbeauty) | |
793 | return motherCharmLabel; | |
794 | else | |
795 | if(TMath::Abs(pdgCode%10000)/100==kPDGcharm || // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4) | |
796 | TMath::Abs(pdgCode)/1000==kPDGcharm || | |
797 | TMath::Abs(pdgCode)/100==kPDGcharm) | |
798 | return CharmFromBeauty(stack, motherCharmLabel); // loop over to see if charm is from beauty -- if yes, return the label of mother of the charm | |
799 | ||
800 | else | |
801 | return -1; | |
802 | } | |
803 | ||
804 | //__________________________________________________________ | |
805 | Int_t AliHFEdisplacedElectrons::ElePhotonDirect(AliStack * const stack, Int_t label) const | |
806 | { | |
807 | // | |
808 | // electron is from photon, and check if this photon is direct | |
809 | // | |
810 | ||
811 | Int_t eleLabel = label; | |
812 | ||
813 | TParticle *particle = stack->Particle(eleLabel); | |
814 | if(particle->GetFirstMother()<0 || particle->GetFirstMother()>stack->GetNtrack()) | |
815 | return -1; | |
816 | ||
817 | TParticle *photonPart = stack->Particle(particle->GetFirstMother()); | |
818 | if(!photonPart) | |
819 | return -1; | |
820 | Int_t motherPhotonLabel = photonPart->GetFirstMother(); | |
821 | if(motherPhotonLabel<0 || motherPhotonLabel>stack->GetNtrack()) | |
822 | return -1; | |
823 | ||
824 | TParticle *motherPhotonPart = stack->Particle(motherPhotonLabel); | |
825 | if(!motherPhotonPart) | |
826 | return -1; | |
827 | ||
828 | Int_t pdgMotherPhoton = motherPhotonPart->GetPdgCode(); | |
829 | if(TMath::Abs(pdgMotherPhoton)<=10 || TMath::Abs(pdgMotherPhoton)==21) | |
830 | return 1; | |
831 | ||
832 | else | |
833 | return -1; | |
834 | ||
835 | ||
836 | } | |
837 | ||
838 | //__________________________________________________________ | |
839 | Int_t AliHFEdisplacedElectrons::GetMCpid(AliStack* const stack, Int_t partLabel) const | |
840 | { | |
841 | ||
842 | // | |
843 | // Simply pdg code | |
844 | // | |
845 | ||
846 | Int_t label = partLabel; | |
847 | // AliStack* stack = fMC->Stack(); | |
848 | if((label < 0) || (label >= stack->GetNtrack())) return -1; | |
849 | ||
850 | // MC Information | |
851 | TParticle * particle = stack->Particle(label); | |
852 | if(!particle) return -1; | |
853 | Int_t pdg = particle->GetPdgCode(); | |
854 | ||
855 | return pdg; | |
856 | ||
857 | } | |
858 | ||
859 | //__________________________________________________________ | |
860 | Int_t AliHFEdisplacedElectrons::GetMotherLabel(AliStack *const stack, Int_t eleLabel) const | |
861 | { | |
862 | // | |
863 | // Simply label of mother | |
864 | // | |
865 | ||
866 | ||
867 | Int_t label = eleLabel; | |
868 | // AliStack* stack = fMC->Stack(); | |
869 | if((label < 0) || (label >= stack->GetNtrack())) return -1; | |
870 | ||
871 | // MC Information | |
872 | TParticle * particle = stack->Particle(label); | |
873 | if(!particle) return -1; | |
874 | ||
875 | return particle->GetFirstMother(); | |
876 | ||
877 | } | |
878 | ||
879 | ||
880 | //__________________________________________________________ | |
881 | Float_t AliHFEdisplacedElectrons::GetRapidity(TParticle *part) const | |
882 | { | |
883 | // return rapidity | |
884 | ||
885 | Float_t rapidity = -999; | |
886 | if((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0) | |
887 | rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); | |
888 | ||
889 | return rapidity; | |
890 | } | |
891 | ||
892 | ||
893 | //__________________________________________________________ | |
894 | Float_t AliHFEdisplacedElectrons::GetTrackRapidity(AliESDtrack * const track) const | |
895 | { | |
896 | // return rapidity of electron | |
897 | ||
898 | Float_t px = track->Px(); | |
899 | Float_t py = track->Py(); | |
900 | Float_t pz = track->Pz(); | |
901 | ||
902 | // electron mass 0.00051099906 GeV/c2 | |
903 | TParticlePDG* electron = TDatabasePDG::Instance()->GetParticle(kPDGelectron); | |
904 | Double_t mass = electron->Mass(); | |
905 | ||
906 | Float_t en = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass); | |
907 | ||
908 | Float_t rapidity = -999; | |
909 | if((en - pz)*(en + pz)>0) | |
910 | rapidity = 0.5*(TMath::Log((en - pz)*(en + pz))); | |
911 | ||
912 | return rapidity; | |
913 | } | |
69ac0e6f | 914 | |
915 | //__________________________________________________________ | |
916 | Int_t AliHFEdisplacedElectrons::CheckCharm(AliStack * const stack, Int_t eleLabel) | |
917 | { | |
918 | // checks the genealogy of the ele from charm for a beauty | |
919 | // this method needs the stack and the label of the electron | |
920 | // warning: it assumes that the ele comes from a charm | |
921 | ||
922 | TParticle *particle = stack->Particle(eleLabel); | |
923 | Int_t label = particle->GetFirstMother(); | |
924 | ||
925 | while(label>0 && label<(stack->GetNtrack())){ | |
926 | particle = stack->Particle(label); | |
927 | if(IsB(TMath::Abs(particle->GetPdgCode()))) return kEleBC; | |
928 | label = particle->GetFirstMother(); | |
929 | } | |
930 | ||
931 | return kEleC; | |
932 | } | |
933 | ||
934 | //__________________________________________________________ | |
6555e2ad | 935 | Bool_t AliHFEdisplacedElectrons::IsB(Int_t pdg) const |
69ac0e6f | 936 | { |
937 | // check if the pdg is that of a beauty particle | |
938 | ||
939 | if((TMath::Abs(pdg)%10000)/100==kPDGbeauty || | |
940 | TMath::Abs(pdg)/1000==kPDGbeauty || | |
941 | TMath::Abs(pdg)/100==kPDGbeauty || | |
942 | TMath::Abs(pdg)==kPDGbeauty) return kTRUE; | |
943 | else return kFALSE; | |
944 | } | |
945 | ||
946 | //__________________________________________________________ | |
6555e2ad | 947 | Bool_t AliHFEdisplacedElectrons::IsC(Int_t pdg) const |
69ac0e6f | 948 | { |
949 | // check if the pdg is that of a charmed particle | |
950 | ||
951 | if((TMath::Abs(pdg)%10000)/100==kPDGcharm || | |
952 | TMath::Abs(pdg)/1000==kPDGcharm || | |
953 | TMath::Abs(pdg)/100==kPDGcharm || | |
954 | TMath::Abs(pdg)==kPDGcharm) return kTRUE; | |
955 | else return kFALSE; | |
956 | } |