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