]>
Commit | Line | Data |
---|---|---|
ffbede40 | 1 | /************************************************************************* |
2 | * Copyright(c) 1998-2009, 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 | // | |
17 | // Dielectron helper functions wrapped in a namespace | |
18 | // | |
19 | // | |
20 | // Authors: | |
21 | // Jens Wiechula <Jens.Wiechula@cern.ch> | |
ba15fdfb | 22 | // Frederick Kramer <Frederick.Kramer@cern.ch> |
41a64bfb | 23 | // Julian Book <Julian.Book@cern.ch> |
ffbede40 | 24 | |
25 | ||
26 | ||
27 | ||
28 | #include <TError.h> | |
29 | #include <TMath.h> | |
30 | #include <TObjString.h> | |
31 | #include <TObjArray.h> | |
32 | #include <TVectorD.h> | |
5720c765 | 33 | #include <TF1.h> |
34 | #include <TRandom.h> | |
35 | #include <TProfile.h> | |
ffbede40 | 36 | |
1201a1a9 | 37 | #include <AliVEvent.h> |
ba15fdfb | 38 | #include <AliVParticle.h> |
1201a1a9 | 39 | #include <AliKFParticle.h> |
ba15fdfb | 40 | #include <AliESDtrackCuts.h> |
41 | #include <AliESDEvent.h> | |
42 | #include <AliMCEvent.h> | |
43 | #include <AliAODEvent.h> | |
44 | #include <AliAODTracklets.h> | |
45 | #include <AliMultiplicity.h> | |
46 | #include <AliStack.h> | |
1201a1a9 | 47 | |
a823f01b | 48 | #include "AliDielectronVarCuts.h" |
49 | #include "AliDielectronTrackCuts.h" | |
5720c765 | 50 | #include "AliDielectronVarManager.h" |
ffbede40 | 51 | #include "AliDielectronHelper.h" |
52 | ||
53 | //_____________________________________________________________________________ | |
54 | TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) | |
55 | { | |
56 | // | |
57 | // Make logarithmic binning | |
58 | // the user has to delete the array afterwards!!! | |
59 | // | |
60 | ||
61 | //check limits | |
62 | if (xmin<1e-20 || xmax<1e-20){ | |
63 | Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!"); | |
64 | return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax); | |
65 | } | |
66 | if (xmax<xmin){ | |
67 | Double_t tmp=xmin; | |
68 | xmin=xmax; | |
69 | xmax=tmp; | |
70 | } | |
71 | TVectorD *binLim=new TVectorD(nbinsX+1); | |
72 | Double_t first=xmin; | |
73 | Double_t last=xmax; | |
74 | Double_t expMax=TMath::Log(last/first); | |
75 | for (Int_t i=0; i<nbinsX+1; ++i){ | |
76 | (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i); | |
77 | } | |
78 | return binLim; | |
79 | } | |
80 | ||
81 | //_____________________________________________________________________________ | |
82 | TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) | |
83 | { | |
84 | // | |
85 | // Make linear binning | |
86 | // the user has to delete the array afterwards!!! | |
87 | // | |
88 | if (xmax<xmin){ | |
89 | Double_t tmp=xmin; | |
90 | xmin=xmax; | |
91 | xmax=tmp; | |
92 | } | |
93 | TVectorD *binLim=new TVectorD(nbinsX+1); | |
94 | Double_t first=xmin; | |
95 | Double_t last=xmax; | |
96 | Double_t binWidth=(last-first)/nbinsX; | |
97 | for (Int_t i=0; i<nbinsX+1; ++i){ | |
98 | (*binLim)[i]=first+binWidth*(Double_t)i; | |
99 | } | |
100 | return binLim; | |
101 | } | |
102 | ||
103 | //_____________________________________________________________________________ | |
104 | TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins) | |
105 | { | |
106 | // | |
107 | // Make arbitrary binning, bins separated by a ',' | |
108 | // | |
109 | TString limits(bins); | |
110 | if (limits.IsNull()){ | |
111 | Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable"); | |
112 | return 0x0; | |
113 | } | |
114 | ||
115 | TObjArray *arr=limits.Tokenize(","); | |
116 | Int_t nLimits=arr->GetEntries(); | |
117 | if (nLimits<2){ | |
118 | Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable"); | |
119 | delete arr; | |
120 | return 0x0; | |
121 | } | |
122 | ||
123 | TVectorD *binLimits=new TVectorD(nLimits); | |
124 | for (Int_t iLim=0; iLim<nLimits; ++iLim){ | |
125 | (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof(); | |
126 | } | |
127 | ||
128 | delete arr; | |
129 | return binLimits; | |
130 | } | |
131 | ||
a94c2e7e | 132 | //_____________________________________________________________________________ |
133 | void AliDielectronHelper::GetMaxPtAndPhi(const AliVEvent *ev, Double_t &ptMax, Double_t &phiOfptMax){ | |
134 | // | |
135 | // find the highest pt and its phi in the event | |
136 | // | |
137 | for(Int_t itrk=0; itrk<ev->GetNumberOfTracks(); itrk++) { | |
138 | AliVParticle *part= ev->GetTrack(itrk); | |
139 | if(part && part->Pt() > ptMax) { | |
140 | ptMax = part->Pt(); | |
141 | phiOfptMax = part->Phi(); | |
142 | } | |
143 | } | |
144 | ||
145 | ||
146 | } | |
ba15fdfb | 147 | |
148 | //_____________________________________________________________________________ | |
149 | Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){ | |
150 | // determination of Nch | |
151 | if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1; | |
152 | ||
153 | AliStack *stack = ((AliMCEvent*)ev)->Stack(); | |
154 | ||
155 | if (!stack) return -1; | |
156 | ||
157 | Int_t nParticles = stack->GetNtrack(); | |
158 | Int_t nCh = 0; | |
159 | ||
160 | // count.. | |
161 | for (Int_t iMc = 0; iMc < nParticles; ++iMc) { | |
162 | if (!stack->IsPhysicalPrimary(iMc)) continue; | |
163 | ||
164 | TParticle* particle = stack->Particle(iMc); | |
165 | if (!particle) continue; | |
166 | if (particle->GetPDG()->Charge() == 0) continue; | |
167 | ||
168 | Float_t eta = particle->Eta(); | |
169 | if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++; | |
170 | } | |
171 | ||
172 | return nCh; | |
173 | } | |
174 | ||
175 | ||
176 | //_____________________________________________________________________________ | |
5720c765 | 177 | Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){ |
ba15fdfb | 178 | // Compute the collision multiplicity based on AOD or ESD tracklets |
179 | // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity() | |
180 | ||
181 | if (!ev) return -1; | |
182 | ||
183 | Int_t nTracklets = 0; | |
184 | Int_t nAcc = 0; | |
ba15fdfb | 185 | |
186 | if (ev->IsA() == AliAODEvent::Class()) { | |
187 | AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets(); | |
0c09cae4 | 188 | if (!tracklets) return -1; |
ba15fdfb | 189 | nTracklets = tracklets->GetNumberOfTracklets(); |
190 | for (Int_t nn = 0; nn < nTracklets; nn++) { | |
191 | Double_t theta = tracklets->GetTheta(nn); | |
192 | Double_t eta = -TMath::Log(TMath::Tan(theta/2.0)); | |
193 | if (TMath::Abs(eta) < etaRange) nAcc++; | |
194 | } | |
195 | } else if (ev->IsA() == AliESDEvent::Class()) { | |
196 | nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets(); | |
197 | for (Int_t nn = 0; nn < nTracklets; nn++) { | |
198 | Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn); | |
199 | if (TMath::Abs(eta) < etaRange) nAcc++; | |
200 | } | |
201 | } else return -1; | |
202 | ||
203 | return nAcc; | |
204 | } | |
205 | ||
206 | ||
5720c765 | 207 | |
208 | //________________________________________________________________ | |
209 | Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) { | |
210 | // | |
211 | // Correct the number of accepted tracklets based on the period and event vertex | |
212 | // | |
213 | ||
214 | Int_t runNo = event->GetRunNumber(); | |
215 | ||
216 | Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e | |
217 | Double_t refMult = 0.0; // reference multiplicity | |
218 | ||
219 | if(runNo>114930 && runNo<117223) period = 0; | |
220 | if(runNo>119158 && runNo<120830) period = 1; | |
221 | if(runNo>122373 && runNo<126438) period = 2; | |
222 | if(runNo>127711 && runNo<130841) period = 3; | |
223 | if(period<0 || period>3) return uncorrectedNacc; | |
224 | ||
225 | if(type<0 || type>8) return uncorrectedNacc; | |
226 | if(type == 0) refMult = 5.0; // SPD tracklets in |eta|<0.5 | |
227 | if(type == 1) refMult = 9.5; // SPD tracklets in |eta|<1.0 | |
228 | if(type == 2) refMult = 13.0; // SPD tracklets in |eta|<1.6 | |
229 | if(type == 3) refMult = 6.0; // ITSTPC+ in |eta|<0.5 | |
230 | if(type == 4) refMult = 12.0; // ITSTPC+ in |eta|<1.0 | |
231 | if(type == 5) refMult = 16.0; // ITSTPC+ in |eta|<1.6 | |
232 | if(type == 6) refMult = 6.0; // ITSSA+ in |eta|<0.5 | |
233 | if(type == 7) refMult = 12.0; // ITSSA+ in |eta|<1.0 | |
234 | if(type == 8) refMult = 15.0; // ITSSA+ in |eta|<1.6 | |
235 | ||
236 | if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc; | |
237 | ||
238 | TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type); | |
239 | if(!estimatorAvg) return uncorrectedNacc; | |
240 | ||
241 | Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ)); | |
242 | ||
243 | Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1); | |
244 | ||
245 | Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM)); | |
246 | ||
247 | return correctedNacc; | |
248 | } | |
249 | ||
ba15fdfb | 250 | //_____________________________________________________________________________ |
a823f01b | 251 | Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){ |
ba15fdfb | 252 | // put a robust Nacc definition here |
253 | ||
a823f01b | 254 | if (!ev) return -1; |
f2d9961b | 255 | |
256 | AliDielectronVarCuts varCuts; | |
257 | varCuts.AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0); | |
258 | varCuts.AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0); | |
259 | varCuts.AddCut(AliDielectronVarManager::kEta, -0.9, 0.9); | |
260 | varCuts.AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0); | |
261 | varCuts.AddCut(AliDielectronVarManager::kNclsTPC, 70.0, 160.0); | |
262 | varCuts.AddCut(AliDielectronVarManager::kKinkIndex0, -0.5, 0.5); //noKinks | |
263 | ||
264 | AliDielectronTrackCuts trkCuts; | |
265 | trkCuts.SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny); | |
266 | trkCuts.SetRequireITSRefit(kTRUE); | |
267 | trkCuts.SetRequireTPCRefit(kTRUE); | |
ba15fdfb | 268 | |
269 | Int_t nRecoTracks = ev->GetNumberOfTracks(); | |
270 | Int_t nAcc = 0; | |
f2d9961b | 271 | |
ba15fdfb | 272 | for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) { |
7c0a2a4d | 273 | AliVTrack *track = static_cast<AliVTrack*>(ev->GetTrack(iTrack)); |
a823f01b | 274 | if (!track) continue; |
f2d9961b | 275 | if (!trkCuts.IsSelected(track)) continue; |
276 | if (!varCuts.IsSelected(track)) continue; | |
41a64bfb | 277 | nAcc++; |
ba15fdfb | 278 | } |
a823f01b | 279 | |
ba15fdfb | 280 | return nAcc; |
281 | } | |
282 | ||
41a64bfb | 283 | //_____________________________________________________________________________ |
284 | Double_t AliDielectronHelper::GetITSTPCMatchEff(const AliVEvent *ev){ | |
285 | // recalulate the its-tpc matching efficiecy | |
286 | ||
287 | if (!ev) return -1; | |
288 | ||
f2d9961b | 289 | AliDielectronVarCuts varCutsTPC; |
290 | varCutsTPC.AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0); | |
291 | varCutsTPC.AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0); | |
292 | varCutsTPC.AddCut(AliDielectronVarManager::kEta, -0.9, 0.9); | |
293 | varCutsTPC.AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0); | |
294 | varCutsTPC.AddCut(AliDielectronVarManager::kNclsTPC, 50.0, 160.0); | |
295 | AliDielectronTrackCuts trkCutsTPC; | |
296 | trkCutsTPC.SetRequireTPCRefit(kTRUE); | |
41a64bfb | 297 | |
f2d9961b | 298 | AliDielectronVarCuts varCutsITS; |
299 | varCutsITS.AddCut(AliDielectronVarManager::kEta, -0.9, 0.9); | |
300 | AliDielectronTrackCuts trkCutsITS; | |
301 | trkCutsITS.SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny); | |
302 | trkCutsITS.SetRequireITSRefit(kTRUE); | |
41a64bfb | 303 | |
304 | ||
305 | Int_t nRecoTracks = ev->GetNumberOfTracks(); | |
306 | Double_t nTPC = 0, nITS = 0; | |
307 | ||
308 | for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) { | |
309 | AliVTrack *track = static_cast<AliVTrack*>(ev->GetTrack(iTrack)); | |
310 | if (!track) continue; | |
311 | ||
f2d9961b | 312 | if(!trkCutsITS.IsSelected(track)) continue; |
313 | if(!varCutsITS.IsSelected(track)) continue; | |
41a64bfb | 314 | nITS+=1.; |
315 | ||
f2d9961b | 316 | if(!trkCutsTPC.IsSelected(track)) continue; |
317 | if(!varCutsTPC.IsSelected(track)) continue; | |
41a64bfb | 318 | nTPC+=1.; |
319 | ||
320 | } | |
321 | ||
f2d9961b | 322 | // printf(" tracks TPC %.3e ITS %.3e = %.5f \n",nTPC,nITS,(nITS>0. ? nTPC/nITS : -1)); |
41a64bfb | 323 | return (nITS>0. ? nTPC/nITS : -1); |
324 | } | |
325 | ||
326 | ||
1201a1a9 | 327 | //_____________________________________________________________________________ |
328 | void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){ | |
329 | // Before rotate needs to be moved to position 0,0,0, ; move back after rotation | |
330 | if (!kfParticle) return; | |
331 | Double_t dx = 0.; | |
332 | Double_t dy = 0.; | |
333 | Double_t dz = 0.; | |
334 | ||
335 | if (ev){ | |
336 | dx = ev->GetPrimaryVertex()->GetX()-0.; | |
337 | dy = ev->GetPrimaryVertex()->GetY()-0.; | |
338 | dz = ev->GetPrimaryVertex()->GetZ()-0.; | |
339 | } | |
340 | ||
341 | kfParticle->X() = kfParticle->GetX() - dx; | |
342 | kfParticle->Y() = kfParticle->GetY() - dy; | |
343 | kfParticle->Z() = kfParticle->GetZ() - dz; | |
344 | ||
345 | ||
346 | // Rotate the kf particle | |
347 | Double_t c = cos(angle); | |
348 | Double_t s = sin(angle); | |
349 | ||
350 | Double_t mA[8][ 8]; | |
351 | for( Int_t i=0; i<8; i++ ){ | |
352 | for( Int_t j=0; j<8; j++){ | |
353 | mA[i][j] = 0; | |
354 | } | |
355 | } | |
356 | for( int i=0; i<8; i++ ){ | |
357 | mA[i][i] = 1; | |
358 | } | |
359 | mA[0][0] = c; mA[0][1] = s; | |
360 | mA[1][0] = -s; mA[1][1] = c; | |
361 | mA[3][3] = c; mA[3][4] = s; | |
362 | mA[4][3] = -s; mA[4][4] = c; | |
363 | ||
364 | Double_t mAC[8][8]; | |
365 | Double_t mAp[8]; | |
366 | ||
367 | for( Int_t i=0; i<8; i++ ){ | |
368 | mAp[i] = 0; | |
369 | for( Int_t k=0; k<8; k++){ | |
370 | mAp[i]+=mA[i][k] * kfParticle->GetParameter(k); | |
371 | } | |
372 | } | |
373 | ||
374 | for( Int_t i=0; i<8; i++){ | |
375 | kfParticle->Parameter(i) = mAp[i]; | |
376 | } | |
377 | ||
378 | for( Int_t i=0; i<8; i++ ){ | |
379 | for( Int_t j=0; j<8; j++ ){ | |
380 | mAC[i][j] = 0; | |
381 | for( Int_t k=0; k<8; k++ ){ | |
382 | mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j); | |
383 | } | |
384 | } | |
385 | } | |
386 | ||
387 | for( Int_t i=0; i<8; i++ ){ | |
388 | for( Int_t j=0; j<=i; j++ ){ | |
389 | Double_t xx = 0; | |
390 | for( Int_t k=0; k<8; k++){ | |
391 | xx+= mAC[i][k]*mA[j][k]; | |
392 | } | |
393 | kfParticle->Covariance(i,j) = xx; | |
394 | } | |
395 | } | |
396 | ||
397 | kfParticle->X() = kfParticle->GetX() + dx; | |
398 | kfParticle->Y() = kfParticle->GetY() + dy; | |
399 | kfParticle->Z() = kfParticle->GetZ() + dz; | |
400 | ||
401 | } | |
402 | ||
5720c765 | 403 | //_____________________________________________________________________________ |
404 | Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){ | |
40875e45 | 405 | // TODO: add AODs |
5720c765 | 406 | // counting number of mother particles generated in given eta range and 2 particle decay |
407 | if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1; | |
408 | ||
409 | AliStack *stack = ((AliMCEvent*)ev)->Stack(); | |
410 | ||
411 | if (!stack) return -1; | |
412 | ||
413 | Int_t nParticles = stack->GetNtrack(); | |
414 | Int_t nMothers = 0; | |
415 | ||
416 | // count.. | |
417 | for (Int_t iMc = 0; iMc < nParticles; ++iMc) { | |
418 | ||
419 | TParticle* particle = stack->Particle(iMc); | |
420 | if (!particle) continue; | |
421 | if (particle->GetPdgCode() != pdgMother) continue; | |
422 | if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue; | |
423 | ||
424 | if (particle->GetNDaughters() != 2) continue; | |
425 | // 1st daugther | |
426 | if (particle->GetFirstDaughter()>=nParticles || | |
427 | particle->GetFirstDaughter()<0 ) continue; | |
428 | ||
429 | TParticle* dau1 = stack->Particle(particle->GetFirstDaughter()); | |
430 | if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter) continue; | |
431 | if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue; | |
432 | ||
433 | // 2nd daughter | |
434 | if (particle->GetLastDaughter()>=nParticles || | |
435 | particle->GetLastDaughter()<0 ) continue; | |
436 | ||
437 | TParticle* dau2 = stack->Particle(particle->GetLastDaughter()); | |
438 | if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter) continue; | |
439 | if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue; | |
440 | ||
441 | // primary | |
442 | if (prim != -1) { | |
443 | if(particle->IsPrimary() != prim) continue; | |
444 | } | |
445 | nMothers++; | |
446 | } | |
447 | return nMothers; | |
448 | } | |
449 | ||
450 |