EVGEN/AliGenFLUKAsource.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.8  2000/02/14 14:49:38  morsch
  19 Correct particle type for gamma and neutrons
  20 More consistent calculation of momentum from kin. energy and mass
  21
  22 Revision 1.7  1999/11/03 17:43:20  fca
  23 New version from G.Martinez & A.Morsch
  24
  25 Revision 1.6  1999/09/29 09:24:12  fca
  26 Introduction of the Copyright and cvs Log
  27
  28 */
  29
  30 #include "AliGenFLUKAsource.h"
  31 #include "AliGenMUONlib.h"
  32 #include "AliMC.h"
  33 #include "AliRun.h"
  34 #include "AliPDG.h"
  35 #include <TDirectory.h>
  36 #include <TFile.h>
  37 #include <TTree.h>
  38 #include <stdlib.h>
  39  ClassImp(AliGenFLUKAsource)
  40      AliGenFLUKAsource::AliGenFLUKAsource()
  41          :AliGenerator(-1)
  42 {
  43     //
  44     fName="FLUKA";
  45     fTitle="FLUKA Boundary Source";
  46     // Read in all particle types by default
  47     fIkine=6;
  48     // Set maximum admitted age of particles to 1.e-05 by default
  49     fAgeMax=1.e-05;
  50     // Do not add weight
  51     fAddWeight=1.;
  52     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  53     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  54     // is no need to shift as it reads boundary source for the whole volume of
  55     // the Muon Arm; the default value corresponds to boundary source for the
  56     // whole volume of the MUON Arm
  57     fZshift=0;
  58     // Set the default file
  59     fFileName="flukasource.root";
  60
  61     fTreeFluka=0;
  62     fTreeChain = new TChain("h1");
  63 //
  64 //  Read all particles
  65     fNpart=-1;
  66
  67 }
  68
  69 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
  70     :AliGenerator(npart)
  71 {
  72     //
  73     fName="FLUKA";
  74     fTitle="FLUKA Boundary Source";
  75     // Read in all particle types by default
  76     fIkine=6;
  77     // Set maximum admitted age of particles to 1.e-05 by default
  78     fAgeMax=1.e-05;
  79     // Do not add weight
  80     fAddWeight=1.;
  81     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  82     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  83     // is no need to shift as it reads boundary source for the whole volume of
  84     // the Muon Arm; the default value corresponds to boundary source for the
  85     // whole volume of the MUON Arm
  86     fZshift=0;
  87     // Set the default file
  88     fFileName="flukasource.root";
  89
  90     fTreeFluka=0;
  91     fTreeChain = new TChain("h1");
  92 }
  93
  94 //____________________________________________________________
  95 AliGenFLUKAsource::~AliGenFLUKAsource()
  96 {
  97  if (fTreeFluka) delete fTreeFluka;
  98  if (fTreeChain) delete fTreeChain;
  99  // if (fFileName)  delete fFileName;
 100 }
 101
 102 void AliGenFLUKAsource::AddFile(const Text_t *filename)
 103 {
 104     fTreeChain->Add(filename);
 105
 106 }
 107
 108
 109 //____________________________________________________________
 110 void AliGenFLUKAsource::FlukaInit()
 111 {
 112     TChain *h2=fTreeChain;
 113 //Set branch addresses
 114     h2->SetBranchAddress("Ip",&Ip);
 115     h2->SetBranchAddress("Ipp",&Ipp);
 116     h2->SetBranchAddress("Xi",&Xi);
 117     h2->SetBranchAddress("Yi",&Yi);
 118     h2->SetBranchAddress("Zi",&Zi);
 119     h2->SetBranchAddress("Px",&Px);
 120     h2->SetBranchAddress("Py",&Py);
 121     h2->SetBranchAddress("Pz",&Pz);
 122     h2->SetBranchAddress("Ekin",&Ekin);
 123     h2->SetBranchAddress("Zv",&Zv);
 124     h2->SetBranchAddress("Rv",&Rv);
 125     h2->SetBranchAddress("Itra",&Itra);
 126     h2->SetBranchAddress("Igas",&Igas);
 127     h2->SetBranchAddress("Wgt",&Wgt);
 128     h2->SetBranchAddress("Etag",&Etag);
 129     h2->SetBranchAddress("Ptg",&Ptg);
 130     h2->SetBranchAddress("Age",&Age);
 131 }
 132
 133 //____________________________________________________________
 134 void AliGenFLUKAsource::Generate()
 135 {
 136     AliMC* gMC = AliMC::GetMC();
 137
 138     const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
 139                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 140                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 141                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 142                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 143                           0,kNuMu,kNuMuBar};
 144     Float_t polar[3]= {0,0,0};
 145   //
 146     Float_t origin[3];
 147     Float_t p[3];
 148     Float_t prwn;
 149     Float_t wgt, fwgt;
 150     Float_t phi;
 151     char name[100];
 152     Float_t amass, charge, tlife;
 153     Int_t itrtyp;
 154     Int_t iwgt;
 155     Int_t i, j, part, nt;
 156     static Int_t irwn=0;
 157     //
 158     Float_t random[2];
 159   //
 160     FlukaInit();
 161     TChain *h2=fTreeChain;
 162     Int_t nentries = (Int_t) h2->GetEntries();
 163     if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
 164
 165   // loop over number of particles
 166     Int_t nb=0;
 167     Int_t ev=gMC->CurrentEvent();
 168     for (i=0; i<fNpart;i++) {
 169         Int_t entry=fNpart*(ev-1)+i;
 170         nb = (Int_t)h2->GetEvent(entry);
 171         if (irwn > nentries) {
 172             printf("No more entries in the FLUKA boundary source file\n");
 173             TFile *File=0;
 174             // Get AliRun object or create it
 175             if (!gAlice) {
 176                 gAlice = (AliRun*)File->Get("gAlice");
 177                 if (gAlice) printf("AliRun object found on file\n");
 178                 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 179             }
 180             TTree *fAli=gAlice->TreeK();
 181             if (fAli) File =fAli->GetCurrentFile();
 182             File->cd();
 183             printf("Generate - I'm out \n");
 184             return;
 185         }
 186         if (Ip > 28 || Ip < 0) {
 187             irwn++;
 188             continue;
 189         }
 190
 191         if ((Ip != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) || Age > fAgeMax){
 192             irwn++;
 193             continue;
 194         } else if (fIkine == 9) {
 195             if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
 196                 irwn++;
 197                 continue;
 198             }
 199         } else if (fIkine == 10) {
 200             if (Ip == 8 || Age > fAgeMax) {
 201                 irwn++;
 202                 continue;
 203             }
 204         }
 205
 206
 207         irwn++;
 208 //
 209 // PDG code from FLUKA particle type (Ip)
 210         part=ifluge[int(Ip)-1];
 211 //
 212 // Calculate momentum from kinetic energy and mass of the particle
 213         gMC->Gfpart(part, name, itrtyp,
 214                     amass, charge, tlife);
 215         prwn=Ekin*sqrt(1. + 2.*amass/Ekin);
 216
 217
 218         origin[0]=Yi;
 219         origin[1]=Xi;
 220         origin[2]=Zi;
 221
 222         p[0]=Py*prwn;
 223         p[1]=Px*prwn;
 224         p[2]=Pz*prwn;
 225
 226         //handle particle weight correctly
 227         wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
 228         iwgt=Int_t(wgt);
 229         fwgt=wgt-Float_t(iwgt);
 230         gMC->Rndm(random,2);
 231         if (random[0] < fwgt) iwgt++;
 232         if (part==1 && iwgt>100) iwgt=100;
 233         Int_t nstack=0;
 234         for (j=0; j<iwgt; j++) {
 235             gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,Age,"Primary",nt);
 236             gMC->Rndm(random,2);
 237             phi=2*random[1]*TMath::Pi();
 238             Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
 239             Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
 240             p[0]=pn1;
 241             p[1]=pn2;
 242             Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
 243             Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
 244             origin[0]=on1;
 245             origin[1]=on2;
 246             nstack++;
 247         }
 248         if (nstack == 0) continue;
 249     }
 250
 251     TFile *File=0;
 252 // Get AliRun object or create it
 253     if (!gAlice) {
 254         gAlice = (AliRun*)File->Get("gAlice");
 255         if (gAlice) printf("AliRun object found on file\n");
 256         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 257     }
 258     TTree *fAli=gAlice->TreeK();
 259     if (fAli) File =fAli->GetCurrentFile();
 260     File->cd();
 261 }
 262
 263
 264
 265
 266
 267
 268
 269
 270