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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
Revision 1.8  2000/02/14 14:49:38  morsch
Correct particle type for gamma and neutrons
More consistent calculation of momentum from kin. energy and mass

Revision 1.7  1999/11/03 17:43:20  fca
New version from G.Martinez & A.Morsch

Revision 1.6  1999/09/29 09:24:12  fca
Introduction of the Copyright and cvs Log

*/

#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;
    fTreeChain = new TChain("h1");
//
//  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;
    fTreeChain = new TChain("h1"); 
}

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

void AliGenFLUKAsource::AddFile(const Text_t *filename)
{
    fTreeChain->Add(filename);
    
}


//____________________________________________________________
void AliGenFLUKAsource::FlukaInit() 
{
    TChain *h2=fTreeChain;
//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()
{ 
    AliMC* gMC = AliMC::GetMC();

    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=0;
    //
    Float_t random[2];
  //
    FlukaInit();
    TChain *h2=fTreeChain;
    Int_t nentries = (Int_t) h2->GetEntries();
    if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
  
  // loop over number of particles
    Int_t nb=0;
    Int_t ev=gMC->CurrentEvent();
    for (i=0; i<fNpart;i++) {
	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 && fIkine != 10) || Age > fAgeMax){
	    irwn++;
	    continue;
	} else if (fIkine == 9) {
	    if (Ip == 7 || Ip == 8 || Age > fAgeMax) { 
		irwn++;
		continue;
	    }
	} else if (fIkine == 10) {
	    if (Ip == 8 || Age > fAgeMax) { 
		irwn++;
		continue;
	    }
	}
    

	irwn++;
//
// PDG code from FLUKA particle type (Ip)
	part=ifluge[int(Ip)-1];	
//
// Calculate momentum from kinetic energy and mass of the particle
	gMC->Gfpart(part, name, itrtyp,  
		    amass, charge, tlife); 
	prwn=Ekin*sqrt(1. + 2.*amass/Ekin);


	origin[0]=Yi;
	origin[1]=Xi;
	origin[2]=Zi;
	
	p[0]=Py*prwn;
	p[1]=Px*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(fTrackIt,-1,part,p,origin,polar,Age,"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
4c039060	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$
e7bafa6d	18	Revision 1.8 2000/02/14 14:49:38 morsch
	19	Correct particle type for gamma and neutrons
	20	More consistent calculation of momentum from kin. energy and mass
	21
02769dba	22	Revision 1.7 1999/11/03 17:43:20 fca
	23	New version from G.Martinez & A.Morsch
	24
886b6f73	25	Revision 1.6 1999/09/29 09:24:12 fca
	26	Introduction of the Copyright and cvs Log
	27
4c039060	28	*/
4c039060	29
fe4da5cc	30	#include "AliGenFLUKAsource.h"
	31	#include "AliGenMUONlib.h"
	32	#include "AliMC.h"
	33	#include "AliRun.h"
c5970876	34	#include "AliPDG.h"
fe4da5cc	35	#include <TDirectory.h>
	36	#include <TFile.h>
	37	#include <TTree.h>
	38	#include <stdlib.h>
	39	ClassImp(AliGenFLUKAsource)
	40	AliGenFLUKAsource::AliGenFLUKAsource()
	41	:AliGenerator(-1)
	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;
886b6f73	62	fTreeChain = new TChain("h1");
fe4da5cc	63	//
	64	// Read all particles
	65	fNpart=-1;
	66
	67	}
	68
	69	AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
	70	:AliGenerator(npart)
	71	{
	72	//
	73	fName="FLUKA";
	74	fTitle="FLUKA Boundary Source";
	75	// Read in all particle types by default
	76	fIkine=6;
	77	// Set maximum admitted age of particles to 1.e-05 by default
	78	fAgeMax=1.e-05;
	79	// Do not add weight
	80	fAddWeight=1.;
	81	// Shift the z-coordinate of the impact point by 4.5 cm only if it reads
	82	// from specific boundary source to the chamber (fZshift=4.5;),else there
	83	// is no need to shift as it reads boundary source for the whole volume of
	84	// the Muon Arm; the default value corresponds to boundary source for the
	85	// whole volume of the MUON Arm
	86	fZshift=0;
	87	// Set the default file
	88	fFileName="flukasource.root";
	89
	90	fTreeFluka=0;
886b6f73	91	fTreeChain = new TChain("h1");
fe4da5cc	92	}
	93
	94	//____________________________________________________________
	95	AliGenFLUKAsource::~AliGenFLUKAsource()
	96	{
886b6f73	97	if (fTreeFluka) delete fTreeFluka;
	98	if (fTreeChain) delete fTreeChain;
	99	// if (fFileName) delete fFileName;
fe4da5cc	100	}
fe4da5cc	101
886b6f73	102	void AliGenFLUKAsource::AddFile(const Text_t *filename)
fe4da5cc	103	{
886b6f73	104	fTreeChain->Add(filename);
fe4da5cc	105
886b6f73	106	}
886b6f73	107
fe4da5cc	108
886b6f73	109	//____________________________________________________________
	110	void AliGenFLUKAsource::FlukaInit()
	111	{
	112	TChain *h2=fTreeChain;
fe4da5cc	113	//Set branch addresses
	114	h2->SetBranchAddress("Ip",&Ip);
	115	h2->SetBranchAddress("Ipp",&Ipp);
	116	h2->SetBranchAddress("Xi",&Xi);
	117	h2->SetBranchAddress("Yi",&Yi);
	118	h2->SetBranchAddress("Zi",&Zi);
	119	h2->SetBranchAddress("Px",&Px);
	120	h2->SetBranchAddress("Py",&Py);
	121	h2->SetBranchAddress("Pz",&Pz);
	122	h2->SetBranchAddress("Ekin",&Ekin);
	123	h2->SetBranchAddress("Zv",&Zv);
	124	h2->SetBranchAddress("Rv",&Rv);
	125	h2->SetBranchAddress("Itra",&Itra);
	126	h2->SetBranchAddress("Igas",&Igas);
	127	h2->SetBranchAddress("Wgt",&Wgt);
	128	h2->SetBranchAddress("Etag",&Etag);
	129	h2->SetBranchAddress("Ptg",&Ptg);
	130	h2->SetBranchAddress("Age",&Age);
	131	}
	132
	133	//____________________________________________________________
	134	void AliGenFLUKAsource::Generate()
886b6f73	135	{
886b6f73	136	AliMC* gMC = AliMC::GetMC();
fe4da5cc	137
886b6f73	138	const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
1578254f	139	kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
	140	kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
	141	kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
	142	kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
	143	0,kNuMu,kNuMuBar};
886b6f73	144	Float_t polar[3]= {0,0,0};
fe4da5cc	145	//
886b6f73	146	Float_t origin[3];
	147	Float_t p[3];
	148	Float_t prwn;
	149	Float_t wgt, fwgt;
	150	Float_t phi;
	151	char name[100];
	152	Float_t amass, charge, tlife;
	153	Int_t itrtyp;
	154	Int_t iwgt;
	155	Int_t i, j, part, nt;
e7bafa6d	156	static Int_t irwn=0;
886b6f73	157	//
886b6f73	158	Float_t random[2];
fe4da5cc	159	//
886b6f73	160	FlukaInit();
	161	TChain *h2=fTreeChain;
	162	Int_t nentries = (Int_t) h2->GetEntries();
	163	if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
fe4da5cc	164
fe4da5cc	165	// loop over number of particles
886b6f73	166	Int_t nb=0;
cfce8870	167	Int_t ev=gMC->CurrentEvent();
886b6f73	168	for (i=0; i<fNpart;i++) {
	169	Int_t entry=fNpart*(ev-1)+i;
	170	nb = (Int_t)h2->GetEvent(entry);
	171	if (irwn > nentries) {
	172	printf("No more entries in the FLUKA boundary source file\n");
	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	printf("Generate - I'm out \n");
	184	return;
	185	}
	186	if (Ip > 28 \|\| Ip < 0) {
	187	irwn++;
	188	continue;
	189	}
	190
	191	if ((Ip != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) \|\| Age > fAgeMax){
	192	irwn++;
	193	continue;
	194	} else if (fIkine == 9) {
	195	if (Ip == 7 \|\| Ip == 8 \|\| Age > fAgeMax) {
	196	irwn++;
	197	continue;
	198	}
	199	} else if (fIkine == 10) {
	200	if (Ip == 8 \|\| Age > fAgeMax) {
	201	irwn++;
	202	continue;
	203	}
	204	}
fe4da5cc	205
fe4da5cc	206
886b6f73	207	irwn++;
02769dba	208	//
	209	// PDG code from FLUKA particle type (Ip)
	210	part=ifluge[int(Ip)-1];
	211	//
	212	// Calculate momentum from kinetic energy and mass of the particle
	213	gMC->Gfpart(part, name, itrtyp,
	214	amass, charge, tlife);
	215	prwn=Ekinsqrt(1. + 2.amass/Ekin);
	216
	217
886b6f73	218	origin[0]=Yi;
	219	origin[1]=Xi;
	220	origin[2]=Zi;
	221
	222	p[0]=Py*prwn;
	223	p[1]=Px*prwn;
	224	p[2]=Pz*prwn;
fe4da5cc	225
886b6f73	226	//handle particle weight correctly
	227	wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
	228	iwgt=Int_t(wgt);
	229	fwgt=wgt-Float_t(iwgt);
cfce8870	230	gMC->Rndm(random,2);
886b6f73	231	if (random[0] < fwgt) iwgt++;
	232	if (part==1 && iwgt>100) iwgt=100;
	233	Int_t nstack=0;
	234	for (j=0; j<iwgt; j++) {
	235	gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,Age,"Primary",nt);
	236	gMC->Rndm(random,2);
	237	phi=2random[1]TMath::Pi();
	238	Float_t pn1=p[0]TMath::Sin(phi) - p[1]TMath::Cos(phi);
	239	Float_t pn2=p[0]TMath::Cos(phi) + p[1]TMath::Sin(phi);
	240	p[0]=pn1;
	241	p[1]=pn2;
	242	Float_t on1=origin[0]TMath::Sin(phi)-origin[1]TMath::Cos(phi);
	243	Float_t on2=origin[0]TMath::Cos(phi)+origin[1]TMath::Sin(phi);
	244	origin[0]=on1;
	245	origin[1]=on2;
	246	nstack++;
	247	}
	248	if (nstack == 0) continue;
fe4da5cc	249	}
886b6f73	250
fe4da5cc	251	TFile *File=0;
	252	// Get AliRun object or create it
	253	if (!gAlice) {
	254	gAlice = (AliRun*)File->Get("gAlice");
	255	if (gAlice) printf("AliRun object found on file\n");
	256	if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
	257	}
	258	TTree *fAli=gAlice->TreeK();
	259	if (fAli) File =fAli->GetCurrentFile();
	260	File->cd();
	261	}
	262
	263
	264
	265
	266
	267
	268
	269
	270