MFT/AliGenParamPionsKaons.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 //
  18 //      Parametric generator of primary pions and kaons
  19 //
  20 //      Contact author: antonio.uras@cern.ch
  21 //
  22 //====================================================================================================================================================
  23
  24 #include "TPDGCode.h"
  25
  26 #include "AliConst.h"
  27 #include "AliRun.h"
  28 #include "AliGenEventHeader.h"
  29 #include "TDatabasePDG.h"
  30 #include "AliPDG.h"
  31 #include "TFile.h"
  32 #include "TROOT.h"
  33 #include "AliGenParamPionsKaons.h"
  34
  35 ClassImp(AliGenParamPionsKaons)
  36
  37 //====================================================================================================================================================
  38
  39 AliGenParamPionsKaons::AliGenParamPionsKaons():
  40   AliGenerator(),
  41   fGeneratePion(kTRUE),
  42   fGenerateKaon(kTRUE),
  43   fPtVsRapidityPrimaryPosPions(0x0),
  44   fPtVsRapidityPrimaryNegPions(0x0),
  45   fPtVsRapidityPrimaryPosKaons(0x0),
  46   fPtVsRapidityPrimaryNegKaons(0x0),
  47   fHistPdgCode(0x0) {
  48
  49   // Default constructor
  50
  51 }
  52
  53 //====================================================================================================================================================
  54
  55 AliGenParamPionsKaons::AliGenParamPionsKaons(Int_t nPart, Char_t *inputFile):
  56   AliGenerator(nPart),
  57   fGeneratePion(kTRUE),
  58   fGenerateKaon(kTRUE),
  59   fPtVsRapidityPrimaryPosPions(0x0),
  60   fPtVsRapidityPrimaryNegPions(0x0),
  61   fPtVsRapidityPrimaryPosKaons(0x0),
  62   fPtVsRapidityPrimaryNegKaons(0x0),
  63   fHistPdgCode(0x0) {
  64
  65   // Standard constructor
  66
  67   fName  = "ParamPionsKaons";
  68   fTitle = "Parametric pions and kaons generator";
  69
  70   LoadInputHistos(inputFile);
  71
  72 }
  73
  74 //====================================================================================================================================================
  75
  76 void AliGenParamPionsKaons::Generate() {
  77
  78   // Generate one trigger
  79
  80   Double_t polar[3]= {0,0,0};
  81
  82   Double_t origin[3];
  83   Double_t p[3];
  84
  85   Double_t mass=0., pt=0., rap=0., mom=0., energy=0, mt=0., phi=0., time=0.;
  86   Int_t nt;
  87   Int_t pdgCode;
  88
  89   for (Int_t j=0; j<3; j++) origin[j] = fOrigin[j];
  90   time = fTimeOrigin;
  91   if (fVertexSmear==kPerEvent) {
  92     Vertex();
  93     for (Int_t j=0; j<3; j++) origin[j] = fVertex[j];
  94     time = fTime;
  95   }
  96
  97   Int_t nPartGenerated = 0;
  98
  99   while (nPartGenerated<fNpart) {
 100
 101     pdgCode = TMath::Nint(fHistPdgCode->GetRandom());
 102     if (TMath::Abs(pdgCode)==321 && !fGenerateKaon) continue;
 103     if (TMath::Abs(pdgCode)==211 && !fGeneratePion) continue;
 104
 105     mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
 106
 107     switch (pdgCode) {
 108     case 211:
 109       fPtVsRapidityPrimaryPosPions->GetRandom2(pt, rap);
 110       break;
 111     case -211:
 112       fPtVsRapidityPrimaryNegPions->GetRandom2(pt, rap);
 113       break;
 114     case 321:
 115       fPtVsRapidityPrimaryPosKaons->GetRandom2(pt, rap);
 116       break;
 117     case -321:
 118       fPtVsRapidityPrimaryNegKaons->GetRandom2(pt, rap);
 119       break;
 120     }
 121
 122     mt = TMath::Sqrt(mass*mass + pt*pt);
 123     energy = mt * TMath::CosH(rap);
 124     mom = TMath::Sqrt(energy*energy - mass*mass);
 125
 126     if (TestBit(kYRange))        if (rap<fYMin || rap>fYMax) continue;
 127     if (TestBit(kMomentumRange)) if (mom<fPMin || rap>fPMax) continue;
 128     if (TestBit(kPtRange))       if (pt<fPtMin || pt>fPtMax) continue;
 129
 130     phi = fPhiMin + gRandom->Rndm()*(fPhiMax-fPhiMin);
 131     p[0] = pt*TMath::Cos(phi);
 132     p[1] = pt*TMath::Sin(phi);
 133     p[2] = mt*TMath::SinH(rap);
 134
 135     PushTrack(1, -1, Int_t(pdgCode),
 136               p[0],p[1],p[2],energy,
 137               origin[0],origin[1],origin[2],Double_t(time),
 138               polar[0],polar[1],polar[2],
 139               kPPrimary, nt, 1., 1);
 140
 141     nPartGenerated++;
 142
 143   }
 144
 145   AliGenEventHeader* header = new AliGenEventHeader("ParamPionsKaons");
 146   header->SetPrimaryVertex(fVertex);
 147   header->SetNProduced(fNpart);
 148   header->SetInteractionTime(fTime);
 149
 150   // Passes header either to the container or to gAlice
 151   if (fContainer) {
 152     fContainer->AddHeader(header);
 153   }
 154   else {
 155     gAlice->SetGenEventHeader(header);
 156   }
 157
 158 }
 159
 160 //====================================================================================================================================================
 161
 162 void AliGenParamPionsKaons::Init() {
 163
 164   // Initialisation, check consistency of selected ranges
 165   if (TestBit(kPtRange) && TestBit(kMomentumRange))
 166     Fatal("Init","You should not set the momentum range and the pt range at the same time!\n");
 167   if ((!TestBit(kPtRange)) && (!TestBit(kMomentumRange)))
 168     Fatal("Init","You should set either the momentum or the pt range!\n");
 169   if ((TestBit(kYRange) && TestBit(kThetaRange)) || (TestBit(kYRange) && TestBit(kEtaRange)) || (TestBit(kEtaRange) && TestBit(kThetaRange)))
 170     Fatal("Init","You should only set the range of one of these variables: y, eta or theta\n");
 171   if ((!TestBit(kYRange)) && (!TestBit(kEtaRange)) && (!TestBit(kThetaRange)))
 172     Fatal("Init","You should set the range of one of these variables: y, eta or theta\n");
 173
 174   AliPDG::AddParticlesToPdgDataBase();
 175
 176 }
 177
 178 //====================================================================================================================================================
 179
 180 void AliGenParamPionsKaons::LoadInputHistos(Char_t *inputFile) {
 181
 182   TFile *fileIn = new TFile(inputFile);
 183
 184   TH2D *myPtVsRapidityPrimaryPosPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosPions");
 185   TH2D *myPtVsRapidityPrimaryNegPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegPions");
 186   TH2D *myPtVsRapidityPrimaryPosKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosKaons");
 187   TH2D *myPtVsRapidityPrimaryNegKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegKaons");
 188   TH1D *myHistPdgCode                 = (TH1D*) fileIn->Get("fHistPdgCode");
 189
 190   myPtVsRapidityPrimaryPosPions -> SetName("myPtVsRapidityPrimaryPosPions");
 191   myPtVsRapidityPrimaryNegPions -> SetName("myPtVsRapidityPrimaryNegPions");
 192   myPtVsRapidityPrimaryPosKaons -> SetName("myPtVsRapidityPrimaryPosKaons");
 193   myPtVsRapidityPrimaryNegKaons -> SetName("myPtVsRapidityPrimaryNegKaons");
 194   myHistPdgCode                 -> SetName("myHistPdgCode");
 195
 196   fPtVsRapidityPrimaryPosPions = new TH2D("fPtVsRapidityPrimaryPosPions","",
 197                                           myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(),
 198                                           myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmin(),
 199                                           myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmax(),
 200                                           myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(),
 201                                           myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmin(),
 202                                           myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmax());
 203   for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(); iBinX++) {
 204     for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(); iBinY++) {
 205       fPtVsRapidityPrimaryPosPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosPions->GetBinContent(iBinX+1,iBinY+1));
 206     }
 207   }
 208
 209   fPtVsRapidityPrimaryNegPions = new TH2D("fPtVsRapidityPrimaryNegPions","",
 210                                           myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(),
 211                                           myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmin(),
 212                                           myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmax(),
 213                                           myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(),
 214                                           myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmin(),
 215                                           myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmax());
 216   for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(); iBinX++) {
 217     for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(); iBinY++) {
 218       fPtVsRapidityPrimaryNegPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegPions->GetBinContent(iBinX+1,iBinY+1));
 219     }
 220   }
 221
 222   fPtVsRapidityPrimaryPosKaons = new TH2D("fPtVsRapidityPrimaryPosKaons","",
 223                                           myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(),
 224                                           myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmin(),
 225                                           myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmax(),
 226                                           myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(),
 227                                           myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmin(),
 228                                           myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmax());
 229   for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(); iBinX++) {
 230     for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(); iBinY++) {
 231       fPtVsRapidityPrimaryPosKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosKaons->GetBinContent(iBinX+1,iBinY+1));
 232     }
 233   }
 234
 235   fPtVsRapidityPrimaryNegKaons = new TH2D("fPtVsRapidityPrimaryNegKaons","",
 236                                           myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(),
 237                                           myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmin(),
 238                                           myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmax(),
 239                                           myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(),
 240                                           myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmin(),
 241                                           myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmax());
 242   for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(); iBinX++) {
 243     for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(); iBinY++) {
 244       fPtVsRapidityPrimaryNegKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegKaons->GetBinContent(iBinX+1,iBinY+1));
 245     }
 246   }
 247
 248   fHistPdgCode = new TH1D("fHistPdgCode","",
 249                           myHistPdgCode->GetXaxis()->GetNbins(),
 250                           myHistPdgCode->GetXaxis()->GetXmin(),
 251                           myHistPdgCode->GetXaxis()->GetXmax());
 252   for (Int_t iBinX=0; iBinX<myHistPdgCode->GetXaxis()->GetNbins(); iBinX++) {
 253     fHistPdgCode->SetBinContent(iBinX+1, myHistPdgCode->GetBinContent(iBinX+1));
 254   }
 255
 256 //   fHistPdgCode = new TH1D("fHistPdgCode", "fHistPdgCode", 10, 0., 10);
 257 //   fHistPdgCode -> Fill(3.);
 258
 259 }
 260
 261 //====================================================================================================================================================