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