[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.12  2000/11/30 07:12:50  alibrary
Introducing new Rndm and QA classes

Revision 1.11  2000/06/14 15:20:40  morsch
Include clean-up (IH)

Revision 1.10  2000/06/09 20:31:34  morsch
All coding rule violations except RS3 corrected

Revision 1.9  2000/03/07 13:52:54  morsch
static Int_t irwn=0;

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

*/


// Read background particles from a FLUKA boundary source file
// This is a very special generator that works for background studies for the muon-spectrometer.
// The input files come from FLUKA simulations.
// Files can be chained. 
// Author: andreas.morsch@cern.ch

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


#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <stdlib.h>
 ClassImp(AliGenFLUKAsource)
     AliGenFLUKAsource::AliGenFLUKAsource()
	 :AliGenerator(-1)
{
    // Constructor
    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)
{
    // Constructor
    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"); 
    fSourceId=-1;
}

AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource)
{
// copy constructor
}


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

void AliGenFLUKAsource::AddFile(const Text_t *filename)
{
// Add a file to the chain
    fTreeChain->Add(filename);
    
}


//____________________________________________________________
void AliGenFLUKAsource::FlukaInit() 
{
// Set branch addresses of data entries
    TChain *h2=fTreeChain;
//Set branch addresses
    h2->SetBranchAddress("Ip",&fIp);
    h2->SetBranchAddress("Ipp",&fIpp);
    h2->SetBranchAddress("Xi",&fXi);
    h2->SetBranchAddress("Yi",&fYi);
    h2->SetBranchAddress("Zi",&fZi);
    h2->SetBranchAddress("Px",&fPx);
    h2->SetBranchAddress("Py",&fPy);
    h2->SetBranchAddress("Pz",&fPz);
    h2->SetBranchAddress("Ekin",&fEkin);
    h2->SetBranchAddress("Zv",&fZv);
    h2->SetBranchAddress("Rv",&fRv);
    h2->SetBranchAddress("Itra",&fItra);
    h2->SetBranchAddress("Igas",&fIgas);
    h2->SetBranchAddress("Wgt",&fWgt);
    h2->SetBranchAddress("Etag",&fEtag);
    h2->SetBranchAddress("Ptg",&fPtg);
    h2->SetBranchAddress("Age",&fAge);
}

//____________________________________________________________
void AliGenFLUKAsource::Generate()
{
// Generate one event 
    AliMC* gMC = AliMC::GetMC();

    const Int_t kIfluge[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 *pFile=0;
	    // Get AliRun object or create it 
	    if (!gAlice) {
		gAlice = (AliRun*)pFile->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) pFile =fAli->GetCurrentFile();
	    pFile->cd();
	    printf("Generate - I'm out \n");
	    return;
	}   

	if (fSourceId != -1 && fIgas !=fSourceId) {
	    irwn++;
	    continue;
	}
	
	if (fIp > 28 || fIp < 0) {
	    irwn++;
	    continue;
	}
	
	if ((fIp != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) || fAge > fAgeMax){
	    irwn++;
	    continue;
	} else if (fIkine == 9) {
	    if (fIp == 7 || fIp == 8 || fAge > fAgeMax) { 
		irwn++;
		continue;
	    }
	} else if (fIkine == 10) {
	    if (fIp == 8 || fAge > fAgeMax) { 
		irwn++;
		continue;
	    }
	}
    

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


	origin[0]=fYi;
	origin[1]=fXi;
	origin[2]=fZi;
	
	p[0]=fPy*prwn;
	p[1]=fPx*prwn;
	p[2]=fPz*prwn;

	//handle particle weight correctly
	wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
	iwgt=Int_t(wgt);
	fwgt=wgt-Float_t(iwgt);
	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,fAge,kPPrimary,nt);
	    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 *pFile=0;
// Get AliRun object or create it 
    if (!gAlice) {
	gAlice = (AliRun*)pFile->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) pFile =fAli->GetCurrentFile();
    pFile->cd();
}


AliGenFLUKAsource& AliGenFLUKAsource::operator=(const  AliGenFLUKAsource& rhs)
{
// Assignment operator
    return *this;
}
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$
675e9664	18	Revision 1.12 2000/11/30 07:12:50 alibrary
	19	Introducing new Rndm and QA classes
	20
