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

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