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

// $Id$

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

#include <Riostream.h>
#include <TFile.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;
	    p->SetWeight(-123); //mark particle
	  }

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