[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"
#include "AliStack.h"


ClassImp(AliGenReaderEMD)

AliGenReaderEMD::AliGenReaderEMD():
    fStartEvent(0),
    fNcurrent(0),  
    fNparticle(0), 
    fTreeNtuple(0),
    fIPSide(0),
    fPcToTrack(0),
    fNnAside(0),
    fEnAside(0),
    fNnCside(0),
    fEnCside(0),
    fNpAside(0),
    fEtapAside(0),
    fNpCside(0),
    fEtapCside(0)
{
// Default constructor
}

AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
    AliGenReader(reader),
    fStartEvent(0),
    fNcurrent(0),  
    fNparticle(0), 
    fTreeNtuple(0),
    fIPSide(0),
    fPcToTrack(0),
    fNnAside(0),
    fEnAside(0),
    fNnCside(0),
    fEnCside(0),
    fNpAside(0),
    fEtapAside(0),
    fNpCside(0),
    fEtapCside(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",&fNnAside);
    Ntu->SetBranchAddress("Eleft",&fEnAside);
    Ntu->SetBranchAddress("Pxl",  fPxnAside);
    Ntu->SetBranchAddress("Pyl",  fPynAside);
    Ntu->SetBranchAddress("Pzl",  fPznAside);
    Ntu->SetBranchAddress("Nright",&fNnCside);
    Ntu->SetBranchAddress("Eright",&fEnCside);
    Ntu->SetBranchAddress("Pxr",   fPxnCside);
    Ntu->SetBranchAddress("Pyr",   fPynCside);
    Ntu->SetBranchAddress("Pzr",   fPznCside);
    // **** protons
    Ntu->SetBranchAddress("Nleft_p",&fNpAside);
    Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
    Ntu->SetBranchAddress("Pxl_p",  fPxpAside);
    Ntu->SetBranchAddress("Pyl_p",  fPypAside);
    Ntu->SetBranchAddress("Pzl_p",  fPzpAside);
    Ntu->SetBranchAddress("Nright_p",&fNpCside);
    Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
    Ntu->SetBranchAddress("Pxr_p",   fPxpCside);
    Ntu->SetBranchAddress("Pyr_p",   fPypCside);
    Ntu->SetBranchAddress("Pzr_p",   fPzpCside);
}

// -----------------------------------------------------------------------------------
Int_t AliGenReaderEMD::NextEvent() 
{
    // Read the next event  
    Int_t nTracks=0;
    fNparticle = 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",fNcurrent);
	//printf("\t A side-> %d neutrons and %d protons emitted\n", fNnAside, fNpAside);
	//printf("\t C side-> %d neutrons and %d protons emitted\n\n", fNnCside, fNpCside);
	//
	// **** IPSIde = 0->both sides, 1->only Cside, 2->only A side  
	// #### fPcToTrack = 0->neutrons, 1->protons
	if(fPcToTrack==0){ // neutrons
	  if(fIPSide==0){
	    printf("\t Tracking %d+%d neutrons emitted on C+A side\n", fNnCside,fNnAside);
	    nTracks = fNnCside+fNnAside;
	  }
	  else if(fIPSide==1){ 
	    printf("\t Tracking %d neutrons emitted on C side\n", fNnCside);
	    nTracks    = fNnCside;
	  }
	  else if(fIPSide==2){
	    printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
	    nTracks    = fNnAside;
	  }
	}
	else if(fPcToTrack==1){ //protons
	  if(fIPSide==0){
	    printf("\t Tracking %d+%d protons emitted on C+A side\n", fNpCside,fNpAside);
	    nTracks    = fNpCside+fNpAside;
	  }
	  else if(fIPSide==1){
	    printf("\t Tracking %d protons emitted on C side\n", fNpCside);
	    nTracks    = fNpCside;
	  }
	  else if(fIPSide==2){
	    printf("\t Tracking %d protons emitted on A side\n", fNpAside);
	    nTracks    = fNpAside;
	  }
	}
	fNcurrent++;
	return nTracks;
    }

    return 0;
}

