[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.                  *
 **************************************************************************/

/* $Id$ */

// 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 <stdlib.h>

#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include <RVersion.h>
#include <TChain.h>
#include <TFile.h>
#include <TTree.h>
#include <TVirtualMC.h>

#include "AliRun.h"

ClassImp(AliGenFLUKAsource)

AliGenFLUKAsource::AliGenFLUKAsource()
    :AliGenerator(-1), 
     fIkine(6),
     fAgeMax(1.e-5), 
     fAddWeight(1.),
     fZshift(0.),
     fFrac(0.),
     fSourceId(-1),
     fFileName(0),
     fTreeChain(0),
     fTreeFluka(0),
     fIp(0.),
     fIpp(0.),
     fXi(0.),
     fYi(0.),
     fZi(0.),
     fPx(0.),
     fPy(0.),
     fPz(0.),
     fEkin(0.),
     fZv(0.),
     fRv(0.),
     fItra(0.),
     fIgas(0.),
     fWgt(0.),
     fEtag(0.),
     fPtg(0.),
     fAge(0.)
{
    // Constructor
    fName="FLUKA";
    fTitle="FLUKA Boundary Source";
    // Read in all particle types by default
    // Set maximum admitted age of particles to 1.e-05 by default 
    // Do not add weight
    // 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 
    // Set the default file 
    fTreeChain = new TChain("h1");
//
//  Read all particles
    fNpart=-1;
}

AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
    :AliGenerator(npart), 
     fIkine(6),
     fAgeMax(1.e-5), 
     fAddWeight(1.),
     fZshift(0.),
     fFrac(0.),
     fSourceId(-1),
     fFileName(""),
     fTreeChain(new TChain("h1")),
     fTreeFluka(0),
     fIp(0.),
     fIpp(0.),
     fXi(0.),
     fYi(0.),
     fZi(0.),
     fPx(0.),
     fPy(0.),
     fPz(0.),
     fEkin(0.),
     fZv(0.),
     fRv(0.),
     fItra(0.),
     fIgas(0.),
     fWgt(0.),
     fEtag(0.),
     fPtg(0.),
     fAge(0.)
{
    // Constructor
    fName  = "FLUKA";
    fTitle = "FLUKA Boundary Source";
}

//____________________________________________________________
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 

    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;
    Float_t amass;
    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;
	}   
	
	Int_t ifip = Int_t(fIp);
	

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

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