EVGEN/AliGenFLUKAsource.cxx

   1 #include "AliGenFLUKAsource.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(AliGenFLUKAsource)
  10      AliGenFLUKAsource::AliGenFLUKAsource()
  11          :AliGenerator(-1)
  12 {
  13     //
  14     fName="FLUKA";
  15     fTitle="FLUKA Boundary Source";
  16     // Read in all particle types by default
  17     fIkine=6;
  18     // Set maximum admitted age of particles to 1.e-05 by default
  19     fAgeMax=1.e-05;
  20     // Do not add weight
  21     fAddWeight=1.;
  22     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  23     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  24     // is no need to shift as it reads boundary source for the whole volume of
  25     // the Muon Arm; the default value corresponds to boundary source for the
  26     // whole volume of the MUON Arm
  27     fZshift=0;
  28     // Set the default file
  29     fFileName="flukasource.root";
  30
  31     fTreeFluka=0;
  32 //
  33 //  Read all particles
  34     fNpart=-1;
  35
  36 }
  37
  38 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
  39     :AliGenerator(npart)
  40 {
  41     //
  42     fName="FLUKA";
  43     fTitle="FLUKA Boundary Source";
  44     // Read in all particle types by default
  45     fIkine=6;
  46     // Set maximum admitted age of particles to 1.e-05 by default
  47     fAgeMax=1.e-05;
  48     // Do not add weight
  49     fAddWeight=1.;
  50     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
  51     // from  specific boundary source to the chamber (fZshift=4.5;),else there
  52     // is no need to shift as it reads boundary source for the whole volume of
  53     // the Muon Arm; the default value corresponds to boundary source for the
  54     // whole volume of the MUON Arm
  55     fZshift=0;
  56     // Set the default file
  57     fFileName="flukasource.root";
  58
  59     fTreeFluka=0;
  60
  61 }
  62
  63 //____________________________________________________________
  64 AliGenFLUKAsource::~AliGenFLUKAsource()
  65 {
  66     delete fTreeFluka;
  67 }
  68
  69 //____________________________________________________________
  70 void AliGenFLUKAsource::FlukaInit()
  71 {
  72 //
  73 // reset the existing file environment and open a new root file if
  74 // the pointer to the Fluka tree is null
  75
  76     TFile *File=0;
  77     if (fTreeFluka==0) {
  78         if (!File) {
  79             File = new TFile(fFileName);
  80             File->cd();
  81             cout<<"I have opened "<<fFileName<<" file "<<endl;
  82         }
  83 // get the tree address in the Fluka boundary source file
  84         fTreeFluka = (TTree*)gDirectory->Get("h1");
  85     } else {
  86         File = fTreeFluka->GetCurrentFile();
  87         File->cd();
  88     }
  89
  90     TTree *h2=fTreeFluka;
  91 //Set branch addresses
  92     h2->SetBranchAddress("Ip",&Ip);
  93     h2->SetBranchAddress("Ipp",&Ipp);
  94     h2->SetBranchAddress("Xi",&Xi);
  95     h2->SetBranchAddress("Yi",&Yi);
  96     h2->SetBranchAddress("Zi",&Zi);
  97     h2->SetBranchAddress("Px",&Px);
  98     h2->SetBranchAddress("Py",&Py);
  99     h2->SetBranchAddress("Pz",&Pz);
 100     h2->SetBranchAddress("Ekin",&Ekin);
 101     h2->SetBranchAddress("Zv",&Zv);
 102     h2->SetBranchAddress("Rv",&Rv);
 103     h2->SetBranchAddress("Itra",&Itra);
 104     h2->SetBranchAddress("Igas",&Igas);
 105     h2->SetBranchAddress("Wgt",&Wgt);
 106     h2->SetBranchAddress("Etag",&Etag);
 107     h2->SetBranchAddress("Ptg",&Ptg);
 108     h2->SetBranchAddress("Age",&Age);
 109 }
 110
 111 //____________________________________________________________
 112 void AliGenFLUKAsource::Generate()
 113 {
 114
 115   AliMC* pMC = AliMC::GetMC();
 116
 117   const Int_t ifluge[28]={14, 15, 3, 2, 4, 4, 1, 13, 25, 5, 6, 10, 8, 9,
 118                          11, 12,18, 26, 16, 21, 19, 20, 7, 16,16,0,0,0};
 119   Float_t polar[3]= {0,0,0};
 120   //
 121   Float_t origin[3];
 122   Float_t p[3];
 123   Float_t prwn;
 124   Float_t wgt, fwgt;
 125   Float_t phi;
 126   char name[100];
 127   Float_t amass, charge, tlife;
 128   Int_t itrtyp;
 129   Int_t iwgt;
 130   Int_t i, j, part, nt;
 131   static Int_t irwn;
 132   //
 133   Float_t random[2];
 134   //
 135   FlukaInit();
 136   TTree *h2=fTreeFluka;
 137   Int_t nentries = (Int_t) h2->GetEntries();
 138   if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
 139
 140   // loop over number of particles
 141   Int_t nb=0;
 142   for (i=0; i<fNpart;i++) {
 143     Int_t ev=pMC->CurrentEvent();
 144     Int_t entry=fNpart*(ev-1)+i;
 145     nb = (Int_t)h2->GetEvent(entry);
 146     if (irwn > nentries) {
 147       printf("No more entries in the FLUKA boundary source file\n");
 148       TFile *File=0;
 149       // Get AliRun object or create it
 150       if (!gAlice) {
 151         gAlice = (AliRun*)File->Get("gAlice");
 152         if (gAlice) printf("AliRun object found on file\n");
 153         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 154       }
 155       TTree *fAli=gAlice->TreeK();
 156       if (fAli) File =fAli->GetCurrentFile();
 157       File->cd();
 158       printf("Generate - I'm out \n");
 159       return;
 160     }
 161     if (Ip > 28 || Ip < 0) {
 162       irwn++;
 163       continue;
 164     }
 165
 166     if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
 167        irwn++;
 168        continue;
 169     } else if (fIkine == 9) {
 170        if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
 171           irwn++;
 172           continue;
 173        }
 174     }
 175
 176
 177     irwn++;
 178     printf("\n Particle type: %f \n \n ", Ip);
 179
 180     if (Ip ==7){
 181        prwn=Ekin;
 182        part=1;
 183     } else if (Ip == 8) {
 184        prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
 185        part=13;
 186     } else {
 187        part=ifluge[int(Ip)-1];
 188        pMC->Gfpart(part, name, itrtyp,
 189                    amass, charge, tlife);
 190        prwn=sqrt(Ekin*Ekin + 2.*amass);
 191     }
 192     origin[0]=Xi;
 193     origin[1]=Yi;
 194     origin[2]=Zi;
 195
 196     p[0]=Px*prwn;
 197     p[1]=Py*prwn;
 198     p[2]=Pz*prwn;
 199
 200     //handle particle weight correctly
 201     wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
 202     iwgt=Int_t(wgt);
 203     fwgt=wgt-Float_t(iwgt);
 204     pMC->Rndm(random,2);
 205     if (random[0] < fwgt) iwgt++;
 206     if (part==1 && iwgt>100) iwgt=100;
 207     Int_t nstack=0;
 208     for (j=0; j<iwgt; j++) {
 209         gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
 210         pMC->Rndm(random,2);
 211         phi=2*random[1]*TMath::Pi();
 212         Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
 213         Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
 214         p[0]=pn1;
 215         p[1]=pn2;
 216         Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
 217         Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
 218         origin[0]=on1;
 219         origin[1]=on2;
 220         nstack++;
 221     }
 222     if (nstack == 0) continue;
 223   }
 224
 225     TFile *File=0;
 226 // Get AliRun object or create it
 227     if (!gAlice) {
 228         gAlice = (AliRun*)File->Get("gAlice");
 229         if (gAlice) printf("AliRun object found on file\n");
 230         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 231     }
 232     TTree *fAli=gAlice->TreeK();
 233     if (fAli) File =fAli->GetCurrentFile();
 234     File->cd();
 235 }
 236
 237
 238
 239
 240
 241
 242
 243
 244