EVGEN/AliGenFLUKAsource.cxx

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