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.5  2000/06/09 20:35:32  morsch
  19 All coding rule violations except RS3 corrected
  20
  21 Revision 1.4  2000/03/20 18:03:24  morsch
  22 Change muon particle code to PDG code.
  23
  24 Revision 1.3  1999/09/29 09:24:08  fca
  25 Introduction of the Copyright and cvs Log
  26
  27 */
  28
  29 //
  30 //
  31 //
  32 //
  33 #include "AliDimuCombinator.h"
  34 #include "AliPDG.h"
  35 #include <TRandom.h>
  36 #include <TClonesArray.h>
  37
  38 //
  39 ClassImp(AliDimuCombinator)
  40     AliDimuCombinator::AliDimuCombinator(TClonesArray* Partarray)
  41 {
  42 // Constructor
  43     fPartArray=Partarray;
  44     fNParticle=fPartArray->GetEntriesFast();
  45
  46     fimuon1 =0;
  47     fimuon2 =0;
  48     fmuon1  =0;
  49     fmuon2  =0;
  50     fimin1  = 0;
  51     fimin2  = 0;
  52     fimax1  = fNParticle;
  53     fimax2  = fNParticle;
  54     fPtMin  =0;
  55     fEtaMin =-10;
  56     fEtaMax =-10;
  57     fRate1=1.;
  58     fRate2=1.;
  59 }
  60
  61 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
  62 {
  63 // copy constructor
  64 }
  65
  66
  67 //
  68 //                       Iterators
  69 //
  70 TParticle* AliDimuCombinator::FirstMuon()
  71 {
  72 // Single muon iterator: initialisation
  73     fimuon1=fimin1;
  74     fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
  75     while(Type(fmuon1)!=kMuonPlus && Type(fmuon1)!=kMuonMinus) {
  76         fimuon1++;
  77         if (fimuon1 >= fimax1) {fmuon1=0; break;}
  78         fmuon1 = (TParticle*) fPartArray->UncheckedAt(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 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
  99     while(Type(fmuon1)!=kMuonPlus && Type(fmuon1)!=kMuonMinus) {
 100         fimuon1++;
 101         if (fimuon1>=fimax1) {fmuon1 = 0; break;}
 102         fmuon1 = (TParticle*) fPartArray->UncheckedAt(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 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
 126     while(Type(fmuon2)!=kMuonPlus && Type(fmuon2)!=kMuonMinus) {
 127         fimuon2++;
 128         if (fimuon2 >= fimax2) {fmuon2=0; break;}
 129         fmuon2 = (TParticle*) fPartArray->UncheckedAt(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 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
 149
 150     while(Type(fmuon2)!=kMuonPlus && Type(fmuon2)!=kMuonMinus) {
 151         fimuon2++;
 152         if (fimuon2>=fimax2) {fmuon2 = 0; break;}
 153         fmuon2 = (TParticle*) fPartArray->UncheckedAt(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     Float_t d, h, theta, cTau;
 321     TParticle* parent = Parent(part);
 322     Int_t ipar=Type(parent);
 323     if (ipar==kPiPlus || ipar==kPiMinus) {
 324         cTau=780.4;
 325     } else if (ipar==kKPlus || ipar==kKMinus) {
 326         cTau=370.9;
 327     } else {
 328         cTau=0;
 329     }
 330
 331
 332     Float_t gammaBeta=(parent->P())/(parent->GetMass());
 333 //
 334 // this part is still very ALICE muon-arm specific
 335 //
 336     theta=parent->Theta();
 337     h=90*TMath::Tan(theta);
 338
 339     if (h<4) {
 340         d=4/TMath::Sin(theta);
 341     } else {
 342         d=90/TMath::Cos(theta);
 343     }
 344
 345     if (cTau > 0) {
 346         return 1-TMath::Exp(-d/cTau/gammaBeta);
 347     } else {
 348         return 1;
 349     }
 350 }
 351
 352 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
 353 {
 354 // Dimuon weight
 355
 356     Float_t wgt=(part1->GetWeight())*(part2->GetWeight());
 357
 358     if (Correlated(part1, part2)) {
 359         return wgt/(Parent(part1)->GetWeight())*fRate1;
 360     } else {
 361         return wgt*fRate1*fRate2;
 362     }
 363 }
 364
 365
 366 Float_t AliDimuCombinator::Weight(TParticle* part)
 367 {
 368 // Single muon weight
 369     return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
 370 }
 371
 372 Bool_t  AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
 373 {
 374 // Check if muons are correlated
 375 //
 376     if (Origin(part1) == Origin(part2)) {
 377         return kTRUE;
 378     } else {
 379         return kFALSE;
 380     }
 381 }
 382
 383 TParticle* AliDimuCombinator::Parent(TParticle* part)
 384 {
 385 // Return pointer to parent
 386 //
 387     return (TParticle*) (fPartArray->UncheckedAt(part->GetFirstMother()));
 388 }
 389
 390 Int_t AliDimuCombinator::Origin(TParticle* part)
 391 {
 392 // Return pointer to primary particle
 393 //
 394     Int_t iparent= part->GetFirstMother();
 395     if (iparent < 0) return iparent;
 396     Int_t ip;
 397     while(1) {
 398         ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
 399         if (ip < 0) {
 400             break;
 401         } else {
 402             iparent=ip;
 403         }
 404     }
 405     return iparent;
 406 }
 407
 408 AliDimuCombinator& AliDimuCombinator::operator=(const  AliDimuCombinator& rhs)
 409 {
 410 // Assignment operator
 411     return *this;
 412 }
 413
 414
 415
 416
 417
 418
 419
 420