65fb704d	21	Revision 1.11 2000/06/14 15:20:40 morsch
	22	Include clean-up (IH)
	23
5c3fd7ea	24	Revision 1.10 2000/06/09 20:31:34 morsch
	25	All coding rule violations except RS3 corrected
	26
f87cfe57	27	Revision 1.9 2000/03/07 13:52:54 morsch
	28	static Int_t irwn=0;
	29
e7bafa6d	30	Revision 1.8 2000/02/14 14:49:38 morsch
	31	Correct particle type for gamma and neutrons
	32	More consistent calculation of momentum from kin. energy and mass
	33
02769dba	34	Revision 1.7 1999/11/03 17:43:20 fca
	35	New version from G.Martinez & A.Morsch
	36
886b6f73	37	Revision 1.6 1999/09/29 09:24:12 fca
	38	Introduction of the Copyright and cvs Log
	39
4c039060	40	*/
4c039060	41
675e9664	42
	43
	44	// Read background particles from a FLUKA boundary source file
	45	// This is a very special generator that works for background studies for the muon-spectrometer.
	46	// The input files come from FLUKA simulations.
	47	// Files can be chained.
	48	// Author: andreas.morsch@cern.ch
	49
fe4da5cc	50	#include "AliGenFLUKAsource.h"
fe4da5cc	51	#include "AliMC.h"
fe4da5cc	52	#include "AliRun.h"
c5970876	53	#include "AliPDG.h"
5c3fd7ea	54
5c3fd7ea	55
fe4da5cc	56	#include <TFile.h>
fe4da5cc	57	#include <TTree.h>
f87cfe57	58	#include <TChain.h>
fe4da5cc	59	#include <stdlib.h>
	60	ClassImp(AliGenFLUKAsource)
	61	AliGenFLUKAsource::AliGenFLUKAsource()
	62	:AliGenerator(-1)
	63	{
f87cfe57	64	// Constructor
fe4da5cc	65	fName="FLUKA";
	66	fTitle="FLUKA Boundary Source";
	67	// Read in all particle types by default
	68	fIkine=6;
	69	// Set maximum admitted age of particles to 1.e-05 by default
	70	fAgeMax=1.e-05;
	71	// Do not add weight
	72	fAddWeight=1.;
	73	// Shift the z-coordinate of the impact point by 4.5 cm only if it reads
	74	// from specific boundary source to the chamber (fZshift=4.5;),else there
	75	// is no need to shift as it reads boundary source for the whole volume of
	76	// the Muon Arm; the default value corresponds to boundary source for the
	77	// whole volume of the MUON Arm
	78	fZshift=0;
	79	// Set the default file
	80	fFileName="flukasource.root";
	81
	82	fTreeFluka=0;
886b6f73	83	fTreeChain = new TChain("h1");
fe4da5cc	84	//
	85	// Read all particles
	86	fNpart=-1;
f87cfe57	87
fe4da5cc	88
	89	}
	90
	91	AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
	92	:AliGenerator(npart)
	93	{
f87cfe57	94	// Constructor
675e9664	95	fName = "FLUKA";
675e9664	96	fTitle = "FLUKA Boundary Source";
fe4da5cc	97	// Read in all particle types by default
	98	fIkine=6;
	99	// Set maximum admitted age of particles to 1.e-05 by default
	100	fAgeMax=1.e-05;
	101	// Do not add weight
	102	fAddWeight=1.;
	103	// Shift the z-coordinate of the impact point by 4.5 cm only if it reads
	104	// from specific boundary source to the chamber (fZshift=4.5;),else there
	105	// is no need to shift as it reads boundary source for the whole volume of
	106	// the Muon Arm; the default value corresponds to boundary source for the
	107	// whole volume of the MUON Arm
	108	fZshift=0;
	109	// Set the default file
	110	fFileName="flukasource.root";
	111
	112	fTreeFluka=0;
886b6f73	113	fTreeChain = new TChain("h1");
f87cfe57	114	fSourceId=-1;
fe4da5cc	115	}
fe4da5cc	116
f87cfe57	117	AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource)
	118	{
	119	// copy constructor
	120	}
	121
	122
