[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderKine.cxx

// $Id$

//_______________________________________________________________________
/////////////////////////////////////////////////////////////////////////
//
// class AliJetParticlesReaderKine
//
// Reader for Kinematics
//
// loizides@ikf.uni-frankfurt.de
//
/////////////////////////////////////////////////////////////////////////

#include <Riostream.h>
#include <TString.h>
#include <TParticle.h>
#include <TLorentzVector.h>
#include <AliRunLoader.h>
#include <AliStack.h>
#include <AliHeader.h>
#include <AliGenEventHeader.h>
#include <AliGenPythiaEventHeader.h>
#include <AliGenHijingEventHeader.h>
#include "AliJetParticle.h"
#include "AliJetEventParticles.h"
#include "AliJetParticlesReaderKine.h"

ClassImp(AliJetParticlesReaderKine)

AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
  AliJetParticlesReader(),
  fFileName("galice.root"),
  fRunLoader(0),
  fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  fUseTracks(kFALSE)
{
  //constructor
}

AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
  AliJetParticlesReader(),
  fFileName(fname),
  fRunLoader(0),
  fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  fUseTracks(kFALSE)
{
  //constructor
}

AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
  AliJetParticlesReader(dirs),
  fFileName(filename),
  fRunLoader(0),
  fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  fUseTracks(kFALSE)
{
  //constructor
}

AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
{
  //destructor
  if(fRunLoader) delete fRunLoader;
}

void AliJetParticlesReaderKine::Rewind()
{
  //Rewinds to the beginning
  if(fRunLoader) delete fRunLoader;
  fRunLoader   = 0;
  fCurrentDir  = 0;
  fNEventsRead = 0;  
}

Int_t AliJetParticlesReaderKine::ReadNext()
{
  //Reads Kinematics Tree

  if((!fOwner) || (fEventParticles == 0)) 
    fEventParticles = new AliJetEventParticles();

  while(fCurrentDir < GetNumberOfDirs())
    { 

      if (!OpenFile(fCurrentDir)) 
      { 
	delete fRunLoader; //close current session
	fRunLoader = 0;    //assure pointer is null
	fCurrentDir++;
	continue;
      }
    
      if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
	{
	  //read next directory
	  delete fRunLoader; //close current session
	  fRunLoader = 0;    //assure pointer is null
	  fCurrentDir++;     //go to next dir
	  continue; 
	}
     
      Info("ReadNext","Reading Event %d",fCurrentEvent);
      fRunLoader->GetEvent(fCurrentEvent);
      AliStack* stack = fRunLoader->Stack();
      if (!stack)
	{
	  Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
	  continue;
	}

      //clear old event
      Int_t nprim = stack->GetNprimary();
      Int_t npart = nprim;
      if(fUseTracks)
	npart = stack->GetNtrack();
      fEventParticles->Reset(npart);

      TString headdesc="";
      Int_t gentype=-1;
      AliHeader *header=fRunLoader->GetHeader();
      if(!header) {
	Warning("ReadNext","Header not found in event %d",fCurrentEvent);
      } else {
	AliGenEventHeader *genheader=header->GenEventHeader();
	TString hname=genheader->GetName();
	if(hname.CompareTo("HIJINGparam")==0) {
	  gentype=1;	  ;
	  headdesc+="HIJINGparam";
	} else if(hname.CompareTo("Hijing")==0) {
	  gentype=2;
	  AliGenHijingEventHeader *hheader=(AliGenHijingEventHeader*)header->GenEventHeader();
	  headdesc+="Hijing";
	  if(!hheader) {
	    Warning("ReadNext","Hijing-Header not found in event %d",fCurrentEvent);
	  } else {
	    Int_t ntrials=hheader->Trials();
	    fEventParticles->SetTrials(ntrials);
	    fEventParticles->SetImpact(hheader->ImpactParameter());
	    fEventParticles->SetNhard(hheader->HardScatters());
	    fEventParticles->SetNpart(hheader->NwNw());

	    TLorentzVector  jet1; 
	    TLorentzVector  jet2; 
	    TLorentzVector  jetFsr1;
	    TLorentzVector  jetFsr2;
	    hheader->GetJets(jet1,jet2,jetFsr1,jetFsr2);
	    fEventParticles->AddHard(0,jet1.X(),jet1.Y(),jet1.Z(),jet1.T(),0);
	    fEventParticles->AddHard(1,jet2.X(),jet2.Y(),jet2.Z(),jet2.T(),0);
	  }
	} else if(hname.CompareTo("Pythia")==0){
	  gentype=3;
	  headdesc+="Pythia";
	  AliGenPythiaEventHeader *pheader=(AliGenPythiaEventHeader*)header->GenEventHeader();
	  if(!pheader) {
	    Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
	  } else {
	    Int_t ntruq=0;
#ifndef NOUQHEADERINFO
	    ntruq=pheader->NUQTriggerJets();
	    if(ntruq){
	      Double_t x0=pheader->GetXJet();
	      Double_t y0=pheader->GetYJet();
	      if(x0==y0==-1){
		x0=y0=0.;
	      }
	      Double_t zquench[4];
	      pheader->GetZQuench(zquench);
	      fEventParticles->SetXYJet(x0,y0);
	      fEventParticles->SetZQuench(zquench);
	      for(Int_t j=0;j<ntruq;j++){
		Float_t pjet[4];
		pheader->UQJet(j,pjet);
		fEventParticles->AddUQJet(pjet);
	      }
	    }
#endif
	    //Int_t ptyp=pheader->ProcessType();
	    Int_t ntrials=pheader->Trials();
	    fEventParticles->SetTrials(ntrials);

	    Int_t ntr=pheader->NTriggerJets();
	    if(ntr){
	      for(Int_t j=0;j<ntr;j++){
		Float_t pjet[4];
		pheader->TriggerJet(j,pjet);
		fEventParticles->AddJet(pjet);
		if(!ntruq) fEventParticles->AddUQJet(pjet);
	      }
	    }
	    for(Int_t i=6;i<=7;i++){
	      TParticle *MP = stack->Particle(i);
	      if(!MP) break;
	      Int_t type=0;
	      if(MP->GetPdgCode()==21) type=2;
	      else type=1;
	      fEventParticles->AddHard(i-6,MP->Px(),MP->Py(),MP->Pz(),MP->Energy(),type);
	    }
	  }
	}
	headdesc+="Run ";
	headdesc+=header->GetRun();
	headdesc+=": Ev ";
	headdesc+=header->GetEventNrInRun();
      }
      fEventParticles->SetHeader(headdesc);

      //get vertex
      const TParticle *kv = stack->Particle(0);
      if(kv) {
	fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
      }

      //loop over particles
      for (Int_t i = 0;i<npart; i++)
	{
	  TParticle *p = stack->Particle(i);
	  if(!p) continue;
	  Int_t child1 = p->GetFirstDaughter();
	  //Int_t child2 = p->GetLastDaughter();	
	  //Int_t mother = p->GetFirstMother();	   
	  //cout << child1 << " " << child2 << " " << mother << endl;
	  if((child1>=0) && (child1<nprim)) continue; 

	  //check status code depending on gentype
	  if(gentype==1){ // Hijing Param
	    if(p->GetStatusCode()!=0) continue;
	  } else if(gentype==2){ //Hijing
	    //if( p->GetStatusCode()!=0){
	    //  cout <<p->GetStatusCode() << " ";p->Print(); 
	    //}
	  } else if(gentype==3){ //Pythia
	    if(p->GetStatusCode()!=1) continue;
	  }

	  //kinematic cuts
	  if(IsAcceptedParticle(p)){ //put particle in event
	    //if(p->Pt()>20){cout <<p->GetStatusCode() << " ";p->Print(); }
	    fEventParticles->AddParticle(p,i); 
	  }
	}
      fCurrentEvent++;
      fNEventsRead++;
      return kTRUE;
    }

  //end of loop over directories specified in fDirs ObjArray
  return kFALSE;
}

Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
{
  //opens file with kine tree

  if(fRunLoader){
    if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
    else return kFALSE;
  }

  const TString& dirname = GetDirName(n);
  if (dirname == "")
    { 
      Error("OpenNextFile","Can't get directory name with index %d",n);
      return kFALSE;
    }

  TString filename = dirname +"/"+ fFileName;
  fRunLoader = AliRunLoader::Open(filename.Data()); 

  if (fRunLoader == 0)
    {
      Error("OpenNextFile","Can't open session from file %s",filename.Data());
      return kFALSE;
    }
  
  if (fRunLoader->GetNumberOfEvents() <= 0)
    {
      Error("OpenNextFile","There is no events in this directory.");
      delete fRunLoader;
      fRunLoader = 0;
      return kFALSE;
    }

  if (fRunLoader->LoadKinematics())
    {
      Error("OpenNextFile","Error occured while loading kinematics.");
      return kFALSE;
    }

  fCurrentEvent = 0;
  return kTRUE;
}
Commit	Line	Data
d7c6ab14	1	// $Id$
	2
	3	//_______________________________________________________________________
	4	/////////////////////////////////////////////////////////////////////////
	5	//
	6	// class AliJetParticlesReaderKine
	7	//
	8	// Reader for Kinematics
	9	//
	10	// loizides@ikf.uni-frankfurt.de
	11	//
	12	/////////////////////////////////////////////////////////////////////////
	13
	14	#include <Riostream.h>
	15	#include <TString.h>
	16	#include <TParticle.h>
d7073a83	17	#include <TLorentzVector.h>
d7c6ab14	18	#include <AliRunLoader.h>
d7c6ab14	19	#include <AliStack.h>
d7073a83	20	#include <AliHeader.h>
	21	#include <AliGenEventHeader.h>
	22	#include <AliGenPythiaEventHeader.h>
	23	#include <AliGenHijingEventHeader.h>
d7c6ab14	24	#include "AliJetParticle.h"
	25	#include "AliJetEventParticles.h"
	26	#include "AliJetParticlesReaderKine.h"
	27
	28	ClassImp(AliJetParticlesReaderKine)
	29
d7c6ab14	30	AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
	31	AliJetParticlesReader(),
	32	fFileName("galice.root"),
	33	fRunLoader(0),
	34	fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
	35	fUseTracks(kFALSE)
	36	{
	37	//constructor
	38	}
	39
d7c6ab14	40	AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
	41	AliJetParticlesReader(),
	42	fFileName(fname),
	43	fRunLoader(0),
	44	fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
	45	fUseTracks(kFALSE)
	46	{
	47	//constructor
	48	}
	49
d7c6ab14	50	AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
	51	AliJetParticlesReader(dirs),
	52	fFileName(filename),
	53	fRunLoader(0),
	54	fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
	55	fUseTracks(kFALSE)
	56	{
	57	//constructor
	58	}
	59
d7c6ab14	60	AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
	61	{
	62	//destructor
	63	if(fRunLoader) delete fRunLoader;
	64	}
	65
d7c6ab14	66	void AliJetParticlesReaderKine::Rewind()
	67	{
	68	//Rewinds to the beginning
	69	if(fRunLoader) delete fRunLoader;
04a02430	70	fRunLoader = 0;
	71	fCurrentDir = 0;
	72	fNEventsRead = 0;
d7c6ab14	73	}
d7c6ab14	74
d7c6ab14	75	Int_t AliJetParticlesReaderKine::ReadNext()
	76	{
	77	//Reads Kinematics Tree
	78
	79	if((!fOwner) \|\| (fEventParticles == 0))
	80	fEventParticles = new AliJetEventParticles();
	81
b2760c9e	82	while(fCurrentDir < GetNumberOfDirs())
d7c6ab14	83	{
04a02430	84
04a02430	85	if (!OpenFile(fCurrentDir))
d7c6ab14	86	{
04a02430	87	delete fRunLoader; //close current session
04a02430	88	fRunLoader = 0; //assure pointer is null
d7c6ab14	89	fCurrentDir++;
	90	continue;
	91	}
	92
	93	if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
	94	{
	95	//read next directory
	96	delete fRunLoader; //close current session
	97	fRunLoader = 0; //assure pointer is null
	98	fCurrentDir++; //go to next dir
	99	continue;
	100	}
	101
04a02430	102	Info("ReadNext","Reading Event %d",fCurrentEvent);
d7c6ab14	103	fRunLoader->GetEvent(fCurrentEvent);
	104	AliStack* stack = fRunLoader->Stack();
	105	if (!stack)
	106	{
04a02430	107	Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
d7c6ab14	108	continue;
	109	}
	110
04a02430	111	//clear old event
	112	Int_t nprim = stack->GetNprimary();
	113	Int_t npart = nprim;
	114	if(fUseTracks)
	115	npart = stack->GetNtrack();
	116	fEventParticles->Reset(npart);
	117
	118	TString headdesc="";
d7073a83	119	Int_t gentype=-1;
04a02430	120	AliHeader *header=fRunLoader->GetHeader();
	121	if(!header) {
	122	Warning("ReadNext","Header not found in event %d",fCurrentEvent);
	123	} else {
d7073a83	124	AliGenEventHeader *genheader=header->GenEventHeader();
	125	TString hname=genheader->GetName();
	126	if(hname.CompareTo("HIJINGparam")==0) {
	127	gentype=1; ;
	128	headdesc+="HIJINGparam";
	129	} else if(hname.CompareTo("Hijing")==0) {
	130	gentype=2;
	131	AliGenHijingEventHeader hheader=(AliGenHijingEventHeader)header->GenEventHeader();
	132	headdesc+="Hijing";
	133	if(!hheader) {
	134	Warning("ReadNext","Hijing-Header not found in event %d",fCurrentEvent);
	135	} else {
	136	Int_t ntrials=hheader->Trials();
	137	fEventParticles->SetTrials(ntrials);
	138	fEventParticles->SetImpact(hheader->ImpactParameter());
	139	fEventParticles->SetNhard(hheader->HardScatters());
	140	fEventParticles->SetNpart(hheader->NwNw());
	141
	142	TLorentzVector jet1;
	143	TLorentzVector jet2;
	144	TLorentzVector jetFsr1;
	145	TLorentzVector jetFsr2;
	146	hheader->GetJets(jet1,jet2,jetFsr1,jetFsr2);
	147	fEventParticles->AddHard(0,jet1.X(),jet1.Y(),jet1.Z(),jet1.T(),0);
	148	fEventParticles->AddHard(1,jet2.X(),jet2.Y(),jet2.Z(),jet2.T(),0);
	149	}
	150	} else if(hname.CompareTo("Pythia")==0){
	151	gentype=3;
	152	headdesc+="Pythia";
	153	AliGenPythiaEventHeader pheader=(AliGenPythiaEventHeader)header->GenEventHeader();
	154	if(!pheader) {
	155	Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
	156	} else {
	157	Int_t ntruq=0;
042d3031	158	#ifndef NOUQHEADERINFO
d7073a83	159	ntruq=pheader->NUQTriggerJets();
	160	if(ntruq){
	161	Double_t x0=pheader->GetXJet();
	162	Double_t y0=pheader->GetYJet();
	163	if(x0==y0==-1){
	164	x0=y0=0.;
	165	}
	166	Double_t zquench[4];
	167	pheader->GetZQuench(zquench);
	168	fEventParticles->SetXYJet(x0,y0);
	169	fEventParticles->SetZQuench(zquench);
	170	for(Int_t j=0;j<ntruq;j++){
	171	Float_t pjet[4];
	172	pheader->UQJet(j,pjet);
	173	fEventParticles->AddUQJet(pjet);
	174	}
04a02430	175	}
042d3031	176	#endif
d7073a83	177	//Int_t ptyp=pheader->ProcessType();
	178	Int_t ntrials=pheader->Trials();
	179	fEventParticles->SetTrials(ntrials);
	180
	181	Int_t ntr=pheader->NTriggerJets();
	182	if(ntr){
	183	for(Int_t j=0;j<ntr;j++){
	184	Float_t pjet[4];
	185	pheader->TriggerJet(j,pjet);
	186	fEventParticles->AddJet(pjet);
	187	if(!ntruq) fEventParticles->AddUQJet(pjet);
	188	}
	189	}
	190	for(Int_t i=6;i<=7;i++){
	191	TParticle *MP = stack->Particle(i);
	192	if(!MP) break;
	193	Int_t type=0;
	194	if(MP->GetPdgCode()==21) type=2;
	195	else type=1;
	196	fEventParticles->AddHard(i-6,MP->Px(),MP->Py(),MP->Pz(),MP->Energy(),type);
04a02430	197	}
	198	}
	199	}
d7073a83	200	headdesc+="Run ";
	201	headdesc+=header->GetRun();
	202	headdesc+=": Ev ";
	203	headdesc+=header->GetEventNrInRun();
04a02430	204	}
	205	fEventParticles->SetHeader(headdesc);
	206
d7c6ab14	207	//get vertex
	208	const TParticle *kv = stack->Particle(0);
	209	if(kv) {
	210	fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
	211	}
	212
04a02430	213	//loop over particles
d7c6ab14	214	for (Int_t i = 0;i<npart; i++)
	215	{
	216	TParticle *p = stack->Particle(i);
	217	if(!p) continue;
b2760c9e	218	Int_t child1 = p->GetFirstDaughter();
	219	//Int_t child2 = p->GetLastDaughter();
	220	//Int_t mother = p->GetFirstMother();
b2760c9e	221	//cout << child1 << " " << child2 << " " << mother << endl;
04a02430	222	if((child1>=0) && (child1<nprim)) continue;
d7073a83	223
	224	//check status code depending on gentype
	225	if(gentype==1){ // Hijing Param
	226	if(p->GetStatusCode()!=0) continue;
	227	} else if(gentype==2){ //Hijing
	228	//if( p->GetStatusCode()!=0){
	229	// cout <<p->GetStatusCode() << " ";p->Print();
	230	//}
	231	} else if(gentype==3){ //Pythia
	232	if(p->GetStatusCode()!=1) continue;
04a02430	233	}
d7073a83	234
	235	//kinematic cuts
	236	if(IsAcceptedParticle(p)){ //put particle in event
	237	//if(p->Pt()>20){cout <<p->GetStatusCode() << " ";p->Print(); }
d7c6ab14	238	fEventParticles->AddParticle(p,i);
d7073a83	239	}
d7c6ab14	240	}
d7c6ab14	241	fCurrentEvent++;
	242	fNEventsRead++;
	243	return kTRUE;
b2760c9e	244	}
04a02430	245
04a02430	246	//end of loop over directories specified in fDirs ObjArray
d7c6ab14	247	return kFALSE;
	248	}
	249
d7c6ab14	250	Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
	251	{
	252	//opens file with kine tree
	253
04a02430	254	if(fRunLoader){
	255	if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
	256	else return kFALSE;
	257	}
	258
d7c6ab14	259	const TString& dirname = GetDirName(n);
	260	if (dirname == "")
	261	{
04a02430	262	Error("OpenNextFile","Can't get directory name with index %d",n);
d7c6ab14	263	return kFALSE;
	264	}
	265
	266	TString filename = dirname +"/"+ fFileName;
	267	fRunLoader = AliRunLoader::Open(filename.Data());
	268
	269	if (fRunLoader == 0)
	270	{
	271	Error("OpenNextFile","Can't open session from file %s",filename.Data());
	272	return kFALSE;
	273	}
	274
	275	if (fRunLoader->GetNumberOfEvents() <= 0)
	276	{
	277	Error("OpenNextFile","There is no events in this directory.");
	278	delete fRunLoader;
	279	fRunLoader = 0;
	280	return kFALSE;
	281	}
04a02430	282
d7c6ab14	283	if (fRunLoader->LoadKinematics())
	284	{
	285	Error("OpenNextFile","Error occured while loading kinematics.");
	286	return kFALSE;
	287	}
04a02430	288
d7c6ab14	289	fCurrentEvent = 0;
	290	return kTRUE;
	291	}