[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEMD.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.                  *
 **************************************************************************/


// Class to read events from external (TNtupla) file
// Events -> neutron removal by EM dissociation of Pb nuclei
// Data from RELDIS code (by I. Pshenichov)

#include <TFile.h>
#include <TParticle.h>
#include <TTree.h>
#include <TVirtualMC.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>

#include "AliGenReaderEMD.h"

ClassImp(AliGenReaderEMD)

AliGenReaderEMD::AliGenReaderEMD():
    fStartEvent(0),
    fNcurrent(0),  
    fNparticle(0), 
    fTreeNtuple(0),
    fIPSide(0),
    fPcToTrack(0),
    fNnLeft(0),
    fEnLeft(0),
    fNnRight(0),
    fEnRight(0),
    fNpLeft(0),
    fEtapLeft(0),
    fNpRight(0),
    fEtapRight(0)
{
// Default constructor
}

AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
    AliGenReader(reader),
    fStartEvent(0),
    fNcurrent(0),  
    fNparticle(0), 
    fTreeNtuple(0),
    fIPSide(0),
    fPcToTrack(0),
    fNnLeft(0),
    fEnLeft(0),
    fNnRight(0),
    fEnRight(0),
    fNpLeft(0),
    fEtapLeft(0),
    fNpRight(0),
    fEtapRight(0)
{
    // Copy Constructor
    reader.Copy(*this);
}
  // -----------------------------------------------------------------------------------
AliGenReaderEMD::~AliGenReaderEMD()
{
    delete fTreeNtuple;
}

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

// -----------------------------------------------------------------------------------
void AliGenReaderEMD::Copy(TObject&) const
{
    //
    // Copy 
    //
    Fatal("Copy","Not implemented!\n");
}

// -----------------------------------------------------------------------------------
void AliGenReaderEMD::Init() 
{
//
// Reset the existing file environment and open a new root file
    
    TFile *pFile=0;
    if (!pFile) {
	pFile = new TFile(fFileName);
	pFile->cd();
	printf("\n %s file opened\n\n", fFileName);
    }
    fTreeNtuple = (TTree*)gDirectory->Get("h2032");
    fNcurrent = fStartEvent;

    TTree *Ntu=fTreeNtuple;
    //
    // Set branch addresses
    // **** neutrons
    Ntu->SetBranchAddress("Nleft",&fNnLeft);
    Ntu->SetBranchAddress("Eleft",&fEnLeft);
    Ntu->SetBranchAddress("Pxl",  fPxnLeft);
    Ntu->SetBranchAddress("Pyl",  fPynLeft);
    Ntu->SetBranchAddress("Pzl",  fPznLeft);
    Ntu->SetBranchAddress("Nright",&fNnRight);
    Ntu->SetBranchAddress("Eright",&fEnRight);
    Ntu->SetBranchAddress("Pxr",   fPxnRight);
    Ntu->SetBranchAddress("Pyr",   fPynRight);
    Ntu->SetBranchAddress("Pzr",   fPznRight);
    // **** protons
    Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
    Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
    Ntu->SetBranchAddress("Pxl_p",  fPxpLeft);
    Ntu->SetBranchAddress("Pyl_p",  fPypLeft);
    Ntu->SetBranchAddress("Pzl_p",  fPzpLeft);
    Ntu->SetBranchAddress("Nright_p",&fNpRight);
    Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
    Ntu->SetBranchAddress("Pxr_p",   fPxpRight);
    Ntu->SetBranchAddress("Pyr_p",   fPypRight);
    Ntu->SetBranchAddress("Pzr_p",   fPzpRight);
}

// -----------------------------------------------------------------------------------
Int_t AliGenReaderEMD::NextEvent() 
{
    // Read the next event  
    Int_t nTracks=0;
    
    TFile* pFile = fTreeNtuple->GetCurrentFile();
    pFile->cd();
    

    Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
    if(fNcurrent < nentries) {
	fTreeNtuple->GetEvent(fNcurrent);
	printf("\n *** Event %d read ***\n\n",fNcurrent);
	fNcurrent++;
	//printf("\n \t \t LEFT-> %d neutrons and %d protons emitted\n", fNnLeft, fNpLeft);
	//printf("\t \t RIGHT-> %d neutrons and %d protons emitted\n\n", fNnRight, fNpRight);
	//
	// **** IPSIde 		=0->RIGHT side, =1->LEFT side  
	// #### fPcToTrack 	=0->neutrons, =1->protons
	if(fIPSide==0){
	  if(fPcToTrack==0){
	    printf("\n \t \t Tracking %d neutrons emitted on RIGHT side\n\n", fNnRight);
	    nTracks    = fNnRight;
	  }
	  else if(fPcToTrack==1){
	    printf("\n \t \t Tracking %d protons emitted on RIGHT side\n\n", fNpRight);
	    nTracks    = fNpRight;
	  }
	}
	else if(fIPSide==1){
	  if(fPcToTrack==0){
	    printf("\n \t \t Tracking %d neutrons emitted on LEFT side\n", fNnLeft);
	    nTracks    = fNnLeft;
	  }
	  else if(fPcToTrack==1){
	    printf("\n \t \t Tracking %d protons emitted on LEFT side\n", fNpLeft);
	    nTracks    = fNpLeft;
	  }
	}
	fNparticle = 0;
	return nTracks;
    }

    return 0;
}

// -----------------------------------------------------------------------------------
TParticle* AliGenReaderEMD::NextParticle() 
{
    // Read the next particle
    Float_t p[4];
    Int_t ipart = kNeutron;
    Double_t amass = TDatabasePDG::Instance()->GetParticle(kNeutron)->Mass();
    p[0] = fPxnRight[fNparticle];
    p[1] = fPynRight[fNparticle];
    p[2] = fPznRight[fNparticle];
    Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
    p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
    //printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
    
    if(p[3]<=amass) 
      Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);

    TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1, 
    	p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
    //fNcurrent++;
    fNparticle++;
    return particle;
}
Commit	Line	Data
f1d1affa	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	// Class to read events from external (TNtupla) file
	18	// Events -> neutron removal by EM dissociation of Pb nuclei
	19	// Data from RELDIS code (by I. Pshenichov)
	20
	21	#include <TFile.h>
	22	#include <TParticle.h>
	23	#include <TTree.h>
	24	#include <TVirtualMC.h>
	25	#include <TDatabasePDG.h>
	26	#include <TPDGCode.h>
	27
	28	#include "AliGenReaderEMD.h"
	29
