EVGEN/AliGenPairFlat.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 // Generator for particle pairs in a preset
  18 // kinematic range
  19 // ranges can be set for invariant mass, pair pT, pair rapidity
  20 // and pair azimuth
  21
  22 //
  23 // Comments and suggestions: markus.konrad.kohler@cern.ch
  24 //
  25
  26 #include "TPDGCode.h"
  27 #include <TDatabasePDG.h>
  28 #include "AliGenEventHeader.h"
  29 #include "TF1.h"
  30
  31 #include "AliGenPairFlat.h"
  32
  33 ClassImp(AliGenPairFlat)
  34
  35 //_____________________________________________________________________________
  36 AliGenPairFlat::AliGenPairFlat()
  37     :AliGenerator(),
  38      fPairNpart(10),
  39      fPairYMin(-10.),
  40      fPairYMax(10.),
  41      fPairPhiMin(0.),
  42      fPairPhiMax(TMath::TwoPi()),
  43      fPairPtMin(0.001),
  44      fPairPtMax(50.),
  45      fPairMassMin(0.),
  46      fPairMassMax(10.),
  47      fLegPdg1(11),
  48      fLegPdg2(-11),
  49      fAlpha(0.),
  50      fDebug(0)
  51 {
  52   //
  53   // Default constructor
  54   //
  55 }
  56
  57 //_____________________________________________________________________________
  58 void AliGenPairFlat::Generate()
  59 {
  60   //
  61   // Generate a random pair of particles
  62   //
  63
  64     Float_t polar[3]= {0,0,0};
  65     Float_t origin[3];
  66     Float_t time;
  67
  68     Float_t p[3];
  69     Int_t i, j, nt;
  70     Double_t phi, y, mass, mt, e, pt;
  71     Float_t weight = 1.;
  72     Double_t  pt1, pt2, y1, y2;
  73
  74     TLorentzVector mother, dau1, dau2;
  75     Float_t random[6];
  76
  77     TF1* pol = new TF1("Pol","1.+[0]*x*x",-1.,1.);
  78     pol->SetParameter(0, fAlpha);
  79
  80
  81     for (j=0;j<3;j++) origin[j]=fOrigin[j];
  82     time = fTimeOrigin;
  83     if(fVertexSmear==kPerEvent) {
  84         Vertex();
  85         for (j=0;j<3;j++) origin[j]=fVertex[j];
  86         time = fTime;
  87      }
  88
  89         printf("\n\n------------------GENERATOR SETTINGS------------------\n\n");
  90         printf("You choosed for the mother the Mass range %f - %f     \n",fPairMassMin,fPairMassMax);
  91         printf("You choosed for the mother the transverse Momentum range %f - %f \n",fPairPtMin,fPairPtMax);
  92         printf("You choosed for the mother the Phi range %f - %f \n",fPairPhiMin,fPairPhiMax);
  93         printf("The Particle will decay in (pdg) %i and %i\n \n",fLegPdg1, fLegPdg2);
  94         printf("The rest Mass of first daughter (%s) is %f \n",
  95         TDatabasePDG::Instance()->GetParticle(TMath::Abs(fLegPdg1))->GetTitle(),
  96         TDatabasePDG::Instance()->GetParticle(TMath::Abs(fLegPdg1))->Mass());
  97         printf("The rest Mass of second daughter (%s) is %f \n",
  98         TDatabasePDG::Instance()->GetParticle(TMath::Abs(fLegPdg2))->GetTitle(),
  99         TDatabasePDG::Instance()->GetParticle(TMath::Abs(fLegPdg2))->Mass());
 100         printf("polarization factor is alpha == %lf \n",fAlpha);
 101         printf("vertex is at x == %f || y == %f || z == %f   \n",origin[0],origin[1],origin[2]);
 102         printf("\n----------------------------------------------------------\n");
 103
 104
 105     for(i=0;i<fPairNpart;i++) {
 106
 107         // mother properties
 108         Rndm(random,4);
 109         mass    = fPairMassMin+random[0]*(fPairMassMax-fPairMassMin);
 110         pt      = fPairPtMin+random[1]*(fPairPtMax-fPairPtMin);
 111         y       = fPairYMin+random[2]*(fPairYMax-fPairYMin);
 112         phi     = fPairPhiMin+random[3]*(fPairPhiMax-fPairPhiMin);
 113
 114         mt      = TMath::Sqrt(pt*pt + mass*mass);
 115         p[0]    = pt*TMath::Cos(phi);
 116         p[1]    = pt*TMath::Sin(phi);
 117         p[2]    = mt*TMath::SinH(y);
 118         e       = mt*TMath::CosH(y);
 119
 120         mother.SetPxPyPzE(p[0],p[1],p[2],e);
 121
 122
 123          if (fDebug == 2)    printf("p = (%+11.4e,%+11.4e,%+11.4e) GeV || pt = %+11.4e || y =  %+11.4e || weight=%+11.4e || E= %+11.4e\n",p[0],p[1],p[2],pt, y,weight,e);
 124
 125         //decay procedure
 126         if(!Decay(mother,dau1,dau2,pol))continue;
 127
 128         pt1 = dau1.Pt();
 129         pt2 = dau2.Pt();
 130         y1 = dau1.Rapidity();
 131         y2 = dau2.Rapidity();
 132
 133         //first daughter
 134         p[0] = dau1.Px();
 135         p[1] = dau1.Py();
 136         p[2] = dau1.Pz();
 137         e = dau1.E();
 138
 139         if(fVertexSmear==kPerTrack) {
 140             Rndm(random,6);
 141             for (j=0;j<3;j++) {
 142                 origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 143                     TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 144             }
 145
 146             Rndm(random,2);
 147             time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
 148               TMath::Cos(2*random[0]*TMath::Pi())*
 149               TMath::Sqrt(-2*TMath::Log(random[1]));
 150         }
 151
 152         PushTrack(fTrackIt,-1,fLegPdg1,p,origin,polar,time,kPPrimary,nt, weight);
 153         if (fDebug == 2) printf("dau1-->id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV || pt = %+11.4e || weight=%+11.4e || E= %+11.4e\n",fLegPdg1,p[0],p[1],p[2],pt1,weight,e);
 154
 155         //second daughter
 156         p[0] = dau2.Px();
 157         p[1] = dau2.Py();
 158         p[2] = dau2.Pz();
 159         e = dau2.E();
 160
 161         if(fVertexSmear==kPerTrack) {
 162             Rndm(random,6);
 163             for (j=0;j<3;j++) {
 164                 origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 165                     TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 166             }
 167
 168             Rndm(random,2);
 169             time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
 170               TMath::Cos(2*random[0]*TMath::Pi())*
 171               TMath::Sqrt(-2*TMath::Log(random[1]));
 172         }
 173
 174         PushTrack(fTrackIt,-1,fLegPdg2,p,origin,polar,time,kPPrimary,nt,weight);
 175       if (fDebug == 2) printf("dau2-->id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV || pt = %+11.4e || weight=%+11.4e || E= %+11.4e\n",fLegPdg2,p[0],p[1],p[2],pt2, weight,e);
 176
 177     }//particle loop
 178
 179 }
 180
 181 //_____________________________________________________________________________
 182
 183 void AliGenPairFlat::Init()
 184 {
 185         //
 186         // Init
 187         //
 188
 189         printf("AliGenPairFlat::Init() not implemented!!!\n");
 190
 191 }
 192
 193
 194 //_____________________________________________________________________________
 195
 196
 197 Bool_t AliGenPairFlat::Decay(TLorentzVector& mother, TLorentzVector &dau1, TLorentzVector &dau2, TF1* polfactor)
 198 {
 199         //
 200         // decay procedure
 201         //
 202
 203         Double_t mp, md1, md2, ed1, ed2, pd1, pd2;
 204         Double_t costheta, sintheta, phi;
 205
 206         mp = mother.M();
 207         md1 = TDatabasePDG::Instance()->GetParticle(TMath::Abs(fLegPdg1))->Mass();
 208         md2 = TDatabasePDG::Instance()->GetParticle(TMath::Abs(fLegPdg2))->Mass();
 209
 210         if(mp < md1+md2){
 211         printf("Daughters are heavier than Mother!! Check Kinematics!! \n");
 212         return 0;
 213         }
 214
 215   ed1 = (mp*mp + md1*md1 - md2*md2)/(2.*mp);
 216   ed2 = mp-ed1;
 217   pd1 = TMath::Sqrt((ed1+md1)*(ed1-md1));
 218   pd2 = TMath::Sqrt((mp*mp-(md1+md2)*(md1+md2))*(mp*mp -(md1-md2)*(md1-md2)))/(2.*mp);
 219
 220   costheta = polfactor->GetRandom(); //polarization
 221   sintheta = TMath::Sqrt((1.+costheta)*(1.-costheta));
 222   phi      = 2.0*TMath::Pi()*gRandom->Rndm();
 223
 224         dau1.SetPxPyPzE(pd1*sintheta*TMath::Cos(phi),
 225                         pd1*sintheta*TMath::Sin(phi),
 226                         pd1*costheta,
 227                         ed1);
 228
 229         dau2.SetPxPyPzE((-1.)*pd1*sintheta*TMath::Cos(phi),
 230                         (-1.)*pd1*sintheta*TMath::Sin(phi),
 231                         (-1.)*pd1*costheta,
 232                         ed2);
 233
 234         TVector3 boost = mother.BoostVector();
 235
 236         dau1.Boost(boost);
 237         dau2.Boost(boost);
 238
 239   return 1;
 240
 241 }
 242
 243
 244