EVGEN/AliDimuCombinator.cxx

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