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