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