// -----------------------------------------------------------------------------------
TParticle* AliGenReaderEMD::NextParticle() 
{
    // Read the next particle
    Float_t p[4]={0.,0.,0.,0.};
    
    Int_t ipart=0;
    if(fPcToTrack==0) ipart = kNeutron;
    else  if(fPcToTrack==1) ipart = kProton;
    Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();

    if(fPcToTrack==0){
      if(fIPSide==0){
        if(fNparticle<fNnCside){
          p[0] = fPxnCside[fNparticle];
          p[1] = fPynCside[fNparticle];
          p[2] = fPznCside[fNparticle];
	}
	else{
          p[0] = fPxnAside[fNparticle-fNnCside];
          p[1] = fPynAside[fNparticle-fNnCside];
          p[2] = fPznAside[fNparticle-fNnCside];
	}
      }
      else if(fIPSide==1){
        p[0] = fPxnCside[fNparticle];
        p[1] = fPynCside[fNparticle];
        p[2] = fPznCside[fNparticle];
      }
      else if(fIPSide==2){
        p[0] = fPxnAside[fNparticle-fNpCside];
        p[1] = fPynAside[fNparticle-fNpCside];
        p[2] = fPznAside[fNparticle-fNpCside];
      }
    }
    else if(fPcToTrack==1){
      if(fIPSide==0){
        if(fNparticle<fNnCside){
          p[0] = fPxpCside[fNparticle];
          p[1] = fPypCside[fNparticle];
          p[2] = fPzpCside[fNparticle];
	}
	else{
          p[0] = fPxpAside[fNparticle];
          p[1] = fPypAside[fNparticle];
          p[2] = fPzpAside[fNparticle];
	}
      }
      else if(fIPSide==1){
        p[0] = fPxpCside[fNparticle];
        p[1] = fPypCside[fNparticle];
        p[2] = fPzpCside[fNparticle];
      }
      else if(fIPSide==2){
        p[0] = fPxpAside[fNparticle];
        p[1] = fPypAside[fNparticle];
        p[2] = fPzpAside[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(" * pc %d: momentum (%f, %f, %f) \n", fNparticle, 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.);
    particle->SetBit(kTransportBit);
    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>
f1d1affa	27	#include "AliGenReaderEMD.h"
2078b0f7	28	#include "AliStack.h"
2078b0f7	29
f1d1affa	30
706938e6	31	ClassImp(AliGenReaderEMD)
f1d1affa	32
1c56e311	33	AliGenReaderEMD::AliGenReaderEMD():
	34	fStartEvent(0),
	35	fNcurrent(0),
	36	fNparticle(0),
	37	fTreeNtuple(0),
	38	fIPSide(0),
	39	fPcToTrack(0),
2078b0f7	40	fNnAside(0),
	41	fEnAside(0),
	42	fNnCside(0),
	43	fEnCside(0),
	44	fNpAside(0),
	45	fEtapAside(0),
	46	fNpCside(0),
	47	fEtapCside(0)
f1d1affa	48	{
f1d1affa	49	// Default constructor
f1d1affa	50	}
f1d1affa	51
1c56e311	52	AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
	53	AliGenReader(reader),
	54	fStartEvent(0),
	55	fNcurrent(0),
	56	fNparticle(0),
	57	fTreeNtuple(0),
	58	fIPSide(0),
	59	fPcToTrack(0),
2078b0f7	60	fNnAside(0),
	61	fEnAside(0),
	62	fNnCside(0),
	63	fEnCside(0),
	64	fNpAside(0),
	65	fEtapAside(0),
	66	fNpCside(0),
	67	fEtapCside(0)
1c56e311	68	{
	69	// Copy Constructor
	70	reader.Copy(*this);
	71	}
f1d1affa	72	// -----------------------------------------------------------------------------------
	73	AliGenReaderEMD::~AliGenReaderEMD()
	74	{
	75	delete fTreeNtuple;
	76	}
	77
	78	// -----------------------------------------------------------------------------------
	79	AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
	80	{
	81	// Assignment operator
	82	rhs.Copy(*this);
	83	return *this;
	84	}
	85
	86	// -----------------------------------------------------------------------------------
	87	void AliGenReaderEMD::Copy(TObject&) const
	88	{
	89	//
	90	// Copy
	91	//
	92	Fatal("Copy","Not implemented!\n");
	93	}
	94
	95	// -----------------------------------------------------------------------------------
	96	void AliGenReaderEMD::Init()
	97	{
	98	//
	99	// Reset the existing file environment and open a new root file
	100
	101	TFile *pFile=0;
	102	if (!pFile) {
	103	pFile = new TFile(fFileName);
	104	pFile->cd();
	105	printf("\n %s file opened\n\n", fFileName);
	106	}
	107	fTreeNtuple = (TTree*)gDirectory->Get("h2032");
	108	fNcurrent = fStartEvent;
	109
	110	TTree *Ntu=fTreeNtuple;
	111	//
	112	// Set branch addresses
	113	// **** neutrons
2078b0f7	114	Ntu->SetBranchAddress("Nleft",&fNnAside);
	115	Ntu->SetBranchAddress("Eleft",&fEnAside);
	116	Ntu->SetBranchAddress("Pxl", fPxnAside);
	117	Ntu->SetBranchAddress("Pyl", fPynAside);
	118	Ntu->SetBranchAddress("Pzl", fPznAside);
	119	Ntu->SetBranchAddress("Nright",&fNnCside);
	120	Ntu->SetBranchAddress("Eright",&fEnCside);
	121	Ntu->SetBranchAddress("Pxr", fPxnCside);
	122	Ntu->SetBranchAddress("Pyr", fPynCside);
	123	Ntu->SetBranchAddress("Pzr", fPznCside);
f1d1affa	124	// **** protons
2078b0f7	125	Ntu->SetBranchAddress("Nleft_p",&fNpAside);
	126	Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
	127	Ntu->SetBranchAddress("Pxl_p", fPxpAside);
	128	Ntu->SetBranchAddress("Pyl_p", fPypAside);
	129	Ntu->SetBranchAddress("Pzl_p", fPzpAside);
	130	Ntu->SetBranchAddress("Nright_p",&fNpCside);
	131	Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
	132	Ntu->SetBranchAddress("Pxr_p", fPxpCside);
	133	Ntu->SetBranchAddress("Pyr_p", fPypCside);
	134	Ntu->SetBranchAddress("Pzr_p", fPzpCside);
f1d1affa	135	}
	136
	137	// -----------------------------------------------------------------------------------
	138	Int_t AliGenReaderEMD::NextEvent()
	139	{
	140	// Read the next event
	141	Int_t nTracks=0;
2078b0f7	142	fNparticle = 0;
f1d1affa	143
	144	TFile* pFile = fTreeNtuple->GetCurrentFile();
	145	pFile->cd();
	146
	147
	148	Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
	149	if(fNcurrent < nentries) {
	150	fTreeNtuple->GetEvent(fNcurrent);
2078b0f7	151	printf("\n * Event %d read *\n",fNcurrent);
	152	//printf("\t A side-> %d neutrons and %d protons emitted\n", fNnAside, fNpAside);
	153	//printf("\t C side-> %d neutrons and %d protons emitted\n\n", fNnCside, fNpCside);
f1d1affa	154	//
2078b0f7	155	// **** IPSIde = 0->both sides, 1->only Cside, 2->only A side
	156	// #### fPcToTrack = 0->neutrons, 1->protons
	157	if(fPcToTrack==0){ // neutrons
	158	if(fIPSide==0){
	159	printf("\t Tracking %d+%d neutrons emitted on C+A side\n", fNnCside,fNnAside);
	160	nTracks = fNnCside+fNnAside;
	161	}
	162	else if(fIPSide==1){
	163	printf("\t Tracking %d neutrons emitted on C side\n", fNnCside);
	164	nTracks = fNnCside;
f1d1affa	165	}
2078b0f7	166	else if(fIPSide==2){
	167	printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
	168	nTracks = fNnAside;
f1d1affa	169	}
f1d1affa	170	}
2078b0f7	171	else if(fPcToTrack==1){ //protons
	172	if(fIPSide==0){
	173	printf("\t Tracking %d+%d protons emitted on C+A side\n", fNpCside,fNpAside);
	174	nTracks = fNpCside+fNpAside;
f1d1affa	175	}
2078b0f7	176	else if(fIPSide==1){
	177	printf("\t Tracking %d protons emitted on C side\n", fNpCside);
	178	nTracks = fNpCside;
	179	}
	180	else if(fIPSide==2){
	181	printf("\t Tracking %d protons emitted on A side\n", fNpAside);
	182	nTracks = fNpAside;
f1d1affa	183	}
f1d1affa	184	}
2078b0f7	185	fNcurrent++;
f1d1affa	186	return nTracks;
	187	}
	188
	189	return 0;
	190	}
	191
	192	// -----------------------------------------------------------------------------------
	193	TParticle* AliGenReaderEMD::NextParticle()
	194	{
	195	// Read the next particle
2078b0f7	196	Float_t p[4]={0.,0.,0.,0.};
2078b0f7	197
ef064345	198	Int_t ipart=0;
	199	if(fPcToTrack==0) ipart = kNeutron;
	200	else if(fPcToTrack==1) ipart = kProton;
	201	Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
	202
2078b0f7	203	if(fPcToTrack==0){
	204	if(fIPSide==0){
	205	if(fNparticle<fNnCside){
	206	p[0] = fPxnCside[fNparticle];
	207	p[1] = fPynCside[fNparticle];
	208	p[2] = fPznCside[fNparticle];
	209	}
	210	else{
	211	p[0] = fPxnAside[fNparticle-fNnCside];
	212	p[1] = fPynAside[fNparticle-fNnCside];
	213	p[2] = fPznAside[fNparticle-fNnCside];
	214	}
	215	}
	216	else if(fIPSide==1){
	217	p[0] = fPxnCside[fNparticle];
	218	p[1] = fPynCside[fNparticle];
	219	p[2] = fPznCside[fNparticle];
	220	}
	221	else if(fIPSide==2){
	222	p[0] = fPxnAside[fNparticle-fNpCside];
	223	p[1] = fPynAside[fNparticle-fNpCside];
	224	p[2] = fPznAside[fNparticle-fNpCside];
	225	}
	226	}
	227	else if(fPcToTrack==1){
	228	if(fIPSide==0){
	229	if(fNparticle<fNnCside){
	230	p[0] = fPxpCside[fNparticle];
	231	p[1] = fPypCside[fNparticle];
	232	p[2] = fPzpCside[fNparticle];
	233	}
	234	else{
	235	p[0] = fPxpAside[fNparticle];
	236	p[1] = fPypAside[fNparticle];
	237	p[2] = fPzpAside[fNparticle];
	238	}
	239	}
	240	else if(fIPSide==1){
	241	p[0] = fPxpCside[fNparticle];
	242	p[1] = fPypCside[fNparticle];
	243	p[2] = fPzpCside[fNparticle];
	244	}
	245	else if(fIPSide==2){
	246	p[0] = fPxpAside[fNparticle];
	247	p[1] = fPypAside[fNparticle];
	248	p[2] = fPzpAside[fNparticle];
	249	}
	250	}
f1d1affa	251	Float_t ptot = TMath::Sqrt(p[0]p[0]+p[1]p[1]+p[2]*p[2]);
f1d1affa	252	p[3] = TMath::Sqrt(ptotptot+amassamass);
2078b0f7	253	//printf(" * pc %d: momentum (%f, %f, %f) \n", fNparticle, p[0],p[1],p[2]);
f1d1affa	254
	255	if(p[3]<=amass)
	256	Warning("Generate","Particle %d E = %f mass = %f \n",ipart,p[3],amass);
	257
	258	TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
	259	p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
2078b0f7	260	particle->SetBit(kTransportBit);
f1d1affa	261	fNparticle++;
	262	return particle;
	263	}