fe4da5cc	123	//____________________________________________________________
	124	AliGenFLUKAsource::~AliGenFLUKAsource()
	125	{
f87cfe57	126	// Destructor
886b6f73	127	if (fTreeFluka) delete fTreeFluka;
886b6f73	128	if (fTreeChain) delete fTreeChain;
fe4da5cc	129	}
fe4da5cc	130
886b6f73	131	void AliGenFLUKAsource::AddFile(const Text_t *filename)
fe4da5cc	132	{
f87cfe57	133	// Add a file to the chain
886b6f73	134	fTreeChain->Add(filename);
fe4da5cc	135
886b6f73	136	}
886b6f73	137
fe4da5cc	138
886b6f73	139	//____________________________________________________________
	140	void AliGenFLUKAsource::FlukaInit()
	141	{
f87cfe57	142	// Set branch addresses of data entries
886b6f73	143	TChain *h2=fTreeChain;
fe4da5cc	144	//Set branch addresses
f87cfe57	145	h2->SetBranchAddress("Ip",&fIp);
	146	h2->SetBranchAddress("Ipp",&fIpp);
	147	h2->SetBranchAddress("Xi",&fXi);
	148	h2->SetBranchAddress("Yi",&fYi);
	149	h2->SetBranchAddress("Zi",&fZi);
	150	h2->SetBranchAddress("Px",&fPx);
	151	h2->SetBranchAddress("Py",&fPy);
	152	h2->SetBranchAddress("Pz",&fPz);
	153	h2->SetBranchAddress("Ekin",&fEkin);
	154	h2->SetBranchAddress("Zv",&fZv);
	155	h2->SetBranchAddress("Rv",&fRv);
	156	h2->SetBranchAddress("Itra",&fItra);
	157	h2->SetBranchAddress("Igas",&fIgas);
	158	h2->SetBranchAddress("Wgt",&fWgt);
	159	h2->SetBranchAddress("Etag",&fEtag);
	160	h2->SetBranchAddress("Ptg",&fPtg);
	161	h2->SetBranchAddress("Age",&fAge);
fe4da5cc	162	}
	163
	164	//____________________________________________________________
	165	void AliGenFLUKAsource::Generate()
