/*
$Log$
+Revision 1.7 1999/11/03 17:43:20 fca
+New version from G.Martinez & A.Morsch
+
+Revision 1.6 1999/09/29 09:24:12 fca
+Introduction of the Copyright and cvs Log
+
*/
#include "AliGenFLUKAsource.h"
#include "AliMC.h"
#include "AliRun.h"
#include "AliPDG.h"
-
#include <TDirectory.h>
#include <TFile.h>
#include <TTree.h>
fFileName="flukasource.root";
fTreeFluka=0;
+ fTreeChain = new TChain("h1");
//
// Read all particles
fNpart=-1;
fFileName="flukasource.root";
fTreeFluka=0;
-
+ fTreeChain = new TChain("h1");
}
//____________________________________________________________
AliGenFLUKAsource::~AliGenFLUKAsource()
{
- delete fTreeFluka;
+ if (fTreeFluka) delete fTreeFluka;
+ if (fTreeChain) delete fTreeChain;
+ // if (fFileName) delete fFileName;
}
-//____________________________________________________________
-void AliGenFLUKAsource::FlukaInit()
+void AliGenFLUKAsource::AddFile(const Text_t *filename)
{
-//
-// reset the existing file environment and open a new root file if
-// the pointer to the Fluka tree is null
+ fTreeChain->Add(filename);
- TFile *File=0;
- if (fTreeFluka==0) {
- if (!File) {
- File = new TFile(fFileName);
- File->cd();
- cout<<"I have opened "<<fFileName<<" file "<<endl;
- }
-// get the tree address in the Fluka boundary source file
- fTreeFluka = (TTree*)gDirectory->Get("h1");
- } else {
- File = fTreeFluka->GetCurrentFile();
- File->cd();
- }
+}
- TTree *h2=fTreeFluka;
+
+//____________________________________________________________
+void AliGenFLUKAsource::FlukaInit()
+{
+ TChain *h2=fTreeChain;
//Set branch addresses
h2->SetBranchAddress("Ip",&Ip);
h2->SetBranchAddress("Ipp",&Ipp);
//____________________________________________________________
void AliGenFLUKAsource::Generate()
-{
+{
+ AliMC* gMC = AliMC::GetMC();
- const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
+ const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
0,kNuMu,kNuMuBar};
- Float_t polar[3]= {0,0,0};
+ Float_t polar[3]= {0,0,0};
//
- Float_t origin[3];
- Float_t p[3];
- Float_t prwn;
- Float_t wgt, fwgt;
- Float_t phi;
- char name[100];
- Float_t amass, charge, tlife;
- Int_t itrtyp;
- Int_t iwgt;
- Int_t i, j, part, nt;
- static Int_t irwn;
- //
- Float_t random[2];
+ Float_t origin[3];
+ Float_t p[3];
+ Float_t prwn;
+ Float_t wgt, fwgt;
+ Float_t phi;
+ char name[100];
+ Float_t amass, charge, tlife;
+ Int_t itrtyp;
+ Int_t iwgt;
+ Int_t i, j, part, nt;
+ static Int_t irwn;
+ //
+ Float_t random[2];
//
- FlukaInit();
- TTree *h2=fTreeFluka;
- Int_t nentries = (Int_t) h2->GetEntries();
- if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
+ FlukaInit();
+ TChain *h2=fTreeChain;
+ Int_t nentries = (Int_t) h2->GetEntries();
+ if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
// loop over number of particles
- Int_t nb=0;
- for (i=0; i<fNpart;i++) {
+ Int_t nb=0;
Int_t ev=gMC->CurrentEvent();
- Int_t entry=fNpart*(ev-1)+i;
- nb = (Int_t)h2->GetEvent(entry);
- if (irwn > nentries) {
- printf("No more entries in the FLUKA boundary source file\n");
- TFile *File=0;
- // Get AliRun object or create it
- if (!gAlice) {
- gAlice = (AliRun*)File->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
- TTree *fAli=gAlice->TreeK();
- if (fAli) File =fAli->GetCurrentFile();
- File->cd();
- printf("Generate - I'm out \n");
- return;
- }
- if (Ip > 28 || Ip < 0) {
- irwn++;
- continue;
- }
-
- if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
- irwn++;
- continue;
- } else if (fIkine == 9) {
- if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
- irwn++;
- continue;
- }
- }
+ for (i=0; i<fNpart;i++) {
+ Int_t entry=fNpart*(ev-1)+i;
+ nb = (Int_t)h2->GetEvent(entry);
+ if (irwn > nentries) {
+ printf("No more entries in the FLUKA boundary source file\n");
+ TFile *File=0;
+ // Get AliRun object or create it
+ if (!gAlice) {
+ gAlice = (AliRun*)File->Get("gAlice");
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+ }
+ TTree *fAli=gAlice->TreeK();
+ if (fAli) File =fAli->GetCurrentFile();
+ File->cd();
+ printf("Generate - I'm out \n");
+ return;
+ }
+ if (Ip > 28 || Ip < 0) {
+ irwn++;
+ continue;
+ }
+
+ if ((Ip != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) || Age > fAgeMax){
+ irwn++;
+ continue;
+ } else if (fIkine == 9) {
+ if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
+ irwn++;
+ continue;
+ }
+ } else if (fIkine == 10) {
+ if (Ip == 8 || Age > fAgeMax) {
+ irwn++;
+ continue;
+ }
+ }
- irwn++;
- printf("\n Particle type: %f \n \n ", Ip);
-
- if (Ip ==7){
- prwn=Ekin;
- part=1;
- } else if (Ip == 8) {
- prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
- part=13;
- } else {
- part=ifluge[int(Ip)-1];
- gMC->Gfpart(part, name, itrtyp,
- amass, charge, tlife);
- prwn=sqrt(Ekin*Ekin + 2.*amass);
- }
- origin[0]=Xi;
- origin[1]=Yi;
- origin[2]=Zi;
+ irwn++;
+//
+// PDG code from FLUKA particle type (Ip)
+ part=ifluge[int(Ip)-1];
+//
+// Calculate momentum from kinetic energy and mass of the particle
+ gMC->Gfpart(part, name, itrtyp,
+ amass, charge, tlife);
+ prwn=Ekin*sqrt(1. + 2.*amass/Ekin);
+
- p[0]=Px*prwn;
- p[1]=Py*prwn;
- p[2]=Pz*prwn;
+ origin[0]=Yi;
+ origin[1]=Xi;
+ origin[2]=Zi;
+
+ p[0]=Py*prwn;
+ p[1]=Px*prwn;
+ p[2]=Pz*prwn;
- //handle particle weight correctly
- wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
- iwgt=Int_t(wgt);
- fwgt=wgt-Float_t(iwgt);
- gMC->Rndm(random,2);
- if (random[0] < fwgt) iwgt++;
- if (part==1 && iwgt>100) iwgt=100;
- Int_t nstack=0;
- for (j=0; j<iwgt; j++) {
- gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
+ //handle particle weight correctly
+ wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
+ iwgt=Int_t(wgt);
+ fwgt=wgt-Float_t(iwgt);
gMC->Rndm(random,2);
- phi=2*random[1]*TMath::Pi();
- Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
- Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
- p[0]=pn1;
- p[1]=pn2;
- Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
- Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
- origin[0]=on1;
- origin[1]=on2;
- nstack++;
+ if (random[0] < fwgt) iwgt++;
+ if (part==1 && iwgt>100) iwgt=100;
+ Int_t nstack=0;
+ for (j=0; j<iwgt; j++) {
+ gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,Age,"Primary",nt);
+ gMC->Rndm(random,2);
+ phi=2*random[1]*TMath::Pi();
+ Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
+ Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
+ p[0]=pn1;
+ p[1]=pn2;
+ Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
+ Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
+ origin[0]=on1;
+ origin[1]=on2;
+ nstack++;
+ }
+ if (nstack == 0) continue;
}
- if (nstack == 0) continue;
- }
-
+
TFile *File=0;
// Get AliRun object or create it
if (!gAlice) {