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