EVGEN/AliGenExtFile.cxx

   1 #include "AliGenExtFile.h"
   2 #include "AliGenMUONlib.h"
   3 #include "AliMC.h"
   4 #include "AliRun.h"
   5 #include <TDirectory.h>
   6 #include <TFile.h>
   7 #include <TTree.h>
   8 #include <stdlib.h>
   9  ClassImp(AliGenExtFile)
  10      AliGenExtFile::AliGenExtFile()
  11          :AliGenerator(-1)
  12 {
  13     //
  14     fName="ExtFile";
  15     fTitle="Primaries from ext. File";
  16     fFileName="dtujet93.root";
  17     fTreeNtuple=0;
  18     fNcurrent=0;
  19 //
  20 //  Read all particles
  21     fNpart=-1;
  22 }
  23
  24 AliGenExtFile::AliGenExtFile(Int_t npart)
  25     :AliGenerator(npart)
  26 {
  27     //
  28     fName="ExtFile";
  29     fTitle="Primaries from ext. File";
  30     fFileName="dtujet93.root";
  31     fTreeNtuple=0;
  32     fNcurrent=0;
  33 }
  34
  35 //____________________________________________________________
  36 AliGenExtFile::~AliGenExtFile()
  37 {
  38     delete fTreeNtuple;
  39 }
  40
  41 //____________________________________________________________
  42 void AliGenExtFile::NtupleInit()
  43 {
  44 //
  45 // reset the existing file environment and open a new root file if
  46 // the pointer to the Fluka tree is null
  47
  48     TFile *File=0;
  49     if (fTreeNtuple==0) {
  50         if (!File) {
  51             File = new TFile(fFileName);
  52             File->cd();
  53             cout<<"I have opened "<<fFileName<<" file "<<endl;
  54         }
  55 // get the tree address in the Fluka boundary source file
  56         fTreeNtuple = (TTree*)gDirectory->Get("h888");
  57     } else {
  58         File = fTreeNtuple->GetCurrentFile();
  59         File->cd();
  60     }
  61
  62     TTree *h2=fTreeNtuple;
  63 //Set branch addresses
  64 //Set branch addresses
  65     h2->SetBranchAddress("Nihead",&Nihead);
  66     h2->SetBranchAddress("Ihead",Ihead);
  67     h2->SetBranchAddress("Nrhead",&Nrhead);
  68     h2->SetBranchAddress("Rhead",Rhead);
  69     h2->SetBranchAddress("Idpart",&Idpart);
  70     h2->SetBranchAddress("Theta",&Theta);
  71     h2->SetBranchAddress("Phi",&Phi);
  72     h2->SetBranchAddress("P",&P);
  73     h2->SetBranchAddress("E",&E);
  74 }
  75
  76
  77 //____________________________________________________________
  78 void AliGenExtFile::Generate()
  79 {
  80
  81   Float_t polar[3]= {0,0,0};
  82   //
  83   Float_t origin[3]={0,0,0};
  84   Float_t p[3];
  85   Float_t random[6];
  86   Float_t prwn;
  87   Int_t i, j, nt, Ntracks=0;
  88   //
  89   NtupleInit();
  90   TTree *h2=fTreeNtuple;
  91   Int_t nentries = (Int_t) h2->GetEntries();
  92   // loop over number of particles
  93   Int_t nb = (Int_t)h2->GetEvent(fNcurrent);
  94   Int_t i5=Ihead[4];
  95   Int_t i6=Ihead[5];
  96
  97   for (j=0;j<3;j++) origin[j]=fOrigin[j];
  98   if(fVertexSmear==perEvent) {
  99     gMC->Rndm(random,6);
 100     for (j=0;j<3;j++) {
 101         origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 102             TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 103     }
 104   }
 105
 106   if (fNcurrent >= nentries) {
 107       printf("\n No more entries !!! !\n");
 108       return;
 109   }
 110
 111
 112   if (i5==0) {
 113       printf("\n This should never happen !\n");
 114   } else {
 115       printf("\n Next event contains %d tracks! \n", i6);
 116       Ntracks=i6;
 117   }
 118   for (i=0; i<Ntracks; i++) {
 119
 120       Double_t amass = TDatabasePDG::Instance()->GetParticle(Idpart)->Mass();
 121       if(E<=amass) {
 122         Warning("Generate","Particle %d no %d E = %f mass = %f\n",Idpart,i,E,amass);
 123         prwn=0;
 124       } else {
 125         prwn=sqrt((E+amass)*(E-amass));
 126       }
 127
 128       Theta *= TMath::Pi()/180.;
 129       Phi    = (Phi-180)*TMath::Pi()/180.;
 130       if(Theta<fThetaMin || Theta>fThetaMax ||
 131          Phi<fPhiMin || Phi>fPhiMax         ||
 132          prwn<fPMin || prwn>fPMax)
 133       {
 134           ;
 135       } else {
 136           p[0]=prwn*TMath::Sin(Theta)*TMath::Cos(Phi);
 137           p[1]=prwn*TMath::Sin(Theta)*TMath::Sin(Phi);
 138           p[2]=prwn*TMath::Cos(Theta);
 139
 140           if(fVertexSmear==perTrack) {
 141               gMC->Rndm(random,6);
 142               for (j=0;j<3;j++) {
 143                   origin[j]=fOrigin[j]
 144                       +fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 145                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 146               }
 147           }
 148           gAlice->SetTrack(fTrackIt,-1,Idpart,p,origin,polar,0,"Primary",nt);
 149       }
 150       fNcurrent++;
 151       nb = (Int_t)h2->GetEvent(fNcurrent);
 152   }
 153
 154     TFile *File=0;
 155 // Get AliRun object or create it
 156     if (!gAlice) {
 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");
 160     }
 161     TTree *fAli=gAlice->TreeK();
 162     if (fAli) File =fAli->GetCurrentFile();
 163     File->cd();
 164 }
 165
 166
 167
 168
 169
 170
 171
 172
 173