EVGEN/AliGenExtFile.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   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  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.6  1999/11/09 07:38:48  fca
  19 Changes for compatibility with version 2.23 of ROOT
  20
  21 Revision 1.5  1999/09/29 09:24:12  fca
  22 Introduction of the Copyright and cvs Log
  23
  24 */
  25
  26 #include <iostream.h>
  27
  28 #include "AliGenExtFile.h"
  29 #include "AliGenMUONlib.h"
  30 #include "AliMC.h"
  31 #include "AliRun.h"
  32 #include <TDirectory.h>
  33 #include <TFile.h>
  34 #include <TTree.h>
  35 #include <stdlib.h>
  36  ClassImp(AliGenExtFile)
  37      AliGenExtFile::AliGenExtFile()
  38          :AliGenerator(-1)
  39 {
  40     //
  41     fName="ExtFile";
  42     fTitle="Primaries from ext. File";
  43     fFileName="dtujet93.root";
  44     fTreeNtuple=0;
  45     fNcurrent=0;
  46 //
  47 //  Read all particles
  48     fNpart=-1;
  49 }
  50
  51 AliGenExtFile::AliGenExtFile(Int_t npart)
  52     :AliGenerator(npart)
  53 {
  54     //
  55     fName="ExtFile";
  56     fTitle="Primaries from ext. File";
  57     fFileName="dtujet93.root";
  58     fTreeNtuple=0;
  59     fNcurrent=0;
  60 }
  61
  62 //____________________________________________________________
  63 AliGenExtFile::~AliGenExtFile()
  64 {
  65     delete fTreeNtuple;
  66 }
  67
  68 //____________________________________________________________
  69 void AliGenExtFile::NtupleInit()
  70 {
  71 //
  72 // reset the existing file environment and open a new root file if
  73 // the pointer to the Fluka tree is null
  74
  75     TFile *File=0;
  76     if (fTreeNtuple==0) {
  77         if (!File) {
  78             File = new TFile(fFileName);
  79             File->cd();
  80             cout<<"I have opened "<<fFileName<<" file "<<endl;
  81         }
  82 // get the tree address in the Fluka boundary source file
  83         fTreeNtuple = (TTree*)gDirectory->Get("h888");
  84     } else {
  85         File = fTreeNtuple->GetCurrentFile();
  86         File->cd();
  87     }
  88
  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);
 101 }
 102
 103
 104 //____________________________________________________________
 105 void AliGenExtFile::Generate()
 106 {
 107
 108   Float_t polar[3]= {0,0,0};
 109   //
 110   Float_t origin[3]={0,0,0};
 111   Float_t p[3];
 112   Float_t random[6];
 113   Float_t prwn;
 114   Int_t i, j, nt, Ntracks=0;
 115   //
 116   NtupleInit();
 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);
 121   Int_t i5=Ihead[4];
 122   Int_t i6=Ihead[5];
 123
 124   for (j=0;j<3;j++) origin[j]=fOrigin[j];
 125   if(fVertexSmear==perEvent) {
 126     gMC->Rndm(random,6);
 127     for (j=0;j<3;j++) {
 128         origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 129             TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 130     }
 131   }
 132
 133   if (fNcurrent >= nentries) {
 134       printf("\n No more entries !!! !\n");
 135       return;
 136   }
 137
 138
 139   if (i5==0) {
 140       printf("\n This should never happen !\n");
 141   } else {
 142       printf("\n Next event contains %d tracks! \n", i6);
 143       Ntracks=i6;
 144   }
 145   for (i=0; i<Ntracks; i++) {
 146       Idpart=gMC->PDGFromId(Idpart);
 147       Double_t amass = TDatabasePDG::Instance()->GetParticle(Idpart)->Mass();
 148       if(E<=amass) {
 149         Warning("Generate","Particle %d no %d E = %f mass = %f\n",
 150                 Idpart,i,E,amass);
 151         prwn=0;
 152       } else {
 153         prwn=sqrt((E+amass)*(E-amass));
 154       }
 155
 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)
 161       {
 162           ;
 163       } else {
 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);
 167
 168           if(fVertexSmear==perTrack) {
 169               gMC->Rndm(random,6);
 170               for (j=0;j<3;j++) {
 171                   origin[j]=fOrigin[j]
 172                       +fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 173                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 174               }
 175           }
 176           gAlice->SetTrack(fTrackIt,-1,Idpart,p,origin,polar,0,"Primary",nt);
 177       }
 178       fNcurrent++;
 179       nb = (Int_t)h2->GetEvent(fNcurrent);
 180   }
 181
 182     TFile *File=0;
 183 // Get AliRun object or create it
 184     if (!gAlice) {
 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");
 188     }
 189     TTree *fAli=gAlice->TreeK();
 190     if (fAli) File =fAli->GetCurrentFile();
 191     File->cd();
 192 }
 193
 194
 195
 196
 197
 198
 199
 200
 201