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