1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.6 1999/11/09 07:38:48 fca
19 Changes for compatibility with version 2.23 of ROOT
21 Revision 1.5 1999/09/29 09:24:12 fca
22 Introduction of the Copyright and cvs Log
28 #include "AliGenExtFile.h"
29 #include "AliGenMUONlib.h"
32 #include <TDirectory.h>
36 ClassImp(AliGenExtFile)
37 AliGenExtFile::AliGenExtFile()
42 fTitle="Primaries from ext. File";
43 fFileName="dtujet93.root";
51 AliGenExtFile::AliGenExtFile(Int_t npart)
56 fTitle="Primaries from ext. File";
57 fFileName="dtujet93.root";
62 //____________________________________________________________
63 AliGenExtFile::~AliGenExtFile()
68 //____________________________________________________________
69 void AliGenExtFile::NtupleInit()
72 // reset the existing file environment and open a new root file if
73 // the pointer to the Fluka tree is null
78 File = new TFile(fFileName);
80 cout<<"I have opened "<<fFileName<<" file "<<endl;
82 // get the tree address in the Fluka boundary source file
83 fTreeNtuple = (TTree*)gDirectory->Get("h888");
85 File = fTreeNtuple->GetCurrentFile();
89 TTree *h2=fTreeNtuple;
90 //Set branch addresses
91 //Set branch addresses
92 h2->SetBranchAddress("Nihead",&Nihead);
93 h2->SetBranchAddress("Ihead",Ihead);
94 h2->SetBranchAddress("Nrhead",&Nrhead);
95 h2->SetBranchAddress("Rhead",Rhead);
96 h2->SetBranchAddress("Idpart",&Idpart);
97 h2->SetBranchAddress("Theta",&Theta);
98 h2->SetBranchAddress("Phi",&Phi);
99 h2->SetBranchAddress("P",&P);
100 h2->SetBranchAddress("E",&E);
104 //____________________________________________________________
105 void AliGenExtFile::Generate()
108 Float_t polar[3]= {0,0,0};
110 Float_t origin[3]={0,0,0};
114 Int_t i, j, nt, Ntracks=0;
117 TTree *h2=fTreeNtuple;
118 Int_t nentries = (Int_t) h2->GetEntries();
119 // loop over number of particles
120 Int_t nb = (Int_t)h2->GetEvent(fNcurrent);
124 for (j=0;j<3;j++) origin[j]=fOrigin[j];
125 if(fVertexSmear==perEvent) {
128 origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
129 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
133 if (fNcurrent >= nentries) {
134 printf("\n No more entries !!! !\n");
140 printf("\n This should never happen !\n");
142 printf("\n Next event contains %d tracks! \n", i6);
145 for (i=0; i<Ntracks; i++) {
146 Idpart=gMC->PDGFromId(Idpart);
147 Double_t amass = TDatabasePDG::Instance()->GetParticle(Idpart)->Mass();
149 Warning("Generate","Particle %d no %d E = %f mass = %f\n",
153 prwn=sqrt((E+amass)*(E-amass));
156 Theta *= TMath::Pi()/180.;
157 Phi = (Phi-180)*TMath::Pi()/180.;
158 if(Theta<fThetaMin || Theta>fThetaMax ||
159 Phi<fPhiMin || Phi>fPhiMax ||
160 prwn<fPMin || prwn>fPMax)
164 p[0]=prwn*TMath::Sin(Theta)*TMath::Cos(Phi);
165 p[1]=prwn*TMath::Sin(Theta)*TMath::Sin(Phi);
166 p[2]=prwn*TMath::Cos(Theta);
168 if(fVertexSmear==perTrack) {
172 +fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
173 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
176 gAlice->SetTrack(fTrackIt,-1,Idpart,p,origin,polar,0,"Primary",nt);
179 nb = (Int_t)h2->GetEvent(fNcurrent);
183 // Get AliRun object or create it
185 gAlice = (AliRun*)File->Get("gAlice");
186 if (gAlice) printf("AliRun object found on file\n");
187 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
189 TTree *fAli=gAlice->TreeK();
190 if (fAli) File =fAli->GetCurrentFile();