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 /* $Id$ */
  17
  18 // Read background particles from a FLUKA boundary source file
  19 // This is a very special generator that works for background studies for the muon-spectrometer.
  20 // The input files come from FLUKA simulations.
  21 // Files can be chained.
  22 // Author: andreas.morsch@cern.ch
  23
  24 #include <RVersion.h>
  25 #include "TPDGCode.h"
  26 #include "TDatabasePDG.h"
  27 #include <TVirtualMC.h>
  28
  29 #include "AliGenFLUKAsource.h"
  30 #include "AliRun.h"
  31
  32
  33 #include <TFile.h>
  34 #include <TTree.h>
  35 #include <TChain.h>
  36 #include <stdlib.h>
  37 ClassImp(AliGenFLUKAsource)
  38
  39 AliGenFLUKAsource::AliGenFLUKAsource()
  40          :AliGenerator(-1)
  41 {
  42     // Constructor
  43     fName="FLUKA";
  44     fTitle="FLUKA Boundary Source";
  45     // Read in all particle types by default
  46     fIkine=6;
  47     // Set maximum admitted age of particles to 1.e-05 by default
  48     fAgeMax=1.e-05;
  49     // Do not add weight
  50     fAddWeight=1.;
  51     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  52     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  53     // is no need to shift as it reads boundary source for the whole volume of
  54     // the Muon Arm; the default value corresponds to boundary source for the
  55     // whole volume of the MUON Arm
  56     fZshift=0;
  57     // Set the default file
  58     fFileName="";
  59
  60     fTreeFluka=0;
  61     fTreeChain = new TChain("h1");
  62 //
  63 //  Read all particles
  64     fNpart=-1;
  65 }
  66
  67 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
  68     :AliGenerator(npart)
  69 {
  70     // Constructor
  71     fName  = "FLUKA";
  72     fTitle = "FLUKA Boundary Source";
  73     // Read in all particle types by default
  74     fIkine=6;
  75     // Set maximum admitted age of particles to 1.e-05 by default
  76     fAgeMax=1.e-05;
  77     // Do not add weight
  78     fAddWeight=1.;
  79     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  80     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  81     // is no need to shift as it reads boundary source for the whole volume of
  82     // the Muon Arm; the default value corresponds to boundary source for the
  83     // whole volume of the MUON Arm
  84     fZshift=0;
  85     // Set the default file
  86     fFileName="";
  87
  88     fTreeFluka=0;
  89     fTreeChain = new TChain("h1");
  90     fSourceId=-1;
  91 }
  92
  93 AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource):
  94     AliGenerator(FLUKAsource)
  95 {
  96 // Copy constructor
  97     FLUKAsource.Copy(*this);
  98 }
  99
 100
 101 //____________________________________________________________
 102 AliGenFLUKAsource::~AliGenFLUKAsource()
 103 {
 104 // Destructor
 105  if (fTreeFluka) delete fTreeFluka;
 106  if (fTreeChain) delete fTreeChain;
 107 }
 108
 109 void AliGenFLUKAsource::AddFile(const Text_t *filename)
 110 {
 111 // Add a file to the chain
 112     fTreeChain->Add(filename);
 113
 114 }
 115
 116
 117 //____________________________________________________________
 118 void AliGenFLUKAsource::FlukaInit()
 119 {
 120 // Set branch addresses of data entries
 121     TChain *h2=fTreeChain;
 122 //Set branch addresses
 123     h2->SetBranchAddress("Ip",&fIp);
 124     h2->SetBranchAddress("Ipp",&fIpp);
 125     h2->SetBranchAddress("Xi",&fXi);
 126     h2->SetBranchAddress("Yi",&fYi);
 127     h2->SetBranchAddress("Zi",&fZi);
 128     h2->SetBranchAddress("Px",&fPx);
 129     h2->SetBranchAddress("Py",&fPy);
 130     h2->SetBranchAddress("Pz",&fPz);
 131     h2->SetBranchAddress("Ekin",&fEkin);
 132     h2->SetBranchAddress("Zv",&fZv);
 133     h2->SetBranchAddress("Rv",&fRv);
 134     h2->SetBranchAddress("Itra",&fItra);
 135     h2->SetBranchAddress("Igas",&fIgas);
 136     h2->SetBranchAddress("Wgt",&fWgt);
 137     h2->SetBranchAddress("Etag",&fEtag);
 138     h2->SetBranchAddress("Ptg",&fPtg);
 139     h2->SetBranchAddress("Age",&fAge);
 140 }
 141
 142 //____________________________________________________________
 143 void AliGenFLUKAsource::Generate()
 144 {
 145 // Generate one event
 146
 147     const Int_t kIfluge[28]={kProton, kProtonBar, kElectron, kPositron,
 148                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 149                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 150                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 151                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 152                           0,kNuMu,kNuMuBar};
 153     Float_t polar[3]= {0,0,0};
 154   //
 155     Float_t origin[3];
 156     Float_t p[3];
 157     Float_t prwn;
 158     Float_t wgt, fwgt;
 159     Float_t phi;
 160     Float_t amass;
 161     Int_t iwgt;
 162     Int_t i, j, part, nt;
 163     static Int_t irwn=0;
 164     //
 165     Float_t random[2];
 166   //
 167     FlukaInit();
 168     TChain *h2=fTreeChain;
 169     Int_t nentries = (Int_t) h2->GetEntries();
 170     if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
 171
 172
 173   // loop over number of particles
 174     Int_t nb=0;
 175     Int_t ev=gMC->CurrentEvent();
 176     for (i=0; i<fNpart;i++) {
 177         Int_t entry=fNpart*(ev-1)+i;
 178         nb = (Int_t)h2->GetEvent(entry);
 179         if (irwn > nentries) {
 180             printf("No more entries in the FLUKA boundary source file\n");
 181             TFile *pFile=0;
 182             // Get AliRun object or create it
 183             if (!gAlice) {
 184                 gAlice = (AliRun*)pFile->Get("gAlice");
 185                 if (gAlice) printf("AliRun object found on file\n");
 186                 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 187             }
 188             TTree *fAli=gAlice->TreeK();
 189             if (fAli) pFile =fAli->GetCurrentFile();
 190             pFile->cd();
 191             printf("Generate - I'm out \n");
 192             return;
 193         }
 194
 195         Int_t ifip = Int_t(fIp);
 196
 197
 198         if (fSourceId != -1 && fIgas !=fSourceId) {
 199             irwn++;
 200             continue;
 201         }
 202
 203         if (ifip > 28 || ifip < 0) {
 204             irwn++;
 205             continue;
 206         }
 207
 208         if ((ifip != fIkine && fIkine != kAll && fIkine != kCharged
 209              && fIkine != 10) || fAge > fAgeMax){
 210             irwn++;
 211             continue;
 212         } else if (fIkine == kCharged) {
 213             if (ifip == 7 || ifip == 8 || fAge > fAgeMax) {
 214                 irwn++;
 215                 continue;
 216             }
 217         } else if (fIkine == kNoNeutron) {
 218             if (ifip == 8 || fAge > fAgeMax) {
 219                 irwn++;
 220                 continue;
 221             }
 222         }
 223
 224
 225         irwn++;
 226 //
 227 // PDG code from FLUKA particle type (ifip)
 228         part=kIfluge[int(ifip)-1];
 229 //
 230 // Calculate momentum from kinetic energy and mass of the particle
 231 #if ROOT_VERSION_CODE > 197895
 232         amass = gMC->ParticleMass(part);
 233 #else
 234         amass = (TDatabasePDG::Instance())->GetParticle(part)->Mass();
 235 #endif
 236         prwn=fEkin*sqrt(1. + 2.*amass/fEkin);
 237
 238
 239         origin[0]=fYi;
 240         origin[1]=fXi;
 241         origin[2]=fZi;
 242
 243         p[0]=fPy*prwn;
 244         p[1]=fPx*prwn;
 245         p[2]=fPz*prwn;
 246
 247         //handle particle weight correctly
 248         wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
 249         iwgt=Int_t(wgt);
 250         fwgt=wgt-Float_t(iwgt);
 251         Rndm(random,2);
 252         if (random[0] < fwgt) iwgt++;
 253         if (part==1 && iwgt>100) iwgt=100;
 254         Int_t nstack=0;
 255         for (j=0; j<iwgt; j++) {
 256             PushTrack(fTrackIt,-1,part,p,origin,polar,fAge,kPPrimary,nt);
 257             Rndm(random,2);
 258             phi=2*random[1]*TMath::Pi();
 259             Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
 260             Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
 261             p[0]=pn1;
 262             p[1]=pn2;
 263             Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
 264             Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
 265             origin[0]=on1;
 266             origin[1]=on2;
 267             nstack++;
 268         }
 269         if (nstack == 0) continue;
 270     }
 271
 272     TFile *pFile=0;
 273 // Get AliRun object or create it
 274     if (!gAlice) {
 275         gAlice = (AliRun*)pFile->Get("gAlice");
 276         if (gAlice) printf("AliRun object found on file\n");
 277         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 278     }
 279     TTree *fAli=gAlice->TreeK();
 280     if (fAli) pFile =fAli->GetCurrentFile();
 281     pFile->cd();
 282 }
 283
 284
 285 AliGenFLUKAsource& AliGenFLUKAsource::operator=(const  AliGenFLUKAsource& rhs)
 286 {
 287 // Assignment operator
 288     rhs.Copy(*this);
 289     return (*this);
 290 }
 291
 292
 293 void AliGenFLUKAsource::Copy(TObject &) const
 294 {
 295     Fatal("Copy","Not implemented!\n");
 296 }
 297
 298
 299
 300
 301