706938e6	30	ClassImp(AliGenReaderEMD)
f1d1affa	31
1c56e311	32	AliGenReaderEMD::AliGenReaderEMD():
	33	fStartEvent(0),
	34	fNcurrent(0),
	35	fNparticle(0),
	36	fTreeNtuple(0),
	37	fIPSide(0),
	38	fPcToTrack(0),
	39	fNnLeft(0),
	40	fEnLeft(0),
	41	fNnRight(0),
	42	fEnRight(0),
	43	fNpLeft(0),
	44	fEtapLeft(0),
	45	fNpRight(0),
	46	fEtapRight(0)
f1d1affa	47	{
f1d1affa	48	// Default constructor
f1d1affa	49	}
f1d1affa	50
1c56e311	51	AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
	52	AliGenReader(reader),
	53	fStartEvent(0),
	54	fNcurrent(0),
	55	fNparticle(0),
	56	fTreeNtuple(0),
	57	fIPSide(0),
	58	fPcToTrack(0),
	59	fNnLeft(0),
	60	fEnLeft(0),
	61	fNnRight(0),
	62	fEnRight(0),
	63	fNpLeft(0),
	64	fEtapLeft(0),
	65	fNpRight(0),
	66	fEtapRight(0)
	67	{
	68	// Copy Constructor
	69	reader.Copy(*this);
	70	}
