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