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