]>
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> |
ffbede40 | 23 | // |
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 | ||
ba15fdfb | 132 | |
133 | //_____________________________________________________________________________ | |
134 | Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){ | |
135 | // determination of Nch | |
136 | if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1; | |
137 | ||
138 | AliStack *stack = ((AliMCEvent*)ev)->Stack(); | |
139 | ||
140 | if (!stack) return -1; | |
141 | ||
142 | Int_t nParticles = stack->GetNtrack(); | |
143 | Int_t nCh = 0; | |
144 | ||
145 | // count.. | |
146 | for (Int_t iMc = 0; iMc < nParticles; ++iMc) { | |
147 | if (!stack->IsPhysicalPrimary(iMc)) continue; | |
148 | ||
149 | TParticle* particle = stack->Particle(iMc); | |
150 | if (!particle) continue; | |
151 | if (particle->GetPDG()->Charge() == 0) continue; | |
152 | ||
153 | Float_t eta = particle->Eta(); | |
154 | if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++; | |
155 | } | |
156 | ||
157 | return nCh; | |
158 | } | |
159 | ||
160 | ||
161 | //_____________________________________________________________________________ | |
5720c765 | 162 | Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){ |
ba15fdfb | 163 | // Compute the collision multiplicity based on AOD or ESD tracklets |
164 | // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity() | |
165 | ||
166 | if (!ev) return -1; | |
167 | ||
168 | Int_t nTracklets = 0; | |
169 | Int_t nAcc = 0; | |
ba15fdfb | 170 | |
171 | if (ev->IsA() == AliAODEvent::Class()) { | |
172 | AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets(); | |
0c09cae4 | 173 | if (!tracklets) return -1; |
ba15fdfb | 174 | nTracklets = tracklets->GetNumberOfTracklets(); |
175 | for (Int_t nn = 0; nn < nTracklets; nn++) { | |
176 | Double_t theta = tracklets->GetTheta(nn); | |
177 | Double_t eta = -TMath::Log(TMath::Tan(theta/2.0)); | |
178 | if (TMath::Abs(eta) < etaRange) nAcc++; | |
179 | } | |
180 | } else if (ev->IsA() == AliESDEvent::Class()) { | |
181 | nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets(); | |
182 | for (Int_t nn = 0; nn < nTracklets; nn++) { | |
183 | Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn); | |
184 | if (TMath::Abs(eta) < etaRange) nAcc++; | |
185 | } | |
186 | } else return -1; | |
187 | ||
188 | return nAcc; | |
189 | } | |
190 | ||
191 | ||
5720c765 | 192 | |
193 | //________________________________________________________________ | |
194 | Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) { | |
195 | // | |
196 | // Correct the number of accepted tracklets based on the period and event vertex | |
197 | // | |
198 | ||
199 | Int_t runNo = event->GetRunNumber(); | |
200 | ||
201 | Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e | |
202 | Double_t refMult = 0.0; // reference multiplicity | |
203 | ||
204 | if(runNo>114930 && runNo<117223) period = 0; | |
205 | if(runNo>119158 && runNo<120830) period = 1; | |
206 | if(runNo>122373 && runNo<126438) period = 2; | |
207 | if(runNo>127711 && runNo<130841) period = 3; | |
208 | if(period<0 || period>3) return uncorrectedNacc; | |
209 | ||
210 | if(type<0 || type>8) return uncorrectedNacc; | |
211 | if(type == 0) refMult = 5.0; // SPD tracklets in |eta|<0.5 | |
212 | if(type == 1) refMult = 9.5; // SPD tracklets in |eta|<1.0 | |
213 | if(type == 2) refMult = 13.0; // SPD tracklets in |eta|<1.6 | |
214 | if(type == 3) refMult = 6.0; // ITSTPC+ in |eta|<0.5 | |
215 | if(type == 4) refMult = 12.0; // ITSTPC+ in |eta|<1.0 | |
216 | if(type == 5) refMult = 16.0; // ITSTPC+ in |eta|<1.6 | |
217 | if(type == 6) refMult = 6.0; // ITSSA+ in |eta|<0.5 | |
218 | if(type == 7) refMult = 12.0; // ITSSA+ in |eta|<1.0 | |
219 | if(type == 8) refMult = 15.0; // ITSSA+ in |eta|<1.6 | |
220 | ||
221 | if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc; | |
222 | ||
223 | TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type); | |
224 | if(!estimatorAvg) return uncorrectedNacc; | |
225 | ||
226 | Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ)); | |
227 | ||
228 | Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1); | |
229 | ||
230 | Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM)); | |
231 | ||
232 | return correctedNacc; | |
233 | } | |
234 | ||
ba15fdfb | 235 | //_____________________________________________________________________________ |
a823f01b | 236 | Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){ |
ba15fdfb | 237 | // put a robust Nacc definition here |
238 | ||
a823f01b | 239 | if (!ev) return -1; |
e4339752 | 240 | |
a823f01b | 241 | AliDielectronVarCuts *varCuts = new AliDielectronVarCuts("VarCuts","VarCuts"); |
242 | varCuts->AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0); | |
243 | varCuts->AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0); | |
244 | varCuts->AddCut(AliDielectronVarManager::kEta, -0.9, 0.9); | |
e4339752 | 245 | varCuts->AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0); |
a823f01b | 246 | varCuts->AddCut(AliDielectronVarManager::kNclsTPC, 70.0, 160.0); |
a823f01b | 247 | varCuts->AddCut(AliDielectronVarManager::kKinkIndex0, -0.5, 0.5); //noKinks |
248 | ||
249 | AliDielectronTrackCuts *trkCuts = new AliDielectronTrackCuts("TrkCuts","TrkCuts"); | |
e4339752 | 250 | trkCuts->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny); |
a823f01b | 251 | trkCuts->SetRequireITSRefit(kTRUE); |
252 | trkCuts->SetRequireTPCRefit(kTRUE); | |
ba15fdfb | 253 | |
254 | Int_t nRecoTracks = ev->GetNumberOfTracks(); | |
255 | Int_t nAcc = 0; | |
a823f01b | 256 | |
ba15fdfb | 257 | for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) { |
7c0a2a4d | 258 | AliVTrack *track = static_cast<AliVTrack*>(ev->GetTrack(iTrack)); |
a823f01b | 259 | if (!track) continue; |
260 | if (varCuts->IsSelected(track) && trkCuts->IsSelected(track)) | |
261 | nAcc++; | |
ba15fdfb | 262 | } |
a823f01b | 263 | |
264 | delete varCuts; | |
265 | delete trkCuts; | |
266 | ||
ba15fdfb | 267 | return nAcc; |
268 | } | |
269 | ||
1201a1a9 | 270 | //_____________________________________________________________________________ |
271 | void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){ | |
272 | // Before rotate needs to be moved to position 0,0,0, ; move back after rotation | |
273 | if (!kfParticle) return; | |
274 | Double_t dx = 0.; | |
275 | Double_t dy = 0.; | |
276 | Double_t dz = 0.; | |
277 | ||
278 | if (ev){ | |
279 | dx = ev->GetPrimaryVertex()->GetX()-0.; | |
280 | dy = ev->GetPrimaryVertex()->GetY()-0.; | |
281 | dz = ev->GetPrimaryVertex()->GetZ()-0.; | |
282 | } | |
283 | ||
284 | kfParticle->X() = kfParticle->GetX() - dx; | |
285 | kfParticle->Y() = kfParticle->GetY() - dy; | |
286 | kfParticle->Z() = kfParticle->GetZ() - dz; | |
287 | ||
288 | ||
289 | // Rotate the kf particle | |
290 | Double_t c = cos(angle); | |
291 | Double_t s = sin(angle); | |
292 | ||
293 | Double_t mA[8][ 8]; | |
294 | for( Int_t i=0; i<8; i++ ){ | |
295 | for( Int_t j=0; j<8; j++){ | |
296 | mA[i][j] = 0; | |
297 | } | |
298 | } | |
299 | for( int i=0; i<8; i++ ){ | |
300 | mA[i][i] = 1; | |
301 | } | |
302 | mA[0][0] = c; mA[0][1] = s; | |
303 | mA[1][0] = -s; mA[1][1] = c; | |
304 | mA[3][3] = c; mA[3][4] = s; | |
305 | mA[4][3] = -s; mA[4][4] = c; | |
306 | ||
307 | Double_t mAC[8][8]; | |
308 | Double_t mAp[8]; | |
309 | ||
310 | for( Int_t i=0; i<8; i++ ){ | |
311 | mAp[i] = 0; | |
312 | for( Int_t k=0; k<8; k++){ | |
313 | mAp[i]+=mA[i][k] * kfParticle->GetParameter(k); | |
314 | } | |
315 | } | |
316 | ||
317 | for( Int_t i=0; i<8; i++){ | |
318 | kfParticle->Parameter(i) = mAp[i]; | |
319 | } | |
320 | ||
321 | for( Int_t i=0; i<8; i++ ){ | |
322 | for( Int_t j=0; j<8; j++ ){ | |
323 | mAC[i][j] = 0; | |
324 | for( Int_t k=0; k<8; k++ ){ | |
325 | mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j); | |
326 | } | |
327 | } | |
328 | } | |
329 | ||
330 | for( Int_t i=0; i<8; i++ ){ | |
331 | for( Int_t j=0; j<=i; j++ ){ | |
332 | Double_t xx = 0; | |
333 | for( Int_t k=0; k<8; k++){ | |
334 | xx+= mAC[i][k]*mA[j][k]; | |
335 | } | |
336 | kfParticle->Covariance(i,j) = xx; | |
337 | } | |
338 | } | |
339 | ||
340 | kfParticle->X() = kfParticle->GetX() + dx; | |
341 | kfParticle->Y() = kfParticle->GetY() + dy; | |
342 | kfParticle->Z() = kfParticle->GetZ() + dz; | |
343 | ||
344 | } | |
345 | ||
5720c765 | 346 | //_____________________________________________________________________________ |
347 | Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){ | |
40875e45 | 348 | // TODO: add AODs |
5720c765 | 349 | // counting number of mother particles generated in given eta range and 2 particle decay |
350 | if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1; | |
351 | ||
352 | AliStack *stack = ((AliMCEvent*)ev)->Stack(); | |
353 | ||
354 | if (!stack) return -1; | |
355 | ||
356 | Int_t nParticles = stack->GetNtrack(); | |
357 | Int_t nMothers = 0; | |
358 | ||
359 | // count.. | |
360 | for (Int_t iMc = 0; iMc < nParticles; ++iMc) { | |
361 | ||
362 | TParticle* particle = stack->Particle(iMc); | |
363 | if (!particle) continue; | |
364 | if (particle->GetPdgCode() != pdgMother) continue; | |
365 | if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue; | |
366 | ||
367 | if (particle->GetNDaughters() != 2) continue; | |
368 | // 1st daugther | |
369 | if (particle->GetFirstDaughter()>=nParticles || | |
370 | particle->GetFirstDaughter()<0 ) continue; | |
371 | ||
372 | TParticle* dau1 = stack->Particle(particle->GetFirstDaughter()); | |
373 | if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter) continue; | |
374 | if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue; | |
375 | ||
376 | // 2nd daughter | |
377 | if (particle->GetLastDaughter()>=nParticles || | |
378 | particle->GetLastDaughter()<0 ) continue; | |
379 | ||
380 | TParticle* dau2 = stack->Particle(particle->GetLastDaughter()); | |
381 | if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter) continue; | |
382 | if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue; | |
383 | ||
384 | // primary | |
385 | if (prim != -1) { | |
386 | if(particle->IsPrimary() != prim) continue; | |
387 | } | |
388 | nMothers++; | |
389 | } | |
390 | return nMothers; | |
391 | } | |
392 | ||
393 |