Line data Source code
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 : /* $Id$ */
17 :
18 : //
19 : // Class for dimuon analysis and fast dimuon simulation.
20 : // It provides single and dimuon iterators, cuts, weighting, kinematic
21 : // It uses the AliRun particle tree.
22 : // Comments and suggestions to
23 : // andreas.morsch@cern.ch
24 :
25 :
26 : #include <TClonesArray.h>
27 : #include <TParticle.h>
28 : #include <TPDGCode.h>
29 : #include <TRandom.h>
30 : #include <TTree.h>
31 :
32 : #include "AliDimuCombinator.h"
33 : #include "AliRun.h"
34 : #include "AliMC.h"
35 :
36 : //
37 6 : ClassImp(AliDimuCombinator)
38 :
39 0 : AliDimuCombinator::AliDimuCombinator():
40 0 : fNParticle((Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries()),
41 0 : fImuon1(0),
42 0 : fImuon2(0),
43 0 : fImin1(0),
44 0 : fImin2(0),
45 0 : fImax1(fNParticle),
46 0 : fImax2(fNParticle),
47 0 : fRate1(1.),
48 0 : fRate2(1.),
49 0 : fMuon1(0),
50 0 : fMuon2(0),
51 0 : fPtMin(0.),
52 0 : fEtaMin(-10.),
53 0 : fEtaMax(10.)
54 0 : {
55 : // Constructor
56 0 : fNParticle = (Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries();
57 0 : }
58 :
59 : AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
60 0 : :TObject(combinator),
61 0 : fNParticle(0),
62 0 : fImuon1(0),
63 0 : fImuon2(0),
64 0 : fImin1(0),
65 0 : fImin2(0),
66 0 : fImax1(0),
67 0 : fImax2(0),
68 0 : fRate1(0),
69 0 : fRate2(0),
70 0 : fMuon1(0),
71 0 : fMuon2(0),
72 0 : fPtMin(0.),
73 0 : fEtaMin(0.),
74 0 : fEtaMax(0.)
75 0 : {
76 : // Dummy copy constructor
77 0 : combinator.Copy(*this);
78 0 : }
79 :
80 :
81 : //
82 : // Iterators
83 : //
84 : TParticle* AliDimuCombinator::Particle(Int_t i) const
85 : {
86 : // Return next particle
87 : //
88 0 : return gAlice->GetMCApp()->Particle(i);
89 : }
90 :
91 : TParticle* AliDimuCombinator::FirstMuon()
92 : {
93 : // Single muon iterator: initialisation
94 0 : fImuon1 = fImin1;
95 0 : fMuon1 = Particle(fImuon1);
96 0 : while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
97 0 : fImuon1++;
98 0 : if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
99 0 : fMuon1 = Particle(fImuon1);
100 : }
101 0 : return fMuon1;
102 : }
103 :
104 : TParticle* AliDimuCombinator::FirstMuonSelected()
105 : {
106 : // Single selected muon iterator: initialisation
107 0 : TParticle* muon = FirstMuon();
108 0 : while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
109 0 : return muon;
110 : }
111 :
112 :
113 : TParticle* AliDimuCombinator::NextMuon()
114 : {
115 : // Single muon iterator: increment
116 0 : fImuon1++;
117 0 : if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
118 :
119 0 : fMuon1 = Particle(fImuon1);
120 0 : while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
121 0 : fImuon1++;
122 0 : if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
123 0 : fMuon1 = Particle(fImuon1);
124 : }
125 0 : return fMuon1;
126 0 : }
127 :
128 : TParticle* AliDimuCombinator::NextMuonSelected()
129 : {
130 : // Single selected muon iterator: increment
131 0 : TParticle * muon = NextMuon();
132 0 : while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
133 0 : return muon;
134 : }
135 :
136 :
137 : void AliDimuCombinator::FirstPartner()
138 : {
139 : // Helper for dimuon iterator: initialisation
140 0 : if (fImin1 == fImin2) {
141 0 : fImuon2 = fImuon1+1;
142 0 : } else {
143 0 : fImuon2 = fImin2;
144 : }
145 0 : if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
146 0 : fMuon2 = Particle(fImuon2);
147 0 : while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
148 0 : fImuon2++;
149 0 : if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
150 0 : fMuon2 = Particle(fImuon2);
151 : }
152 0 : }
153 :
154 : void AliDimuCombinator::FirstPartnerSelected()
155 : {
156 : // Helper for selected dimuon iterator: initialisation
157 0 : FirstPartner();
158 0 : while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
159 0 : }
160 :
161 :
162 : void AliDimuCombinator::NextPartner()
163 : {
164 : // Helper for dimuon iterator: increment
165 0 : fImuon2++;
166 0 : if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
167 :
168 :
169 0 : fMuon2 = Particle(fImuon2);
170 :
171 0 : while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
172 0 : fImuon2++;
173 0 : if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
174 0 : fMuon2 = Particle(fImuon2);
175 : }
176 0 : }
177 :
178 : void AliDimuCombinator::NextPartnerSelected()
179 : {
180 : // Helper for selected dimuon iterator: increment
181 0 : NextPartner();
182 0 : while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
183 0 : }
184 :
185 :
186 : TParticle* AliDimuCombinator::Partner() const
187 : {
188 : // Returns current partner for muon to form a dimuon
189 0 : return fMuon2;
190 : }
191 :
192 : void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
193 : {
194 : // Dimuon iterator: initialisation
195 0 : FirstMuon();
196 0 : FirstPartner();
197 0 : muon1 = fMuon1;
198 0 : muon2 = fMuon2;
199 0 : }
200 :
201 : void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
202 : {
203 : // Dimuon iterator: increment
204 0 : NextPartner();
205 0 : if (!Partner()) {
206 0 : NextMuon();
207 0 : FirstPartner();
208 0 : }
209 0 : muon1 = fMuon1;
210 0 : muon2 = fMuon2;
211 0 : }
212 : void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
213 : TParticle* & muon2)
214 : {
215 : // Selected dimuon iterator: initialisation
216 0 : FirstMuonSelected();
217 0 : FirstPartnerSelected();
218 0 : muon1 = fMuon1;
219 0 : muon2 = fMuon2;
220 0 : }
221 :
222 : void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
223 : TParticle* & muon2)
224 : {
225 : // Selected dimuon iterator: increment
226 0 : NextPartnerSelected();
227 0 : if (!Partner()) {
228 0 : NextMuonSelected();
229 0 : FirstPartnerSelected();
230 0 : }
231 0 : muon1 = fMuon1;
232 0 : muon2 = fMuon2;
233 0 : }
234 :
235 : void AliDimuCombinator::ResetRange()
236 : {
237 : // Reset index ranges for single muons
238 0 : fImin1 = fImin2 = 0;
239 0 : fImax1 = fImax2 = fNParticle;
240 0 : }
241 :
242 : void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
243 : {
244 : // Reset index range for first muon
245 0 : fImin1 = from;
246 0 : fImax1 = to;
247 0 : if (fImax1 > fNParticle) fImax1 = fNParticle;
248 0 : }
249 :
250 : void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
251 : {
252 : // Reset index range for second muon
253 0 : fImin2 = from;
254 0 : fImax2 = to;
255 0 : if (fImax2 > fNParticle) fImax2 = fNParticle;
256 0 : }
257 : //
258 : // Selection
259 : //
260 :
261 : Bool_t AliDimuCombinator::Selected(const TParticle* part) const
262 : {
263 : // Selection cut for single muon
264 : //
265 0 : if (part == 0) {return 0;}
266 :
267 0 : if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
268 0 : return 1;
269 : } else {
270 0 : return 0;
271 : }
272 0 : }
273 :
274 : Bool_t AliDimuCombinator::Selected(const TParticle* part1, const TParticle* part2) const
275 : {
276 : // Selection cut for dimuons
277 : //
278 0 : return Selected(part1)*Selected(part2);
279 : }
280 : //
281 : // Kinematics
282 : //
283 : Float_t AliDimuCombinator::Mass(const TParticle* part1, const TParticle* part2) const
284 : {
285 : // Invariant mass
286 : //
287 : Float_t px,py,pz,e;
288 0 : px = part1->Px()+part2->Px();
289 0 : py = part1->Py()+part2->Py();
290 0 : pz = part1->Pz()+part2->Pz();
291 0 : e = part1->Energy()+part2->Energy();
292 0 : Float_t p = px*px+py*py+pz*pz;
293 0 : if (e*e < p) {
294 0 : return -1;
295 : } else {
296 0 : return TMath::Sqrt(e*e-p);
297 : }
298 0 : }
299 :
300 : Float_t AliDimuCombinator::PT(const TParticle* part1, const TParticle* part2) const
301 : {
302 : // Transverse momentum of dimuons
303 : //
304 : Float_t px,py;
305 0 : px = part1->Px()+part2->Px();
306 0 : py = part1->Py()+part2->Py();
307 0 : return TMath::Sqrt(px*px+py*py);
308 : }
309 :
310 : Float_t AliDimuCombinator::Pz(const TParticle* part1, const TParticle* part2) const
311 : {
312 : // Pz of dimuon system
313 : //
314 0 : return part1->Pz()+part2->Pz();
315 : }
316 :
317 : Float_t AliDimuCombinator::Y(const TParticle* part1, const TParticle* part2) const
318 : {
319 : // Rapidity of dimuon system
320 : //
321 : Float_t pz,e;
322 0 : pz = part1->Pz()+part2->Pz();
323 0 : e = part1->Energy()+part2->Energy();
324 0 : return 0.5*TMath::Log((e+pz)/(e-pz));
325 : }
326 : // Response
327 : //
328 : void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const
329 : {
330 : // Apply gaussian smearing
331 : //
332 0 : value+=gRandom->Gaus(0, width);
333 0 : }
334 : // Weighting
335 : //
336 :
337 : Float_t AliDimuCombinator::DecayProbability(const TParticle* part) const
338 : {
339 : // Calculate decay probability for muons from pion and kaon decays
340 : //
341 :
342 : Float_t d, h, theta, cTau;
343 0 : TParticle* parent = Parent(part);
344 0 : Int_t ipar = Type(parent);
345 0 : if (ipar == kPiPlus || ipar == kPiMinus) {
346 : cTau=780.4;
347 0 : } else if (ipar == kKPlus || ipar == kKMinus) {
348 : cTau = 370.9;
349 0 : } else {
350 : cTau = 0;
351 : }
352 :
353 :
354 0 : Float_t gammaBeta=(parent->P())/(parent->GetMass());
355 : //
356 : // this part is still very ALICE muon-arm specific
357 : //
358 :
359 :
360 0 : theta = parent->Theta();
361 0 : h = 90*TMath::Tan(theta);
362 :
363 0 : if (h<4) {
364 0 : d=4/TMath::Sin(theta);
365 0 : } else {
366 0 : d=90/TMath::Cos(theta);
367 : }
368 :
369 0 : if (cTau > 0) {
370 0 : return 1-TMath::Exp(-d/cTau/gammaBeta);
371 : } else {
372 0 : return 1;
373 : }
374 0 : }
375 :
376 : //Begin_Html
377 : /*
378 : <p> In the the code above :
379 : <P>If h is less than 4 cm, pions or kaons go in the beam pipe and can have a long way
380 : <BR>If h is greater than 4 cm, pions or kaons crash into the front absorber
381 : <P><IMG SRC="absorbeur.jpg" HEIGHT=292 WIDTH=819>
382 : */
383 : //End_Html
384 :
385 :
386 : Float_t AliDimuCombinator::Weight(const TParticle* part1, const TParticle* part2) const
387 : {
388 : // Dimuon weight
389 :
390 0 : Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
391 :
392 0 : if (Correlated(part1, part2)) {
393 0 : if ( part1->GetFirstMother() == part2->GetFirstMother()) {
394 0 : return part1->GetWeight()*fRate1;
395 : } else {
396 0 : return wgt/(Parent(part1)->GetWeight())*fRate1;
397 : }
398 : } else {
399 0 : return wgt*fRate1*fRate2;
400 : }
401 0 : }
402 :
403 : //Begin_Html
404 : /*
405 : <p>Some clarifications on the calculation of the dimuons weight :
406 : <P>We must keep in mind that if we force the meson decay in muons and we put
407 : lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
408 : obliged to calculate different weights to correct the number
409 : of muons
410 : <BR>
411 : <P>First -->
412 : <BR>The particle weight is given by w=R*M*Br
413 : <BR> with :
414 : <UL>R = the rate by event. This number gives the number
415 : of produced J/psi, upsilon, pion ... in a collision.
416 : <BR>It corresponds of the weight 0.06 given for example in gener->AddGenerator(jpsi,"J/Psi",
417 : 0.06); from the config.C macro.
418 : <BR>In this example R=0.06
419 :
420 : <P>M = the rate of the mother production. This number depend on :
421 : <BR> - the number of generated events --> fParentWeight=1./Float_t(fNpart) in AliGenPythia.cxx . This
422 : is a normalization to 1 of the number of generated particles.
423 : <BR> - the kinematic bias coming
424 : from the y and Pt cuts. Method AliGenPythia::AdjustWeights() in AliGenPythia.cxx
425 : <BR>(in AliGenParam.cxx this 2 things are taken into account in fParentWeight
426 : = fYWgt*fPtWgt*phiWgt/fNpart )
427 :
428 : <P>Br = the branching ratio in muon from the mother decay</UL>
429 :
430 : <P><BR>In this method, part->GetWeight() = M*Br
431 : <UL> </UL>
432 : Next -->
433 : <BR>The weight of the dimuon depends on the correlation between muons
434 : <BR>
435 : <UL>If the muons are correlated and come from a resonance (for example
436 : J/psi -> mu+ mu-) , the weight of the dimuon is the weight of one muon then
437 : <BR>w12= R*M*Br = w1* R1 (in this method this gives part1->GetWeight()*fRate1)
438 :
439 : <P>If the muons are correlated and come from a charm or a bottom pair then
440 : w12 = M*R*Br1*Br2 = w1*w2*R1/M1
441 : <BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
442 : Indeed the 2 muons come from the same mother so the
443 : <BR>weight of a DD~ or BB~ is M*Br and they are no correlation in the decay
444 : (Br1*Br2)
445 :
446 : <P>If the muons are not correlated w12 = M1*M2*R1*R2*Br1*Br2 = w1*w2*R1*R2
447 : (in this method this gives wgt*fRate1*fRate2)
448 : <BR> </UL>
449 : */
450 : //End_Html
451 :
452 :
453 : Float_t AliDimuCombinator::Weight(const TParticle* part) const
454 : {
455 : // Single muon weight
456 0 : return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
457 : }
458 :
459 : Bool_t AliDimuCombinator::Correlated(const TParticle* part1, const TParticle* part2) const
460 : {
461 : // Check if muons are correlated
462 : //
463 0 : if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
464 :
465 0 : return kTRUE;
466 : } else {
467 0 : return kFALSE;
468 : }
469 0 : }
470 :
471 : TParticle* AliDimuCombinator::Parent(const TParticle* part) const
472 : {
473 : // Return pointer to parent
474 : //
475 0 : return Particle(part->GetFirstMother());
476 : }
477 :
478 : Int_t AliDimuCombinator::Origin(const TParticle* part) const
479 : {
480 : // Return pointer to primary particle
481 : //
482 0 : Int_t iparent= part->GetFirstMother();
483 0 : if (iparent < 0) return iparent;
484 : Int_t ip;
485 0 : while(1) {
486 0 : ip = (Particle(iparent))->GetFirstMother();
487 0 : if (ip < 0) {
488 : break;
489 : } else {
490 : iparent = ip;
491 : }
492 : }
493 : return iparent;
494 0 : }
495 :
496 : Int_t AliDimuCombinator::Type(const TParticle *part) const
497 : {
498 : // Return particle type for
499 0 : return part->GetPdgCode();
500 : }
501 :
502 : AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
503 : {
504 : // Assignment operator
505 0 : rhs.Copy(*this);
506 0 : return *this;
507 : }
508 :
509 :
510 : void AliDimuCombinator::Copy(TObject&) const
511 : {
512 : //
513 : // Copy
514 : //
515 0 : Fatal("Copy","Not implemented!\n");
516 0 : }
517 :
518 :
519 :
520 :
521 :
|
