1 #include "AliGenFLUKAsource.h"
2 #include "AliGenMUONlib.h"
7 #include <TDirectory.h>
11 ClassImp(AliGenFLUKAsource)
12 AliGenFLUKAsource::AliGenFLUKAsource()
17 fTitle="FLUKA Boundary Source";
18 // Read in all particle types by default
20 // Set maximum admitted age of particles to 1.e-05 by default
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
30 // Set the default file
31 fFileName="flukasource.root";
40 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
45 fTitle="FLUKA Boundary Source";
46 // Read in all particle types by default
48 // Set maximum admitted age of particles to 1.e-05 by default
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
58 // Set the default file
59 fFileName="flukasource.root";
65 //____________________________________________________________
66 AliGenFLUKAsource::~AliGenFLUKAsource()
71 //____________________________________________________________
72 void AliGenFLUKAsource::FlukaInit()
75 // reset the existing file environment and open a new root file if
76 // the pointer to the Fluka tree is null
81 File = new TFile(fFileName);
83 cout<<"I have opened "<<fFileName<<" file "<<endl;
85 // get the tree address in the Fluka boundary source file
86 fTreeFluka = (TTree*)gDirectory->Get("h1");
88 File = fTreeFluka->GetCurrentFile();
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);
113 //____________________________________________________________
114 void AliGenFLUKAsource::Generate()
117 AliMC* pMC = AliMC::GetMC();
119 const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
120 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
121 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
122 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
123 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
125 Float_t polar[3]= {0,0,0};
133 Float_t amass, charge, tlife;
136 Int_t i, j, part, nt;
142 TTree *h2=fTreeFluka;
143 Int_t nentries = (Int_t) h2->GetEntries();
144 if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
146 // loop over number of particles
148 for (i=0; i<fNpart;i++) {
149 Int_t ev=pMC->CurrentEvent();
150 Int_t entry=fNpart*(ev-1)+i;
151 nb = (Int_t)h2->GetEvent(entry);
152 if (irwn > nentries) {
153 printf("No more entries in the FLUKA boundary source file\n");
155 // Get AliRun object or create it
157 gAlice = (AliRun*)File->Get("gAlice");
158 if (gAlice) printf("AliRun object found on file\n");
159 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
161 TTree *fAli=gAlice->TreeK();
162 if (fAli) File =fAli->GetCurrentFile();
164 printf("Generate - I'm out \n");
167 if (Ip > 28 || Ip < 0) {
172 if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
175 } else if (fIkine == 9) {
176 if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
184 printf("\n Particle type: %f \n \n ", Ip);
189 } else if (Ip == 8) {
190 prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
193 part=ifluge[int(Ip)-1];
194 pMC->Gfpart(part, name, itrtyp,
195 amass, charge, tlife);
196 prwn=sqrt(Ekin*Ekin + 2.*amass);
206 //handle particle weight correctly
207 wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
209 fwgt=wgt-Float_t(iwgt);
211 if (random[0] < fwgt) iwgt++;
212 if (part==1 && iwgt>100) iwgt=100;
214 for (j=0; j<iwgt; j++) {
215 gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
217 phi=2*random[1]*TMath::Pi();
218 Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
219 Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
222 Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
223 Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
228 if (nstack == 0) continue;
232 // Get AliRun object or create it
234 gAlice = (AliRun*)File->Get("gAlice");
235 if (gAlice) printf("AliRun object found on file\n");
236 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
238 TTree *fAli=gAlice->TreeK();
239 if (fAli) File =fAli->GetCurrentFile();