[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)
{
    // 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="";

    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="";

    fTreeFluka=0;
    fTreeChain = new TChain("h1"); 
    fSourceId=-1;
}

AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource):
    AliGenerator(FLUKAsource)
{
// Copy constructor
    FLUKAsource.Copy(*this);
}


//____________________________________________________________
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();
}


AliGenFLUKAsource& AliGenFLUKAsource::operator=(const  AliGenFLUKAsource& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return (*this);
}


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