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