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