EVGEN/DimuonCombinator.cxx

   1 //
   2 //
   3 //
   4 //
   5 #include "DimuonCombinator.h"
   6 #include "AliRun.h"
   7 #include "TRandom.h"
   8 //
   9 ClassImp(DimuonCombinator)
  10 //
  11 //                       Iterators
  12 //
  13  GParticle* DimuonCombinator::FirstMuon()
  14      {
  15          fimuon1=fimin1;
  16          fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
  17          while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
  18              fimuon1++;
  19              if (fimuon1 >= fimax1) {fmuon1=0; break;}
  20              fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
  21          }
  22          return fmuon1;
  23      }
  24
  25  GParticle* DimuonCombinator::FirstMuonSelected()
  26      {
  27          GParticle * muon=FirstMuon();
  28          while(muon!=0 && !Selected(muon)) {muon=NextMuon();}
  29          return muon;
  30      }
  31
  32
  33  GParticle* DimuonCombinator::NextMuon()
  34      {
  35          fimuon1++;
  36          if (fimuon1>=fNParticle) {fmuon1 = 0; return fmuon1;}
  37
  38          fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
  39          while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
  40              fimuon1++;
  41              if (fimuon1>=fimax1) {fmuon1 = 0; break;}
  42              fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
  43          }
  44          return fmuon1;
  45      }
  46
  47 GParticle* DimuonCombinator::NextMuonSelected()
  48      {
  49          GParticle * muon=NextMuon();
  50          while(muon !=0 && !Selected(muon)) {muon=NextMuon();}
  51          return muon;
  52      }
  53
  54
  55  void DimuonCombinator::FirstPartner()
  56      {
  57          if (fimin1==fimin2) {
  58              fimuon2=fimuon1+1;
  59          } else {
  60              fimuon2=fimin2;
  61          }
  62          if (fimuon2 >= fimax2) {fmuon2=0; return;}
  63          fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
  64          while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
  65              fimuon2++;
  66              if (fimuon2 >= fimax2) {fmuon2=0; break;}
  67              fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
  68          }
  69      }
  70 void DimuonCombinator::FirstPartnerSelected()
  71 {
  72          FirstPartner();
  73          while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
  74 }
  75
  76
  77  void DimuonCombinator::NextPartner()
  78      {
  79          fimuon2++;
  80          if (fimuon2>=fimax2) {fmuon2 = 0; return;}
  81
  82
  83          fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
  84
  85          while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
  86              fimuon2++;
  87              if (fimuon2>=fimax2) {fmuon2 = 0; break;}
  88              fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
  89          }
  90
  91      }
  92
  93 void DimuonCombinator::NextPartnerSelected()
  94 {
  95          NextPartner();
  96          while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
  97 }
  98
  99
 100  GParticle*  DimuonCombinator::Partner()
 101      {
 102          return fmuon2;
 103      }
 104
 105 void DimuonCombinator::FirstMuonPair(GParticle* & muon1, GParticle* & muon2)
 106      {
 107          FirstMuon();
 108          FirstPartner();
 109          muon1=fmuon1;
 110          muon2=fmuon2;
 111      }
 112 void DimuonCombinator::NextMuonPair(GParticle* & muon1, GParticle* & muon2)
 113      {
 114          NextPartner();
 115          if (!Partner()) {
 116              NextMuon();
 117              FirstPartner();
 118          }
 119          muon1=fmuon1;
 120          muon2=fmuon2;
 121      }
 122 void DimuonCombinator::FirstMuonPairSelected(GParticle* & muon1, GParticle* & muon2)
 123      {
 124          FirstMuonSelected();
 125          FirstPartnerSelected();
 126          muon1=fmuon1;
 127          muon2=fmuon2;
 128      }
 129 void DimuonCombinator::NextMuonPairSelected(GParticle* & muon1, GParticle* & muon2)
 130      {
 131          NextPartnerSelected();
 132          if (!Partner()) {
 133              NextMuonSelected();
 134              FirstPartnerSelected();
 135          }
 136          muon1=fmuon1;
 137          muon2=fmuon2;
 138      }
 139 void DimuonCombinator::ResetRange()
 140 {
 141     fimin1=fimin2=0;
 142     fimax1=fimax2=fNParticle;
 143 }
 144
 145 void DimuonCombinator::SetFirstRange(Int_t from, Int_t to)
 146 {
 147     fimin1=from;
 148     fimax1=to;
 149     if (fimax1 > fNParticle) fimax1=fNParticle;
 150 }
 151
 152 void DimuonCombinator::SetSecondRange(Int_t from, Int_t to)
 153 {
 154     fimin2=from;
 155     fimax2=to;
 156     if (fimax2 > fNParticle) fimax2=fNParticle;
 157 }
 158 //
 159 //                       Selection
 160 //
 161
 162 Int_t DimuonCombinator::Selected(GParticle* part)
 163 {
 164 //
 165 //
 166     if (part==0) {return 0;}
 167
 168     if (part->GetPT() > fPtMin && part->GetEta()>fEtaMin && part->GetEta()<fEtaMax) {
 169         return 1;
 170     } else {
 171         return 0;
 172     }
 173
 174
 175 }
 176
 177 Int_t DimuonCombinator::Selected(GParticle* part1, GParticle* part2)
 178 {
 179      return Selected(part1)*Selected(part2);
 180 }
 181 //
 182 //                       Kinematics
 183 //
 184 Float_t DimuonCombinator::Mass(GParticle* part1, GParticle* part2)
 185 {
 186     Float_t px,py,pz,e;
 187     px=part1->GetPx()+part2->GetPx();
 188     py=part1->GetPy()+part2->GetPy();
 189     pz=part1->GetPz()+part2->GetPz();
 190     e =part1->GetEnergy()+part2->GetEnergy();
 191     Float_t p=px*px+py*py+pz*pz;
 192     if (e*e < p) {
 193         return -1;
 194     } else {
 195         return TMath::Sqrt(e*e-p);
 196     }
 197 }
 198
 199 Float_t DimuonCombinator::PT(GParticle* part1, GParticle* part2)
 200 {
 201     Float_t px,py;
 202     px=part1->GetPx()+part2->GetPx();
 203     py=part1->GetPy()+part2->GetPy();
 204     return TMath::Sqrt(px*px+py*py);
 205 }
 206
 207 Float_t DimuonCombinator::Pz(GParticle* part1, GParticle* part2)
 208 {
 209     return part1->GetPz()+part2->GetPz();
 210 }
 211
 212 Float_t DimuonCombinator::Y(GParticle* part1, GParticle* part2)
 213 {
 214     Float_t pz,e;
 215     pz=part1->GetPz()+part2->GetPz();
 216     e =part1->GetEnergy()+part2->GetEnergy();
 217     return 0.5*TMath::Log((e+pz)/(e-pz));
 218 }
 219 //                  Response
 220 //
 221 void DimuonCombinator::SmearGauss(Float_t width, Float_t & value)
 222 {
 223     value+=gRandom->Gaus(0, width);
 224 }
 225 //              Weighting
 226 //
 227
 228 Float_t DimuonCombinator::Weight(GParticle* part1, GParticle* part2)
 229 {
 230     if (part1->GetParent() == part2->GetParent()) {
 231         return (part1->GetWgt())*fRate1;
 232     } else {
 233         return (part1->GetWgt())*(part2->GetWgt())*fRate1*fRate2;
 234     }
 235 }
 236
 237 Float_t DimuonCombinator::Weight(GParticle* part)
 238 {
 239     return part->GetWgt()*fRate1;
 240 }
 241