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.8  2000/12/21 16:24:06  morsch
  19 Coding convention clean-up
  20
  21 Revision 1.7  2000/10/02 15:16:37  morsch
  22 Correct coding rule violation for member data names of type fi -> fI.
  23
  24 Revision 1.6  2000/06/14 15:19:47  morsch
  25 Include clean-up (IH)
  26
  27 Revision 1.5  2000/06/09 20:35:32  morsch
  28 All coding rule violations except RS3 corrected
  29
  30 Revision 1.4  2000/03/20 18:03:24  morsch
  31 Change muon particle code to PDG code.
  32
  33 Revision 1.3  1999/09/29 09:24:08  fca
  34 Introduction of the Copyright and cvs Log
  35
  36 */
  37
  38 /*
  39  Class for dimuon analysis and fast dimuon simulation.
  40  It provides single and dimuon iterators, cuts, weighting, kinematic
  41  It uses the AliRun particle tree.
  42  Comments and suggestions to
  43  andreas.morsch@cern.ch
  44 */
  45
  46
  47 #include "AliDimuCombinator.h"
  48 #include "AliPDG.h"
  49 #include "AliRun.h"
  50 #include <TRandom.h>
  51 #include <TClonesArray.h>
  52 #include <TParticle.h>
  53 #include <TTree.h>
  54
  55 //
  56 ClassImp(AliDimuCombinator)
  57     AliDimuCombinator::AliDimuCombinator()
  58 {
  59 // Constructor
  60     fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
  61     fImuon1 = 0;
  62     fImuon2 = 0;
  63     fMuon1  = 0;
  64     fMuon2  = 0;
  65     fImin1  = 0;
  66     fImin2  = 0;
  67     fImax1  = fNParticle;
  68     fImax2  = fNParticle;
  69     fPtMin  = 0;
  70     fEtaMin = -10;
  71     fEtaMax = -10;
  72     fRate1  = 1.;
  73     fRate2  = 1.;
  74 }
  75
  76 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
  77 {
  78 // Dummy copy constructor
  79 }
  80
  81
  82 //
  83 //                       Iterators
  84 //
  85 TParticle* AliDimuCombinator::Particle(Int_t i)
  86 {
  87     return gAlice->Particle(i);
  88 }
  89
  90 TParticle* AliDimuCombinator::FirstMuon()
  91 {
  92 // Single muon iterator: initialisation
  93     fImuon1 = fImin1;
  94     fMuon1  = Particle(fImuon1);
  95     while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
  96         fImuon1++;
  97         if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
  98         fMuon1 = Particle(fImuon1);
  99     }
 100     return fMuon1;
 101 }
 102
 103 TParticle* AliDimuCombinator::FirstMuonSelected()
 104 {
 105 // Single selected muon iterator: initialisation
 106     TParticle* muon = FirstMuon();
 107     while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
 108     return muon;
 109 }
 110
 111
 112 TParticle* AliDimuCombinator::NextMuon()
 113 {
 114 // Single muon iterator: increment
 115     fImuon1++;
 116     if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
 117
 118     fMuon1 = Particle(fImuon1);
 119     while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
 120         fImuon1++;
 121         if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
 122         fMuon1 = Particle(fImuon1);
 123     }
 124     return fMuon1;
 125 }
 126
 127 TParticle* AliDimuCombinator::NextMuonSelected()
 128 {
 129 // Single selected muon iterator: increment
 130     TParticle * muon = NextMuon();
 131     while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
 132     return muon;
 133 }
 134
 135
 136 void AliDimuCombinator::FirstPartner()
 137 {
 138 // Helper for  dimuon iterator: initialisation
 139     if (fImin1 == fImin2) {
 140         fImuon2 = fImuon1+1;
 141     } else {
 142         fImuon2 = fImin2;
 143     }
 144     if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
 145     fMuon2 = Particle(fImuon2);
 146     while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
 147         fImuon2++;
 148         if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
 149         fMuon2 = Particle(fImuon2);
 150     }
 151 }
 152
 153 void AliDimuCombinator::FirstPartnerSelected()
 154 {
 155 // Helper for selected dimuon iterator: initialisation
 156     FirstPartner();
 157     while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
 158 }
 159
 160
 161 void AliDimuCombinator::NextPartner()
 162 {
 163 // Helper for dimuon iterator: increment
 164     fImuon2++;
 165     if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
 166
 167
 168     fMuon2 = Particle(fImuon2);
 169
 170     while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
 171         fImuon2++;
 172         if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
 173         fMuon2 = Particle(fImuon2);
 174     }
 175 }
 176
 177 void AliDimuCombinator::NextPartnerSelected()
 178 {
 179 // Helper for selected dimuon iterator: increment
 180     NextPartner();
 181     while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
 182 }
 183
 184
 185 TParticle*  AliDimuCombinator::Partner()
 186 {
 187 // Returns current partner for muon to form a dimuon
 188     return fMuon2;
 189 }
 190
 191 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
 192 {
 193 // Dimuon iterator: initialisation
 194     FirstMuon();
 195     FirstPartner();
 196     muon1 = fMuon1;
 197     muon2 = fMuon2;
 198 }
 199
 200 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
 201 {
 202 // Dimuon iterator: increment
 203     NextPartner();
 204     if (!Partner()) {
 205         NextMuon();
 206         FirstPartner();
 207     }
 208     muon1 = fMuon1;
 209     muon2 = fMuon2;
 210 }
 211 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
 212                                               TParticle* & muon2)
 213 {
 214 // Selected dimuon iterator: initialisation
 215     FirstMuonSelected();
 216     FirstPartnerSelected();
 217     muon1 = fMuon1;
 218     muon2 = fMuon2;
 219 }
 220
 221 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
 222                                              TParticle* & muon2)
 223 {
 224 // Selected dimuon iterator: increment
 225     NextPartnerSelected();
 226     if (!Partner()) {
 227         NextMuonSelected();
 228         FirstPartnerSelected();
 229     }
 230     muon1 = fMuon1;
 231     muon2 = fMuon2;
 232 }
 233
 234 void AliDimuCombinator::ResetRange()
 235 {
 236 // Reset index ranges for single muons
 237     fImin1 = fImin2 = 0;
 238     fImax1 = fImax2 = fNParticle;
 239 }
 240
 241 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
 242 {
 243 // Reset index range for first muon
 244     fImin1 = from;
 245     fImax1 = to;
 246     if (fImax1 > fNParticle) fImax1 = fNParticle;
 247 }
 248
 249 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
 250 {
 251 // Reset index range for second muon
 252     fImin2 = from;
 253     fImax2 = to;
 254     if (fImax2 > fNParticle) fImax2 = fNParticle;
 255 }
 256 //
 257 //                       Selection
 258 //
 259
 260 Bool_t AliDimuCombinator::Selected(TParticle* part)
 261 {
 262 // Selection cut for single muon
 263 //
 264     if (part == 0) {return 0;}
 265
 266     if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
 267         return 1;
 268     } else {
 269         return 0;
 270     }
 271 }
 272
 273 Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
 274 {
 275 // Selection cut for dimuons
 276 //
 277      return Selected(part1)*Selected(part2);
 278 }
 279 //
 280 //                       Kinematics
 281 //
 282 Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
 283 {
 284 // Invariant mass
 285 //
 286     Float_t px,py,pz,e;
 287     px = part1->Px()+part2->Px();
 288     py = part1->Py()+part2->Py();
 289     pz = part1->Pz()+part2->Pz();
 290     e  = part1->Energy()+part2->Energy();
 291     Float_t p = px*px+py*py+pz*pz;
 292     if (e*e < p) {
 293         return -1;
 294     } else {
 295         return TMath::Sqrt(e*e-p);
 296     }
 297 }
 298
 299 Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
 300 {
 301 // Transverse momentum of dimuons
 302 //
 303     Float_t px,py;
 304     px = part1->Px()+part2->Px();
 305     py = part1->Py()+part2->Py();
 306     return TMath::Sqrt(px*px+py*py);
 307 }
 308
 309 Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
 310 {
 311 // Pz of dimuon system
 312 //
 313     return part1->Pz()+part2->Pz();
 314 }
 315
 316 Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
 317 {
 318 // Rapidity of dimuon system
 319 //
 320     Float_t pz,e;
 321     pz = part1->Pz()+part2->Pz();
 322     e  = part1->Energy()+part2->Energy();
 323     return 0.5*TMath::Log((e+pz)/(e-pz));
 324 }
 325 //                  Response
 326 //
 327 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
 328 {
 329 // Apply gaussian smearing
 330 //
 331     value+=gRandom->Gaus(0, width);
 332 }
 333 //              Weighting
 334 //
 335
 336 Float_t AliDimuCombinator::DecayProbability(TParticle* part)
 337 {
 338 // Calculate decay probability for muons from pion and kaon decays
 339 //
 340     Float_t d, h, theta, cTau;
 341     TParticle* parent = Parent(part);
 342     Int_t ipar = Type(parent);
 343     if (ipar == kPiPlus || ipar == kPiMinus) {
 344         cTau=780.4;
 345     } else if (ipar == kKPlus || ipar == kKMinus) {
 346         cTau = 370.9;
 347     } else {
 348         cTau = 0;
 349     }
 350
 351
 352     Float_t gammaBeta=(parent->P())/(parent->GetMass());
 353 //
 354 // this part is still very ALICE muon-arm specific
 355 //
 356     theta = parent->Theta();
 357     h = 90*TMath::Tan(theta);
 358
 359     if (h<4) {
 360         d=4/TMath::Sin(theta);
 361     } else {
 362         d=90/TMath::Cos(theta);
 363     }
 364
 365     if (cTau > 0) {
 366         return 1-TMath::Exp(-d/cTau/gammaBeta);
 367     } else {
 368         return 1;
 369     }
 370 }
 371
 372 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
 373 {
 374 // Dimuon weight
 375
 376     Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
 377
 378     if (Correlated(part1, part2)) {
 379         return wgt/(Parent(part1)->GetWeight())*fRate1;
 380     } else {
 381         return wgt*fRate1*fRate2;
 382     }
 383 }
 384
 385
 386 Float_t AliDimuCombinator::Weight(TParticle* part)
 387 {
 388 // Single muon weight
 389     return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
 390 }
 391
 392 Bool_t  AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
 393 {
 394 // Check if muons are correlated
 395 //
 396     if (Origin(part1) == Origin(part2)) {
 397         return kTRUE;
 398     } else {
 399         return kFALSE;
 400     }
 401 }
 402
 403 TParticle* AliDimuCombinator::Parent(TParticle* part)
 404 {
 405 // Return pointer to parent
 406 //
 407     return Particle(part->GetFirstMother());
 408 }
 409
 410 Int_t AliDimuCombinator::Origin(TParticle* part)
 411 {
 412 // Return pointer to primary particle
 413 //
 414     Int_t iparent= part->GetFirstMother();
 415     if (iparent < 0) return iparent;
 416     Int_t ip;
 417     while(1) {
 418         ip = (Particle(iparent))->GetFirstMother();
 419         if (ip < 0) {
 420             break;
 421         } else {
 422             iparent = ip;
 423         }
 424     }
 425     return iparent;
 426 }
 427
 428 Int_t AliDimuCombinator::Type(TParticle *part)
 429 {
 430 // Return particle type for
 431 return part->GetPdgCode();
 432 }
 433
 434 AliDimuCombinator& AliDimuCombinator::operator=(const  AliDimuCombinator& rhs)
 435 {
 436 // Assignment operator
 437     return *this;
 438 }
 439
 440
 441 void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
 442 {
 443   //
 444   // Copy *this onto lego -- not implemented
 445   //
 446   Fatal("Copy","Not implemented!\n");
 447 }
 448
 449
 450
 451
 452