f87cfe57	166	{
f87cfe57	167	// Generate one event
886b6f73	168	AliMC* gMC = AliMC::GetMC();
fe4da5cc	169
f87cfe57	170	const Int_t kIfluge[28]={kProton, kProtonBar, kElectron, kPositron,
1578254f	171	kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
	172	kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
	173	kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
	174	kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
	175	0,kNuMu,kNuMuBar};
886b6f73	176	Float_t polar[3]= {0,0,0};
fe4da5cc	177	//
886b6f73	178	Float_t origin[3];
	179	Float_t p[3];
	180	Float_t prwn;
	181	Float_t wgt, fwgt;
	182	Float_t phi;
	183	char name[100];
	184	Float_t amass, charge, tlife;
	185	Int_t itrtyp;
	186	Int_t iwgt;
	187	Int_t i, j, part, nt;
e7bafa6d	188	static Int_t irwn=0;
886b6f73	189	//
886b6f73	190	Float_t random[2];
fe4da5cc	191	//
886b6f73	192	FlukaInit();
	193	TChain *h2=fTreeChain;
	194	Int_t nentries = (Int_t) h2->GetEntries();
	195	if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
fe4da5cc	196
fe4da5cc	197	// loop over number of particles
886b6f73	198	Int_t nb=0;
cfce8870	199	Int_t ev=gMC->CurrentEvent();
886b6f73	200	for (i=0; i<fNpart;i++) {
	201	Int_t entry=fNpart*(ev-1)+i;
	202	nb = (Int_t)h2->GetEvent(entry);
	203	if (irwn > nentries) {
	204	printf("No more entries in the FLUKA boundary source file\n");
f87cfe57	205	TFile *pFile=0;
886b6f73	206	// Get AliRun object or create it
886b6f73	207	if (!gAlice) {
f87cfe57	208	gAlice = (AliRun*)pFile->Get("gAlice");
886b6f73	209	if (gAlice) printf("AliRun object found on file\n");
	210	if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
	211	}
	212	TTree *fAli=gAlice->TreeK();
f87cfe57	213	if (fAli) pFile =fAli->GetCurrentFile();
f87cfe57	214	pFile->cd();
886b6f73	215	printf("Generate - I'm out \n");
	216	return;
	217	}
f87cfe57	218
	219	if (fSourceId != -1 && fIgas !=fSourceId) {
	220	irwn++;
	221	continue;
	222	}
	223
	224	if (fIp > 28 \|\| fIp < 0) {
886b6f73	225	irwn++;
	226	continue;
	227	}
	228
f87cfe57	229	if ((fIp != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) \|\| fAge > fAgeMax){
886b6f73	230	irwn++;
	231	continue;
	232	} else if (fIkine == 9) {
f87cfe57	233	if (fIp == 7 \|\| fIp == 8 \|\| fAge > fAgeMax) {
886b6f73	234	irwn++;
	235	continue;
	236	}
	237	} else if (fIkine == 10) {
f87cfe57	238	if (fIp == 8 \|\| fAge > fAgeMax) {
886b6f73	239	irwn++;
	240	continue;
	241	}
	242	}
fe4da5cc	243
fe4da5cc	244
886b6f73	245	irwn++;
02769dba	246	//
f87cfe57	247	// PDG code from FLUKA particle type (fIp)
f87cfe57	248	part=kIfluge[int(fIp)-1];
02769dba	249	//
	250	// Calculate momentum from kinetic energy and mass of the particle
	251	gMC->Gfpart(part, name, itrtyp,
	252	amass, charge, tlife);
f87cfe57	253	prwn=fEkinsqrt(1. + 2.amass/fEkin);
02769dba	254
02769dba	255
f87cfe57	256	origin[0]=fYi;
	257	origin[1]=fXi;
	258	origin[2]=fZi;
886b6f73	259
f87cfe57	260	p[0]=fPy*prwn;
	261	p[1]=fPx*prwn;
	262	p[2]=fPz*prwn;
fe4da5cc	263
886b6f73	264	//handle particle weight correctly
f87cfe57	265	wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
886b6f73	266	iwgt=Int_t(wgt);
886b6f73	267	fwgt=wgt-Float_t(iwgt);
65fb704d	268	Rndm(random,2);
886b6f73	269	if (random[0] < fwgt) iwgt++;
	270	if (part==1 && iwgt>100) iwgt=100;
	271	Int_t nstack=0;
	272	for (j=0; j<iwgt; j++) {
65fb704d	273	gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,fAge,kPPrimary,nt);
65fb704d	274	Rndm(random,2);
886b6f73	275	phi=2random[1]TMath::Pi();
	276	Float_t pn1=p[0]TMath::Sin(phi) - p[1]TMath::Cos(phi);
	277	Float_t pn2=p[0]TMath::Cos(phi) + p[1]TMath::Sin(phi);
	278	p[0]=pn1;
	279	p[1]=pn2;
	280	Float_t on1=origin[0]TMath::Sin(phi)-origin[1]TMath::Cos(phi);
	281	Float_t on2=origin[0]TMath::Cos(phi)+origin[1]TMath::Sin(phi);
	282	origin[0]=on1;
	283	origin[1]=on2;
	284	nstack++;
	285	}
	286	if (nstack == 0) continue;
fe4da5cc	287	}
886b6f73	288
f87cfe57	289	TFile *pFile=0;
fe4da5cc	290	// Get AliRun object or create it
fe4da5cc	291	if (!gAlice) {
f87cfe57	292	gAlice = (AliRun*)pFile->Get("gAlice");
fe4da5cc	293	if (gAlice) printf("AliRun object found on file\n");
	294	if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
	295	}
	296	TTree *fAli=gAlice->TreeK();
f87cfe57	297	if (fAli) pFile =fAli->GetCurrentFile();
f87cfe57	298	pFile->cd();
fe4da5cc	299	}
	300
	301
f87cfe57	302	AliGenFLUKAsource& AliGenFLUKAsource::operator=(const AliGenFLUKAsource& rhs)
	303	{
	304	// Assignment operator
	305	return *this;
	306	}
fe4da5cc	307
	308
	309
	310
	311
	312