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
  39 AliDimuCombinator::AliDimuCombinator():
  40     fNParticle((Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries()),
  41     fImuon1(0),
  42     fImuon2(0),
  43     fImin1(0),
  44     fImin2(0),
  45     fImax1(fNParticle),
  46     fImax2(fNParticle),
  47     fRate1(1.),
  48     fRate2(1.),
  49     fMuon1(0),
  50     fMuon2(0),
  51     fPtMin(0.),
  52     fEtaMin(-10.),
  53     fEtaMax(10.)
  54 {
  55 // Constructor
  56     fNParticle = (Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries();
  57 }
  58
  59 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
  60     :TObject(combinator),
  61      fNParticle(0),
  62      fImuon1(0),
  63      fImuon2(0),
  64      fImin1(0),
  65      fImin2(0),
  66      fImax1(0),
  67      fImax2(0),
  68      fRate1(0),
  69      fRate2(0),
  70      fMuon1(0),
  71      fMuon2(0),
  72      fPtMin(0.),
  73      fEtaMin(0.),
  74      fEtaMax(0.)
  75 {
  76 // Dummy copy constructor
  77     combinator.Copy(*this);
  78 }
  79
  80
  81 //
  82 //                       Iterators
  83 //
  84 TParticle* AliDimuCombinator::Particle(Int_t i) const
  85 {
  86 //  Return next particle
  87 //
  88     return gAlice->GetMCApp()->Particle(i);
  89 }
  90
  91 TParticle* AliDimuCombinator::FirstMuon()
  92 {
  93 // Single muon iterator: initialisation
  94     fImuon1 = fImin1;
  95     fMuon1  = Particle(fImuon1);
  96     while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
  97         fImuon1++;
  98         if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
  99         fMuon1 = Particle(fImuon1);
 100     }
 101     return fMuon1;
 102 }
 103
 104 TParticle* AliDimuCombinator::FirstMuonSelected()
 105 {
 106 // Single selected muon iterator: initialisation
 107     TParticle* muon = FirstMuon();
 108     while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
 109     return muon;
 110 }
 111
 112
 113 TParticle* AliDimuCombinator::NextMuon()
 114 {
 115 // Single muon iterator: increment
 116     fImuon1++;
 117     if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
 118
 119     fMuon1 = Particle(fImuon1);
 120     while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
 121         fImuon1++;
 122         if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
 123         fMuon1 = Particle(fImuon1);
 124     }
 125     return fMuon1;
 126 }
 127
 128 TParticle* AliDimuCombinator::NextMuonSelected()
 129 {
 130 // Single selected muon iterator: increment
 131     TParticle * muon = NextMuon();
 132     while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
 133     return muon;
 134 }
 135
 136
 137 void AliDimuCombinator::FirstPartner()
 138 {
 139 // Helper for  dimuon iterator: initialisation
 140     if (fImin1 == fImin2) {
 141         fImuon2 = fImuon1+1;
 142     } else {
 143         fImuon2 = fImin2;
 144     }
 145     if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
 146     fMuon2 = Particle(fImuon2);
 147     while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
 148         fImuon2++;
 149         if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
 150         fMuon2 = Particle(fImuon2);
 151     }
 152 }
 153
 154 void AliDimuCombinator::FirstPartnerSelected()
 155 {
 156 // Helper for selected dimuon iterator: initialisation
 157     FirstPartner();
 158     while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
 159 }
 160
 161
 162 void AliDimuCombinator::NextPartner()
 163 {
 164 // Helper for dimuon iterator: increment
 165     fImuon2++;
 166     if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
 167
 168
 169     fMuon2 = Particle(fImuon2);
 170
 171     while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
 172         fImuon2++;
 173         if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
 174         fMuon2 = Particle(fImuon2);
 175     }
 176 }
 177
 178 void AliDimuCombinator::NextPartnerSelected()
 179 {
 180 // Helper for selected dimuon iterator: increment
 181     NextPartner();
 182     while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
 183 }
 184
 185
 186 TParticle*  AliDimuCombinator::Partner() const
 187 {
 188 // Returns current partner for muon to form a dimuon
 189     return fMuon2;
 190 }
 191
 192 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
 193 {
 194 // Dimuon iterator: initialisation
 195     FirstMuon();
 196     FirstPartner();
 197     muon1 = fMuon1;
 198     muon2 = fMuon2;
 199 }
 200
 201 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
 202 {
 203 // Dimuon iterator: increment
 204     NextPartner();
 205     if (!Partner()) {
 206         NextMuon();
 207         FirstPartner();
 208     }
 209     muon1 = fMuon1;
 210     muon2 = fMuon2;
 211 }
 212 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
 213                                               TParticle* & muon2)
 214 {
 215 // Selected dimuon iterator: initialisation
 216     FirstMuonSelected();
 217     FirstPartnerSelected();
 218     muon1 = fMuon1;
 219     muon2 = fMuon2;
 220 }
 221
 222 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
 223                                              TParticle* & muon2)
 224 {
 225 // Selected dimuon iterator: increment
 226     NextPartnerSelected();
 227     if (!Partner()) {
 228         NextMuonSelected();
 229         FirstPartnerSelected();
 230     }
 231     muon1 = fMuon1;
 232     muon2 = fMuon2;
 233 }
 234
 235 void AliDimuCombinator::ResetRange()
 236 {
 237 // Reset index ranges for single muons
 238     fImin1 = fImin2 = 0;
 239     fImax1 = fImax2 = fNParticle;
 240 }
 241
 242 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
 243 {
 244 // Reset index range for first muon
 245     fImin1 = from;
 246     fImax1 = to;
 247     if (fImax1 > fNParticle) fImax1 = fNParticle;
 248 }
 249
 250 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
 251 {
 252 // Reset index range for second muon
 253     fImin2 = from;
 254     fImax2 = to;
 255     if (fImax2 > fNParticle) fImax2 = fNParticle;
 256 }
 257 //
 258 //                       Selection
 259 //
 260
 261 Bool_t AliDimuCombinator::Selected(const TParticle* part) const
 262 {
 263 // Selection cut for single muon
 264 //
 265     if (part == 0) {return 0;}
 266
 267     if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
 268         return 1;
 269     } else {
 270         return 0;
 271     }
 272 }
 273
 274 Bool_t AliDimuCombinator::Selected(const TParticle* part1, const TParticle* part2) const
 275 {
 276 // Selection cut for dimuons
 277 //
 278      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     px = part1->Px()+part2->Px();
 289     py = part1->Py()+part2->Py();
 290     pz = part1->Pz()+part2->Pz();
 291     e  = part1->Energy()+part2->Energy();
 292     Float_t p = px*px+py*py+pz*pz;
 293     if (e*e < p) {
 294         return -1;
 295     } else {
 296         return TMath::Sqrt(e*e-p);
 297     }
 298 }
 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     px = part1->Px()+part2->Px();
 306     py = part1->Py()+part2->Py();
 307     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     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     pz = part1->Pz()+part2->Pz();
 323     e  = part1->Energy()+part2->Energy();
 324     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     value+=gRandom->Gaus(0, width);
 333 }
 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     TParticle* parent = Parent(part);
 344     Int_t ipar = Type(parent);
 345     if (ipar == kPiPlus || ipar == kPiMinus) {
 346         cTau=780.4;
 347     } else if (ipar == kKPlus || ipar == kKMinus) {
 348         cTau = 370.9;
 349     } else {
 350         cTau = 0;
 351     }
 352
 353
 354     Float_t gammaBeta=(parent->P())/(parent->GetMass());
 355 //
 356 // this part is still very ALICE muon-arm specific
 357 //
 358
 359
 360     theta = parent->Theta();
 361     h = 90*TMath::Tan(theta);
 362
 363     if (h<4) {
 364         d=4/TMath::Sin(theta);
 365     } else {
 366         d=90/TMath::Cos(theta);
 367     }
 368
 369     if (cTau > 0) {
 370         return 1-TMath::Exp(-d/cTau/gammaBeta);
 371     } else {
 372         return 1;
 373     }
 374 }
 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     Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
 391
 392     if (Correlated(part1, part2)) {
 393         if ( part1->GetFirstMother() == part2->GetFirstMother()) {
 394             return part1->GetWeight()*fRate1;
 395         } else {
 396             return wgt/(Parent(part1)->GetWeight())*fRate1;
 397         }
 398     } else {
 399         return wgt*fRate1*fRate2;
 400     }
 401 }
 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>&nbsp;
 411 <P>First -->
 412 <BR>The particle weight is given by w=R*M*Br
 413 <BR>&nbsp;with&nbsp; :
 414 <UL>R&nbsp;&nbsp; =&nbsp; 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&nbsp; 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&nbsp; = the rate of the mother production. This number depend on :
 421 <BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - 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>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the kinematic bias coming
 424 from the y and Pt cuts.&nbsp; Method&nbsp; 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>&nbsp;</UL>
 432 Next -->
 433 <BR>The weight of the dimuon depends on the correlation between muons
 434 <BR>&nbsp;
 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>&nbsp;</UL>
 449 */
 450 //End_Html
 451
 452
 453 Float_t AliDimuCombinator::Weight(const TParticle* part) const
 454 {
 455 // Single muon weight
 456     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     if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
 464
 465         return kTRUE;
 466     } else {
 467         return kFALSE;
 468     }
 469 }
 470
 471 TParticle* AliDimuCombinator::Parent(const TParticle* part) const
 472 {
 473 // Return pointer to parent
 474 //
 475     return Particle(part->GetFirstMother());
 476 }
 477
 478 Int_t AliDimuCombinator::Origin(const TParticle* part) const
 479 {
 480 // Return pointer to primary particle
 481 //
 482     Int_t iparent= part->GetFirstMother();
 483     if (iparent < 0) return iparent;
 484     Int_t ip;
 485     while(1) {
 486         ip = (Particle(iparent))->GetFirstMother();
 487         if (ip < 0) {
 488             break;
 489         } else {
 490             iparent = ip;
 491         }
 492     }
 493     return iparent;
 494 }
 495
 496 Int_t AliDimuCombinator::Type(const TParticle *part)  const
 497 {
 498 // Return particle type for
 499 return part->GetPdgCode();
 500 }
 501
 502 AliDimuCombinator& AliDimuCombinator::operator=(const  AliDimuCombinator& rhs)
 503 {
 504 // Assignment operator
 505     rhs.Copy(*this);
 506     return *this;
 507 }
 508
 509
 510 void AliDimuCombinator::Copy(TObject&) const
 511 {
 512   //
 513   // Copy
 514   //
 515   Fatal("Copy","Not implemented!\n");
 516 }
 517
 518
 519
 520
 521