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