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 */
  19
  20 #include "AliGenFLUKAsource.h"
  21 #include "AliGenMUONlib.h"
  22 #include "AliMC.h"
  23 #include "AliRun.h"
  24 #include "AliPDG.h"
  25
  26 #include <TDirectory.h>
  27 #include <TFile.h>
  28 #include <TTree.h>
  29 #include <stdlib.h>
  30  ClassImp(AliGenFLUKAsource)
  31      AliGenFLUKAsource::AliGenFLUKAsource()
  32          :AliGenerator(-1)
  33 {
  34     //
  35     fName="FLUKA";
  36     fTitle="FLUKA Boundary Source";
  37     // Read in all particle types by default
  38     fIkine=6;
  39     // Set maximum admitted age of particles to 1.e-05 by default
  40     fAgeMax=1.e-05;
  41     // Do not add weight
  42     fAddWeight=1.;
  43     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  44     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  45     // is no need to shift as it reads boundary source for the whole volume of
  46     // the Muon Arm; the default value corresponds to boundary source for the
  47     // whole volume of the MUON Arm
  48     fZshift=0;
  49     // Set the default file
  50     fFileName="flukasource.root";
  51
  52     fTreeFluka=0;
  53 //
  54 //  Read all particles
  55     fNpart=-1;
  56
  57 }
  58
  59 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
  60     :AliGenerator(npart)
  61 {
  62     //
  63     fName="FLUKA";
  64     fTitle="FLUKA Boundary Source";
  65     // Read in all particle types by default
  66     fIkine=6;
  67     // Set maximum admitted age of particles to 1.e-05 by default
  68     fAgeMax=1.e-05;
  69     // Do not add weight
  70     fAddWeight=1.;
  71     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  72     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  73     // is no need to shift as it reads boundary source for the whole volume of
  74     // the Muon Arm; the default value corresponds to boundary source for the
  75     // whole volume of the MUON Arm
  76     fZshift=0;
  77     // Set the default file
  78     fFileName="flukasource.root";
  79
  80     fTreeFluka=0;
  81
  82 }
  83
  84 //____________________________________________________________
  85 AliGenFLUKAsource::~AliGenFLUKAsource()
  86 {
  87     delete fTreeFluka;
  88 }
  89
  90 //____________________________________________________________
  91 void AliGenFLUKAsource::FlukaInit()
  92 {
  93 //
  94 // reset the existing file environment and open a new root file if
  95 // the pointer to the Fluka tree is null
  96
  97     TFile *File=0;
  98     if (fTreeFluka==0) {
  99         if (!File) {
 100             File = new TFile(fFileName);
 101             File->cd();
 102             cout<<"I have opened "<<fFileName<<" file "<<endl;
 103         }
 104 // get the tree address in the Fluka boundary source file
 105         fTreeFluka = (TTree*)gDirectory->Get("h1");
 106     } else {
 107         File = fTreeFluka->GetCurrentFile();
 108         File->cd();
 109     }
 110
 111     TTree *h2=fTreeFluka;
 112 //Set branch addresses
 113     h2->SetBranchAddress("Ip",&Ip);
 114     h2->SetBranchAddress("Ipp",&Ipp);
 115     h2->SetBranchAddress("Xi",&Xi);
 116     h2->SetBranchAddress("Yi",&Yi);
 117     h2->SetBranchAddress("Zi",&Zi);
 118     h2->SetBranchAddress("Px",&Px);
 119     h2->SetBranchAddress("Py",&Py);
 120     h2->SetBranchAddress("Pz",&Pz);
 121     h2->SetBranchAddress("Ekin",&Ekin);
 122     h2->SetBranchAddress("Zv",&Zv);
 123     h2->SetBranchAddress("Rv",&Rv);
 124     h2->SetBranchAddress("Itra",&Itra);
 125     h2->SetBranchAddress("Igas",&Igas);
 126     h2->SetBranchAddress("Wgt",&Wgt);
 127     h2->SetBranchAddress("Etag",&Etag);
 128     h2->SetBranchAddress("Ptg",&Ptg);
 129     h2->SetBranchAddress("Age",&Age);
 130 }
 131
 132 //____________________________________________________________
 133 void AliGenFLUKAsource::Generate()
 134 {
 135
 136   const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
 137                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 138                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 139                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 140                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 141                           0,kNuMu,kNuMuBar};
 142   Float_t polar[3]= {0,0,0};
 143   //
 144   Float_t origin[3];
 145   Float_t p[3];
 146   Float_t prwn;
 147   Float_t wgt, fwgt;
 148   Float_t phi;
 149   char name[100];
 150   Float_t amass, charge, tlife;
 151   Int_t itrtyp;
 152   Int_t iwgt;
 153   Int_t i, j, part, nt;
 154   static Int_t irwn;
 155   //
 156   Float_t random[2];
 157   //
 158   FlukaInit();
 159   TTree *h2=fTreeFluka;
 160   Int_t nentries = (Int_t) h2->GetEntries();
 161   if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
 162
 163   // loop over number of particles
 164   Int_t nb=0;
 165   for (i=0; i<fNpart;i++) {
 166     Int_t ev=gMC->CurrentEvent();
 167     Int_t entry=fNpart*(ev-1)+i;
 168     nb = (Int_t)h2->GetEvent(entry);
 169     if (irwn > nentries) {
 170       printf("No more entries in the FLUKA boundary source file\n");
 171       TFile *File=0;
 172       // Get AliRun object or create it
 173       if (!gAlice) {
 174         gAlice = (AliRun*)File->Get("gAlice");
 175         if (gAlice) printf("AliRun object found on file\n");
 176         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 177       }
 178       TTree *fAli=gAlice->TreeK();
 179       if (fAli) File =fAli->GetCurrentFile();
 180       File->cd();
 181       printf("Generate - I'm out \n");
 182       return;
 183     }
 184     if (Ip > 28 || Ip < 0) {
 185       irwn++;
 186       continue;
 187     }
 188
 189     if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
 190        irwn++;
 191        continue;
 192     } else if (fIkine == 9) {
 193        if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
 194           irwn++;
 195           continue;
 196        }
 197     }
 198
 199
 200     irwn++;
 201     printf("\n Particle type: %f \n \n ", Ip);
 202
 203     if (Ip ==7){
 204        prwn=Ekin;
 205        part=1;
 206     } else if (Ip == 8) {
 207        prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
 208        part=13;
 209     } else {
 210        part=ifluge[int(Ip)-1];
 211        gMC->Gfpart(part, name, itrtyp,
 212                    amass, charge, tlife);
 213        prwn=sqrt(Ekin*Ekin + 2.*amass);
 214     }
 215     origin[0]=Xi;
 216     origin[1]=Yi;
 217     origin[2]=Zi;
 218
 219     p[0]=Px*prwn;
 220     p[1]=Py*prwn;
 221     p[2]=Pz*prwn;
 222
 223     //handle particle weight correctly
 224     wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
 225     iwgt=Int_t(wgt);
 226     fwgt=wgt-Float_t(iwgt);
 227     gMC->Rndm(random,2);
 228     if (random[0] < fwgt) iwgt++;
 229     if (part==1 && iwgt>100) iwgt=100;
 230     Int_t nstack=0;
 231     for (j=0; j<iwgt; j++) {
 232         gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
 233         gMC->Rndm(random,2);
 234         phi=2*random[1]*TMath::Pi();
 235         Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
 236         Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
 237         p[0]=pn1;
 238         p[1]=pn2;
 239         Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
 240         Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
 241         origin[0]=on1;
 242         origin[1]=on2;
 243         nstack++;
 244     }
 245     if (nstack == 0) continue;
 246   }
 247
 248     TFile *File=0;
 249 // Get AliRun object or create it
 250     if (!gAlice) {
 251         gAlice = (AliRun*)File->Get("gAlice");
 252         if (gAlice) printf("AliRun object found on file\n");
 253         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 254     }
 255     TTree *fAli=gAlice->TreeK();
 256     if (fAli) File =fAli->GetCurrentFile();
 257     File->cd();
 258 }
 259
 260
 261
 262
 263
 264
 265
 266
 267