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