[u/mrichter/AliRoot.git] / EVGEN / AliGenFLUKAsource.cxx

#include "AliGenFLUKAsource.h"
#include "AliGenMUONlib.h"
#include "AliMC.h"
#include "AliRun.h"
#include "AliPDG.h"

#include <TDirectory.h>
#include <TFile.h>
#include <TTree.h>
#include <stdlib.h>
 ClassImp(AliGenFLUKAsource)
     AliGenFLUKAsource::AliGenFLUKAsource()
	 :AliGenerator(-1)
{
    //
    fName="FLUKA";
    fTitle="FLUKA Boundary Source";
    // Read in all particle types by default
    fIkine=6;
    // Set maximum admitted age of particles to 1.e-05 by default 
    fAgeMax=1.e-05;
    // Do not add weight
    fAddWeight=1.;
    // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
    // from  specific boundary source to the chamber (fZshift=4.5;),else there 
    // is no need to shift as it reads boundary source for the whole volume of 
    // the Muon Arm; the default value corresponds to boundary source for the
    // whole volume of the MUON Arm 
    fZshift=0;
    // Set the default file 
    fFileName="flukasource.root";

    fTreeFluka=0;
//
//  Read all particles
    fNpart=-1;
    
}

AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
    :AliGenerator(npart)
{
    //
    fName="FLUKA";
    fTitle="FLUKA Boundary Source";
    // Read in all particle types by default
    fIkine=6;
    // Set maximum admitted age of particles to 1.e-05 by default 
    fAgeMax=1.e-05;
    // Do not add weight
    fAddWeight=1.;
    // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
    // from  specific boundary source to the chamber (fZshift=4.5;),else there 
    // is no need to shift as it reads boundary source for the whole volume of 
    // the Muon Arm; the default value corresponds to boundary source for the
    // whole volume of the MUON Arm 
    fZshift=0;
    // Set the default file 
    fFileName="flukasource.root";

    fTreeFluka=0;

}

//____________________________________________________________
AliGenFLUKAsource::~AliGenFLUKAsource()
{
    delete fTreeFluka;
}

//____________________________________________________________
void AliGenFLUKAsource::FlukaInit() 
{
//
// reset the existing file environment and open a new root file if
// the pointer to the Fluka tree is null
    
    TFile *File=0;
    if (fTreeFluka==0) {
        if (!File) {
	    File = new TFile(fFileName);
	    File->cd();
	    cout<<"I have opened "<<fFileName<<" file "<<endl;
        }
// get the tree address in the Fluka boundary source file
	fTreeFluka = (TTree*)gDirectory->Get("h1");
    } else {
        File = fTreeFluka->GetCurrentFile();
        File->cd();
    }

    TTree *h2=fTreeFluka;
//Set branch addresses
    h2->SetBranchAddress("Ip",&Ip);
    h2->SetBranchAddress("Ipp",&Ipp);
    h2->SetBranchAddress("Xi",&Xi);
    h2->SetBranchAddress("Yi",&Yi);
    h2->SetBranchAddress("Zi",&Zi);
    h2->SetBranchAddress("Px",&Px);
    h2->SetBranchAddress("Py",&Py);
    h2->SetBranchAddress("Pz",&Pz);
    h2->SetBranchAddress("Ekin",&Ekin);
    h2->SetBranchAddress("Zv",&Zv);
    h2->SetBranchAddress("Rv",&Rv);
    h2->SetBranchAddress("Itra",&Itra);
    h2->SetBranchAddress("Igas",&Igas);
    h2->SetBranchAddress("Wgt",&Wgt);
    h2->SetBranchAddress("Etag",&Etag);
    h2->SetBranchAddress("Ptg",&Ptg);
    h2->SetBranchAddress("Age",&Age);
}

//____________________________________________________________
void AliGenFLUKAsource::Generate()
{

  const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
			  kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
			  kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
			  kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
			  kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
			  0,kNuMu,kNuMuBar};
  Float_t polar[3]= {0,0,0};
  //
  Float_t origin[3];
  Float_t p[3];
  Float_t prwn;
  Float_t wgt, fwgt;
  Float_t phi;
  char name[100];
  Float_t amass, charge, tlife;
  Int_t itrtyp;
  Int_t iwgt;
  Int_t i, j, part, nt;
  static Int_t irwn;
  //
  Float_t random[2];
  //
  FlukaInit();
  TTree *h2=fTreeFluka;
  Int_t nentries = (Int_t) h2->GetEntries();
  if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
  
  // loop over number of particles
  Int_t nb=0;
  for (i=0; i<fNpart;i++) {
    Int_t ev=gMC->CurrentEvent();
    Int_t entry=fNpart*(ev-1)+i; 
    nb = (Int_t)h2->GetEvent(entry); 
    if (irwn > nentries) {
      printf("No more entries in the FLUKA boundary source file\n");
      TFile *File=0;
      // Get AliRun object or create it 
      if (!gAlice) {
        gAlice = (AliRun*)File->Get("gAlice");
        if (gAlice) printf("AliRun object found on file\n");
        if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
      }
      TTree *fAli=gAlice->TreeK();
      if (fAli) File =fAli->GetCurrentFile();
      File->cd();
      printf("Generate - I'm out \n");
      return;
    }   
    if (Ip > 28 || Ip < 0) {
      irwn++;
      continue;
    }
 
    if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
       irwn++;
       continue;
    } else if (fIkine == 9) {
       if (Ip == 7 || Ip == 8 || Age > fAgeMax) { 
          irwn++;
          continue;
       }
    }
    

    irwn++;
    printf("\n Particle type: %f \n \n ", Ip);
    
    if (Ip ==7){
       prwn=Ekin;
       part=1;
    } else if (Ip == 8) {
       prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
       part=13;
    } else {
       part=ifluge[int(Ip)-1];
       gMC->Gfpart(part, name, itrtyp,  
		   amass, charge, tlife); 
       prwn=sqrt(Ekin*Ekin + 2.*amass);
    }
    origin[0]=Xi;
    origin[1]=Yi;
    origin[2]=Zi;

    p[0]=Px*prwn;
    p[1]=Py*prwn;
    p[2]=Pz*prwn;

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