f1d1affa	71	// -----------------------------------------------------------------------------------
	72	AliGenReaderEMD::~AliGenReaderEMD()
	73	{
	74	delete fTreeNtuple;
	75	}
	76
	77	// -----------------------------------------------------------------------------------
	78	AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
	79	{
	80	// Assignment operator
	81	rhs.Copy(*this);
	82	return *this;
	83	}
	84
	85	// -----------------------------------------------------------------------------------
	86	void AliGenReaderEMD::Copy(TObject&) const
	87	{
	88	//
	89	// Copy
	90	//
	91	Fatal("Copy","Not implemented!\n");
	92	}
	93
	94	// -----------------------------------------------------------------------------------
	95	void AliGenReaderEMD::Init()
	96	{
	97	//
	98	// Reset the existing file environment and open a new root file
	99
	100	TFile *pFile=0;
	101	if (!pFile) {
	102	pFile = new TFile(fFileName);
	103	pFile->cd();
	104	printf("\n %s file opened\n\n", fFileName);
	105	}
	106	fTreeNtuple = (TTree*)gDirectory->Get("h2032");
	107	fNcurrent = fStartEvent;
	108
	109	TTree *Ntu=fTreeNtuple;
	110	//
	111	// Set branch addresses
	112	// **** neutrons
	113	Ntu->SetBranchAddress("Nleft",&fNnLeft);
	114	Ntu->SetBranchAddress("Eleft",&fEnLeft);
2524c56f	115	Ntu->SetBranchAddress("Pxl", fPxnLeft);
	116	Ntu->SetBranchAddress("Pyl", fPynLeft);
	117	Ntu->SetBranchAddress("Pzl", fPznLeft);
f1d1affa	118	Ntu->SetBranchAddress("Nright",&fNnRight);
f1d1affa	119	Ntu->SetBranchAddress("Eright",&fEnRight);
2524c56f	120	Ntu->SetBranchAddress("Pxr", fPxnRight);
	121	Ntu->SetBranchAddress("Pyr", fPynRight);
	122	Ntu->SetBranchAddress("Pzr", fPznRight);
f1d1affa	123	// **** protons
	124	Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
	125	Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
2524c56f	126	Ntu->SetBranchAddress("Pxl_p", fPxpLeft);
	127	Ntu->SetBranchAddress("Pyl_p", fPypLeft);
	128	Ntu->SetBranchAddress("Pzl_p", fPzpLeft);
f1d1affa	129	Ntu->SetBranchAddress("Nright_p",&fNpRight);
f1d1affa	130	Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
2524c56f	131	Ntu->SetBranchAddress("Pxr_p", fPxpRight);
	132	Ntu->SetBranchAddress("Pyr_p", fPypRight);
	133	Ntu->SetBranchAddress("Pzr_p", fPzpRight);
f1d1affa	134	}
	135
	136	// -----------------------------------------------------------------------------------
	137	Int_t AliGenReaderEMD::NextEvent()
	138	{
	139	// Read the next event
	140	Int_t nTracks=0;
	141
	142	TFile* pFile = fTreeNtuple->GetCurrentFile();
	143	pFile->cd();
	144
	145
	146	Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
	147	if(fNcurrent < nentries) {
	148	fTreeNtuple->GetEvent(fNcurrent);
	149	printf("\n * Event %d read *\n\n",fNcurrent);
	150	fNcurrent++;
	151	//printf("\n \t \t LEFT-> %d neutrons and %d protons emitted\n", fNnLeft, fNpLeft);
	152	//printf("\t \t RIGHT-> %d neutrons and %d protons emitted\n\n", fNnRight, fNpRight);
	153	//
	154	// **** IPSIde =0->RIGHT side, =1->LEFT side
	155	// #### fPcToTrack =0->neutrons, =1->protons
	156	if(fIPSide==0){
	157	if(fPcToTrack==0){
	158	printf("\n \t \t Tracking %d neutrons emitted on RIGHT side\n\n", fNnRight);
	159	nTracks = fNnRight;
	160	}
	161	else if(fPcToTrack==1){
	162	printf("\n \t \t Tracking %d protons emitted on RIGHT side\n\n", fNpRight);
	163	nTracks = fNpRight;
	164	}
	165	}
	166	else if(fIPSide==1){
	167	if(fPcToTrack==0){
	168	printf("\n \t \t Tracking %d neutrons emitted on LEFT side\n", fNnLeft);
	169	nTracks = fNnLeft;
	170	}
	171	else if(fPcToTrack==1){
	172	printf("\n \t \t Tracking %d protons emitted on LEFT side\n", fNpLeft);
	173	nTracks = fNpLeft;
	174	}
	175	}
	176	fNparticle = 0;
	177	return nTracks;
	178	}
	179
	180	return 0;
	181	}
	182
	183	// -----------------------------------------------------------------------------------
	184	TParticle* AliGenReaderEMD::NextParticle()
	185	{
	186	// Read the next particle
	187	Float_t p[4];
	188	Int_t ipart = kNeutron;
	189	Double_t amass = TDatabasePDG::Instance()->GetParticle(kNeutron)->Mass();
	190	p[0] = fPxnRight[fNparticle];
	191	p[1] = fPynRight[fNparticle];
	192	p[2] = fPznRight[fNparticle];
	193	Float_t ptot = TMath::Sqrt(p[0]p[0]+p[1]p[1]+p[2]*p[2]);
	194	p[3] = TMath::Sqrt(ptotptot+amassamass);
	195	//printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
	196
	197	if(p[3]<=amass)
198	Warning("Generate","Particle %d E = %f mass = %f \n",ipart,p[3],amass);
199
200	TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
201	p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
202	//fNcurrent++;
203	fNparticle++;
204	return particle;
205	}