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