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