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 AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource):
 121     AliGenerator(FLUKAsource),
 122     fIkine(6),
 123     fAgeMax(1.e-5),
 124     fAddWeight(1.),
 125     fZshift(0.),
 126     fFrac(0.),
 127     fSourceId(-1),
 128     fFileName(0),
 129     fTreeChain(0),
 130     fTreeFluka(0),
 131     fIp(0.),
 132     fIpp(0.),
 133     fXi(0.),
 134     fYi(0.),
 135     fZi(0.),
 136     fPx(0.),
 137     fPy(0.),
 138     fPz(0.),
 139     fEkin(0.),
 140     fZv(0.),
 141     fRv(0.),
 142     fItra(0.),
 143     fIgas(0.),
 144     fWgt(0.),
 145     fEtag(0.),
 146     fPtg(0.),
 147     fAge(0.)
 148 {
 149 // Copy constructor
 150     FLUKAsource.Copy(*this);
 151 }
 152
 153
 154 //____________________________________________________________
 155 AliGenFLUKAsource::~AliGenFLUKAsource()
 156 {
 157 // Destructor
 158  if (fTreeFluka) delete fTreeFluka;
 159  if (fTreeChain) delete fTreeChain;
 160 }
 161
 162 void AliGenFLUKAsource::AddFile(const Text_t *filename)
 163 {
 164 // Add a file to the chain
 165     fTreeChain->Add(filename);
 166
 167 }
 168
 169
 170 //____________________________________________________________
 171 void AliGenFLUKAsource::FlukaInit()
 172 {
 173 // Set branch addresses of data entries
 174     TChain *h2=fTreeChain;
 175 //Set branch addresses
 176     h2->SetBranchAddress("Ip",&fIp);
 177     h2->SetBranchAddress("Ipp",&fIpp);
 178     h2->SetBranchAddress("Xi",&fXi);
 179     h2->SetBranchAddress("Yi",&fYi);
 180     h2->SetBranchAddress("Zi",&fZi);
 181     h2->SetBranchAddress("Px",&fPx);
 182     h2->SetBranchAddress("Py",&fPy);
 183     h2->SetBranchAddress("Pz",&fPz);
 184     h2->SetBranchAddress("Ekin",&fEkin);
 185     h2->SetBranchAddress("Zv",&fZv);
 186     h2->SetBranchAddress("Rv",&fRv);
 187     h2->SetBranchAddress("Itra",&fItra);
 188     h2->SetBranchAddress("Igas",&fIgas);
 189     h2->SetBranchAddress("Wgt",&fWgt);
 190     h2->SetBranchAddress("Etag",&fEtag);
 191     h2->SetBranchAddress("Ptg",&fPtg);
 192     h2->SetBranchAddress("Age",&fAge);
 193 }
 194
 195 //____________________________________________________________
 196 void AliGenFLUKAsource::Generate()
 197 {
 198 // Generate one event
 199
 200     const Int_t kIfluge[28]={kProton, kProtonBar, kElectron, kPositron,
 201                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 202                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 203                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 204                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 205                           0,kNuMu,kNuMuBar};
 206     Float_t polar[3]= {0,0,0};
 207   //
 208     Float_t origin[3];
 209     Float_t p[3];
 210     Float_t prwn;
 211     Float_t wgt, fwgt;
 212     Float_t phi;
 213     Float_t amass;
 214     Int_t iwgt;
 215     Int_t i, j, part, nt;
 216     static Int_t irwn=0;
 217     //
 218     Float_t random[2];
 219   //
 220     FlukaInit();
 221     TChain *h2=fTreeChain;
 222     Int_t nentries = (Int_t) h2->GetEntries();
 223     if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
 224
 225
 226   // loop over number of particles
 227     Int_t nb=0;
 228     Int_t ev=gMC->CurrentEvent();
 229     for (i=0; i<fNpart;i++) {
 230         Int_t entry=fNpart*(ev-1)+i;
 231         nb = (Int_t)h2->GetEvent(entry);
 232         if (irwn > nentries) {
 233             printf("No more entries in the FLUKA boundary source file\n");
 234             TFile *pFile=0;
 235             // Get AliRun object or create it
 236             if (!gAlice) {
 237                 gAlice = (AliRun*)pFile->Get("gAlice");
 238                 if (gAlice) printf("AliRun object found on file\n");
 239                 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 240             }
 241             TTree *fAli=gAlice->TreeK();
 242             if (fAli) pFile =fAli->GetCurrentFile();
 243             pFile->cd();
 244             printf("Generate - I'm out \n");
 245             return;
 246         }
 247
 248         Int_t ifip = Int_t(fIp);
 249
 250
 251         if (fSourceId != -1 && fIgas !=fSourceId) {
 252             irwn++;
 253             continue;
 254         }
 255
 256         if (ifip > 28 || ifip < 0) {
 257             irwn++;
 258             continue;
 259         }
 260
 261         if ((ifip != fIkine && fIkine != kAll && fIkine != kCharged
 262              && fIkine != 10) || fAge > fAgeMax){
 263             irwn++;
 264             continue;
 265         } else if (fIkine == kCharged) {
 266             if (ifip == 7 || ifip == 8 || fAge > fAgeMax) {
 267                 irwn++;
 268                 continue;
 269             }
 270         } else if (fIkine == kNoNeutron) {
 271             if (ifip == 8 || fAge > fAgeMax) {
 272                 irwn++;
 273                 continue;
 274             }
 275         }
 276
 277
 278         irwn++;
 279 //
 280 // PDG code from FLUKA particle type (ifip)
 281         part=kIfluge[int(ifip)-1];
 282 //
 283 // Calculate momentum from kinetic energy and mass of the particle
 284 #if ROOT_VERSION_CODE > 197895
 285         amass = gMC->ParticleMass(part);
 286 #else
 287         amass = (TDatabasePDG::Instance())->GetParticle(part)->Mass();
 288 #endif
 289         prwn=fEkin*sqrt(1. + 2.*amass/fEkin);
 290
 291
 292         origin[0]=fYi;
 293         origin[1]=fXi;
 294         origin[2]=fZi;
 295
 296         p[0]=fPy*prwn;
 297         p[1]=fPx*prwn;
 298         p[2]=fPz*prwn;
 299
 300         //handle particle weight correctly
 301         wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
 302         iwgt=Int_t(wgt);
 303         fwgt=wgt-Float_t(iwgt);
 304         Rndm(random,2);
 305         if (random[0] < fwgt) iwgt++;
 306         if (part==1 && iwgt>100) iwgt=100;
 307         Int_t nstack=0;
 308         for (j=0; j<iwgt; j++) {
 309             PushTrack(fTrackIt,-1,part,p,origin,polar,fAge,kPPrimary,nt);
 310             Rndm(random,2);
 311             phi=2*random[1]*TMath::Pi();
 312             Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
 313             Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
 314             p[0]=pn1;
 315             p[1]=pn2;
 316             Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
 317             Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
 318             origin[0]=on1;
 319             origin[1]=on2;
 320             nstack++;
 321         }
 322         if (nstack == 0) continue;
 323     }
 324
 325     TFile *pFile=0;
 326 // Get AliRun object or create it
 327     if (!gAlice) {
 328         gAlice = (AliRun*)pFile->Get("gAlice");
 329         if (gAlice) printf("AliRun object found on file\n");
 330         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 331     }
 332     TTree *fAli=gAlice->TreeK();
 333     if (fAli) pFile =fAli->GetCurrentFile();
 334     pFile->cd();
 335 }
 336
 337
 338 AliGenFLUKAsource& AliGenFLUKAsource::operator=(const  AliGenFLUKAsource& rhs)
 339 {
 340 // Assignment operator
 341     rhs.Copy(*this);
 342     return (*this);
 343 }
 344
 345
 346 void AliGenFLUKAsource::Copy(TObject &) const
 347 {
 348     Fatal("Copy","Not implemented!\n");
 349 }
 350
 351
 352
 353
 354