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