]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliDimuCombinator.cxx
Allowing modularity of the MUON geometry during the generation (geant) phase with...
[u/mrichter/AliRoot.git] / EVGEN / AliDimuCombinator.cxx
index f439fd5c01bd30f22b471dfce090e3c8a4c358a3..44001c878c2fb60597eca10af51270422d779960 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.8  2000/12/21 16:24:06  morsch
-Coding convention clean-up
-
-Revision 1.7  2000/10/02 15:16:37  morsch
-Correct coding rule violation for member data names of type fi -> fI.
-
-Revision 1.6  2000/06/14 15:19:47  morsch
-Include clean-up (IH)
-
-Revision 1.5  2000/06/09 20:35:32  morsch
-All coding rule violations except RS3 corrected
-
-Revision 1.4  2000/03/20 18:03:24  morsch
-Change muon particle code to PDG code.
-
-Revision 1.3  1999/09/29 09:24:08  fca
-Introduction of the Copyright and cvs Log
-
-*/
+/* $Id$ */
 
-/*  
- Class for dimuon analysis and fast dimuon simulation.
- It provides single and dimuon iterators, cuts, weighting, kinematic
- It uses the AliRun particle tree.
- Comments and suggestions to 
- andreas.morsch@cern.ch
-*/
+// 
+// Class for dimuon analysis and fast dimuon simulation.
+// It provides single and dimuon iterators, cuts, weighting, kinematic
+// It uses the AliRun particle tree.
+// Comments and suggestions to 
+// andreas.morsch@cern.ch
 
 
-#include "AliDimuCombinator.h" 
-#include "AliPDG.h" 
-#include "AliRun.h" 
-#include <TRandom.h>
 #include <TClonesArray.h>
 #include <TParticle.h>
+#include <TPDGCode.h> 
+#include <TRandom.h>
 #include <TTree.h>
 
+#include "AliDimuCombinator.h" 
+#include "AliRun.h" 
+#include "AliMC.h"
+
 //
 ClassImp(AliDimuCombinator)
     AliDimuCombinator::AliDimuCombinator() 
@@ -74,17 +55,21 @@ ClassImp(AliDimuCombinator)
 }
 
 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
+    :TObject(combinator)
 {
 // Dummy copy constructor
+    combinator.Copy(*this);
 }
 
 
 //
 //                       Iterators
 // 
-TParticle* AliDimuCombinator::Particle(Int_t i)
+TParticle* AliDimuCombinator::Particle(Int_t i) const
 {
-    return gAlice->Particle(i);
+//  Return next particle
+//
+    return gAlice->GetMCApp()->Particle(i);
 }
 
 TParticle* AliDimuCombinator::FirstMuon()
@@ -182,7 +167,7 @@ void AliDimuCombinator::NextPartnerSelected()
 }
 
 
-TParticle*  AliDimuCombinator::Partner()
+TParticle*  AliDimuCombinator::Partner() const
 {
 // Returns current partner for muon to form a dimuon
     return fMuon2;
@@ -257,7 +242,7 @@ void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
 //                       Selection
 //
 
-Bool_t AliDimuCombinator::Selected(TParticle* part)
+Bool_t AliDimuCombinator::Selected(TParticle* part) const
 {
 // Selection cut for single muon 
 //
@@ -270,7 +255,7 @@ Bool_t AliDimuCombinator::Selected(TParticle* part)
     }
 }
 
-Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
+Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2) const
 {
 // Selection cut for dimuons
 //
@@ -279,7 +264,7 @@ Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
 //
 //                       Kinematics
 //
-Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
+Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2) const
 {
 // Invariant mass
 //
@@ -296,7 +281,7 @@ Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
     }
 }
 
-Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
+Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2) const
 {
 // Transverse momentum of dimuons
 //
@@ -306,14 +291,14 @@ Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
     return TMath::Sqrt(px*px+py*py);
 }
 
-Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
+Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2) const
 {
 // Pz of dimuon system
 //
     return part1->Pz()+part2->Pz();
 }
 
-Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
+Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2) const
 {
 // Rapidity of dimuon system
 //
@@ -324,7 +309,7 @@ Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
 }
 //                  Response
 //
-void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
+void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const
 {
 // Apply gaussian smearing
 //
@@ -333,10 +318,11 @@ void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
 //              Weighting
 // 
 
-Float_t AliDimuCombinator::DecayProbability(TParticle* part)
+Float_t AliDimuCombinator::DecayProbability(TParticle* part) const
 {
 // Calculate decay probability for muons from pion and kaon decays
 // 
+
     Float_t d, h, theta, cTau;
     TParticle* parent = Parent(part);
     Int_t ipar = Type(parent);
@@ -353,6 +339,8 @@ Float_t AliDimuCombinator::DecayProbability(TParticle* part)
 //
 // this part is still very ALICE muon-arm specific
 //
+
+
     theta = parent->Theta();
     h = 90*TMath::Tan(theta);
     
@@ -369,45 +357,109 @@ Float_t AliDimuCombinator::DecayProbability(TParticle* part)
     }
 }
 
-Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
+//Begin_Html
+/*
+<p> In the the code above :
+<P>If h is less than 4 cm, pions or kaons go in the beam pipe and can have a long way
+<BR>If h is greater than 4 cm, pions or kaons crash into the front absorber
+<P><IMG SRC="absorbeur.jpg" HEIGHT=292 WIDTH=819>
+*/
+//End_Html
+
+
+Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2) const
 {
 // Dimuon weight
 
     Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
     
     if (Correlated(part1, part2)) {
-       return wgt/(Parent(part1)->GetWeight())*fRate1;
+       if ( part1->GetFirstMother() == part2->GetFirstMother()) {
+           return part1->GetWeight()*fRate1;
+       } else {
+           return wgt/(Parent(part1)->GetWeight())*fRate1;
+       }
     } else {
        return wgt*fRate1*fRate2;
     }
 } 
 
+//Begin_Html
+/*
+<p>Some clarifications on the calculation of the dimuons weight :
+<P>We must keep in mind that if we force the meson decay in muons and we put
+lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
+obliged to calculate different weights to correct the number
+of muons
+<BR>&nbsp;
+<P>First -->
+<BR>The particle weight is given by w=R*M*Br
+<BR>&nbsp;with&nbsp; :
+<UL>R&nbsp;&nbsp; =&nbsp; the rate by event. This number gives the number
+of produced J/psi, upsilon, pion ... in a collision.
+<BR>It corresponds of the weight 0.06 given for example in&nbsp; gener->AddGenerator(jpsi,"J/Psi",
+0.06); from the config.C macro.
+<BR>In this example R=0.06
+
+<P>M&nbsp; = the rate of the mother production. This number depend on :
+<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the number of generated events --> fParentWeight=1./Float_t(fNpart) in AliGenPythia.cxx . This
+is a normalization to 1 of the number of generated particles.
+<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the kinematic bias coming
+from the y and Pt cuts.&nbsp; Method&nbsp; AliGenPythia::AdjustWeights() in AliGenPythia.cxx
+<BR>(in AliGenParam.cxx this 2 things are taken into account in fParentWeight
+= fYWgt*fPtWgt*phiWgt/fNpart )
+
+<P>Br = the branching ratio in muon from the mother decay</UL>
+
+<P><BR>In this method, part->GetWeight() = M*Br
+<UL>&nbsp;</UL>
+Next -->
+<BR>The weight of the dimuon depends on the correlation between muons
+<BR>&nbsp;
+<UL>If the muons are correlated and come from a resonance (for example
+J/psi -> mu+ mu-) , the weight of the dimuon is the weight of one muon then
+<BR>w12= R*M*Br = w1* R1 (in this method this gives part1->GetWeight()*fRate1)
+
+<P>If the muons are correlated and come from a charm or a bottom pair then
+w12 = M*R*Br1*Br2 = w1*w2*R1/M1
+<BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
+Indeed the 2 muons come from the same mother so the
+<BR>weight of a DD~ or BB~ is M*Br and they are no correlation in the decay
+(Br1*Br2)
+
+<P>If the muons are not correlated w12 = M1*M2*R1*R2*Br1*Br2 = w1*w2*R1*R2
+(in this method this gives wgt*fRate1*fRate2)
+<BR>&nbsp;</UL>
+*/
+//End_Html
+
 
-Float_t AliDimuCombinator::Weight(TParticle* part)
+Float_t AliDimuCombinator::Weight(TParticle* part) const
 {
 // Single muon weight
     return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
 }
 
-Bool_t  AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
+Bool_t  AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2) const
 {
 // Check if muons are correlated
 //
-    if (Origin(part1) == Origin(part2)) {
+    if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
+
        return kTRUE;
     } else {
        return kFALSE;
     }
 }
 
-TParticle* AliDimuCombinator::Parent(TParticle* part)
+TParticle* AliDimuCombinator::Parent(TParticle* part) const
 {
 // Return pointer to parent
 //
     return Particle(part->GetFirstMother());
 }
 
-Int_t AliDimuCombinator::Origin(TParticle* part)
+Int_t AliDimuCombinator::Origin(TParticle* part) const
 {
 // Return pointer to primary particle
 //
@@ -425,7 +477,7 @@ Int_t AliDimuCombinator::Origin(TParticle* part)
     return iparent;
 }
 
-Int_t AliDimuCombinator::Type(TParticle *part) 
+Int_t AliDimuCombinator::Type(TParticle *part)  const
 {
 // Return particle type for 
 return part->GetPdgCode();
@@ -434,14 +486,15 @@ return part->GetPdgCode();
 AliDimuCombinator& AliDimuCombinator::operator=(const  AliDimuCombinator& rhs)
 {
 // Assignment operator
+    rhs.Copy(*this);
     return *this;
 }
 
 
-void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
+void AliDimuCombinator::Copy(TObject&) const
 {
   //
-  // Copy *this onto lego -- not implemented
+  // Copy 
   //
   Fatal("Copy","Not implemented!\n");
 }