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