[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESD.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoEventReaderESD - the reader class for the Alice ESD              ///
/// Reads in ESD information and converts it into internal AliFemtoEvent     ///
/// Reads in AliESDfriend to create shared hit/quality information           ///
/// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl                        ///
///          Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////

/*
 *$Id$
 *$Log$
 *Revision 1.2.2.2  2007/10/04 13:10:52  akisiel
 *Add Kink index storageAliFemtoEventReaderESD.cxx AliFemtoTrack.cxx AliFemtoTrack.h
 *
 *Revision 1.2.2.1  2007/09/30 11:38:59  akisiel
 *Adapt the readers to the new AliESDEvent structure
 *
 *Revision 1.2  2007/05/22 09:01:42  akisiel
 *Add the possibiloity to save cut settings in the ROOT file
 *
 *Revision 1.1  2007/05/16 10:22:11  akisiel
 *Making the directory structure of AliFemto flat. All files go into one common directory
 *
 *Revision 1.5  2007/05/03 09:45:20  akisiel
 *Fixing Effective C++ warnings
 *
 *Revision 1.4  2007/04/27 07:28:34  akisiel
 *Remove event number reading due to interface changes
 *
 *Revision 1.3  2007/04/27 07:25:16  akisiel
 *Make revisions needed for compilation from the main AliRoot tree
 *
 *Revision 1.1.1.1  2007/04/25 15:38:41  panos
 *Importing the HBT code dir
 *
 */

#include "AliFemtoEventReaderESD.h"

#include "TFile.h"
#include "TTree.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"

//#include "TSystem.h"

#include "AliFmPhysicalHelixD.h"
#include "AliFmThreeVectorF.h"

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"
#include "AliFemtoModelHiddenInfo.h"

ClassImp(AliFemtoEventReaderESD)

#if !(ST_NO_NAMESPACES)
  using namespace units;
#endif

using namespace std;
//____________________________
//constructor with 0 parameters , look at default settings 
AliFemtoEventReaderESD::AliFemtoEventReaderESD():
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fReadInner(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fTree(0x0),
  fEsdFile(0x0),
  fEvent(0x0)
{
  // default constructor
}

AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
  AliFemtoEventReader(aReader),
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fReadInner(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fTree(0x0),
  fEsdFile(0x0),
  fEvent(0x0)
{
  // copy constructor
  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fReadInner = aReader.fReadInner;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  //  fTree = aReader.fTree->CloneTree();
  //  fEvent = new AliESD(*aReader.fEvent);
  fEvent = new AliESDEvent();
  fEsdFile = new TFile(aReader.fEsdFile->GetName());
}
//__________________
//Destructor
AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
{
  // destructor
  //delete fListOfFiles;
  delete fTree;
  delete fEvent;
  delete fEsdFile;
}

//__________________
AliFemtoEventReaderESD& AliFemtoEventReaderESD::operator=(const AliFemtoEventReaderESD& aReader)
{
  // assignment operator
  if (this == &aReader)
    return *this;

  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fReadInner = aReader.fReadInner;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  if (fTree) delete fTree;
  //  fTree = aReader.fTree->CloneTree();
  if (fEvent) delete fEvent;
  fEvent = new AliESDEvent();
  if (fEsdFile) delete fEsdFile;
  fEsdFile = new TFile(aReader.fEsdFile->GetName());

  return *this;
}
//__________________
AliFemtoString AliFemtoEventReaderESD::Report()
{
  // create reader report
  AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
  return temp;
}

//__________________
void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
{
  //setting the name of file where names of ESD file are written 
  //it takes only this files which have good trees
  char buffer[256];
  fInputFile=string(inputFile);
  cout<<"Input File set on "<<fInputFile<<endl;
  ifstream infile(inputFile);

  fTree = new TChain("esdTree");

  if(infile.good()==true)
    { 
      //checking if all give files have good tree inside
      while (infile.eof()==false)
	{
	  infile.getline(buffer,256);
	  //ifstream test_file(buffer);
	  TFile *esdFile=TFile::Open(buffer,"READ");
	  if (esdFile!=0x0)
	    {	
	      TTree* tree = (TTree*) esdFile->Get("esdTree");
	      if (tree!=0x0)
		{
		  cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
		  fTree->AddFile(buffer);
		  delete tree;
		}
	      esdFile->Close();	
	    }
	  delete esdFile;
	}
    }
}

void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}

bool AliFemtoEventReaderESD::GetConstrained() const
{
  return fConstrained;
}

void AliFemtoEventReaderESD::SetReadTPCInner(const bool readinner)
{
  fReadInner=readinner;
}

bool AliFemtoEventReaderESD::GetReadTPCInner() const
{
  return fReadInner;
}

AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
{
  // read in a next hbt event from the chain
  // convert it to AliFemtoEvent and return
  // for further analysis
  AliFemtoEvent *hbtEvent = 0;

  if (fCurEvent==fNumberofEvent)//open next file  
    {
      if(fNumberofEvent==0)	
	{
	  //	  delete fEvent;//added 1.04.2007
	  fEvent=new AliESDEvent();
	  //	  delete fTree;
	  //fTree=0;
	  //	  delete fEsdFile;
		
	  //ESD data
	  //	  fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	  //	  fTree = (TTree*) fEsdFile->Get("esdTree");			
	  //	  fTree->SetBranchAddress("ESD", &fEvent);			
	  fTree->SetBranchStatus("MuonTracks*",0);
	  fTree->SetBranchStatus("PmdTracks*",0);
	  fTree->SetBranchStatus("TrdTracks*",0);
	  fTree->SetBranchStatus("V0s*",0);
	  fTree->SetBranchStatus("Cascades*",0);
	  fTree->SetBranchStatus("Kinks*",0);
	  fTree->SetBranchStatus("CaloClusters*",0);
	  fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
	  fTree->SetBranchStatus("ESDfriend*",0);
	  fEvent->ReadFromTree(fTree);

	  fNumberofEvent=fTree->GetEntries();
	  cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	  fCurEvent=0;
	  //sim data
	}
      else //no more data to read
	{
	  cout<<"no more files "<<hbtEvent<<endl;
	  fReaderStatus=1;
	  return hbtEvent; 
	}
    }		
  cout<<"starting to read event "<<fCurEvent<<endl;
  fTree->GetEvent(fCurEvent);//getting next event
  cout << "Read event " << fEvent << " from file " << fTree << endl;
  //  vector<int> tLabelTable;//to check labels
	
  hbtEvent = new AliFemtoEvent;
  //setting basic things
  //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
  hbtEvent->SetRunNumber(fEvent->GetRunNumber());
  //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
  hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
  hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
  hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
  hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
  hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
  hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
  hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
  hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
  hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	
  //Vertex
  double fV1[3];
  fEvent->GetVertex()->GetXYZ(fV1);

  AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
  hbtEvent->SetPrimVertPos(vertex);
	
  //starting to reading tracks
  int nofTracks=0;  //number of reconstructed tracks in event
  nofTracks=fEvent->GetNumberOfTracks();
  int realnofTracks=0;//number of track which we use ina analysis
  cout << "Event has " << nofTracks << " tracks " << endl;

  for (int i=0;i<nofTracks;i++)
    {
      bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
		
      AliFemtoTrack* trackCopy = new AliFemtoTrack();	
      const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
      //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();

      trackCopy->SetCharge((short)esdtrack->GetSign());

      //in aliroot we have AliPID 
      //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
      //we use only 5 first
      double esdpid[5];
      esdtrack->GetESDpid(esdpid);
      trackCopy->SetPidProbElectron(esdpid[0]);
      trackCopy->SetPidProbMuon(esdpid[1]);
      trackCopy->SetPidProbPion(esdpid[2]);
      trackCopy->SetPidProbKaon(esdpid[3]);
      trackCopy->SetPidProbProton(esdpid[4]);
						
      double pxyz[3];
      if (fReadInner == true) {
	
	if (esdtrack->GetTPCInnerParam()) {
	  AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
	  param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
	  param->GetPxPyPz(pxyz);//reading noconstarined momentum
	  delete param;

	  AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
	  tInfo->SetPDGPid(211);
	  tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
	  tInfo->SetMass(0.13957);
	  trackCopy->SetHiddenInfo(tInfo);
	}
      }
      if (fConstrained==true)		    
	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
      else
	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
      AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
      trackCopy->SetP(v);//setting momentum
      trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
      const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
      if (ktP.Mag() == 0) {
	delete trackCopy;
	continue;
      }
      const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
      //setting helix I do not if it is ok
      AliFmPhysicalHelixD helix(ktP,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
      trackCopy->SetHelix(helix);
	    	
      trackCopy->SetTrackId(esdtrack->GetID());
      trackCopy->SetFlags(esdtrack->GetStatus());
      trackCopy->SetLabel(esdtrack->GetLabel());
		
      //some stuff which could be useful 
      float impact[2];
      float covimpact[3];
      esdtrack->GetImpactParameters(impact,covimpact);
      trackCopy->SetImpactD(impact[0]);
      trackCopy->SetImpactZ(impact[1]);
      trackCopy->SetCdd(covimpact[0]);
      trackCopy->SetCdz(covimpact[1]);
      trackCopy->SetCzz(covimpact[2]);
      trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
      trackCopy->SetITSncls(esdtrack->GetNcls(0));     
      trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
      trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
      trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
      trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
      trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 

      trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
      trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());

      double pvrt[3];
      fEvent->GetPrimaryVertex()->GetXYZ(pvrt);

      double xtpc[3];
      esdtrack->GetInnerXYZ(xtpc);
      xtpc[2] -= pvrt[2];
      trackCopy->SetNominalTPCEntrancePoint(xtpc);

      esdtrack->GetOuterXYZ(xtpc);
      xtpc[2] -= pvrt[2];
      trackCopy->SetNominalTPCExitPoint(xtpc);

      int indexes[3];
      for (int ik=0; ik<3; ik++) {
	indexes[ik] = esdtrack->GetKinkIndex(ik);
      }
      trackCopy->SetKinkIndexes(indexes);
      //decision if we want this track
      //if we using diffrent labels we want that this label was use for first time 
      //if we use hidden info we want to have match between sim data and ESD
      if (tGoodMomentum==true)
	{
	  hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	  realnofTracks++;//real number of tracks
	}
      else
	{
	  delete  trackCopy;
	}
		
    }

  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event	
  fCurEvent++;	
  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
//   if (fCurEvent== fNumberofEvent)//if end of current file close all
//     {   
//       fTree->Reset(); 
//       delete fTree;
//       fEsdFile->Close();
//     }
  return hbtEvent; 
}
Commit	Line	Data
0215f606	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoEventReaderESD - the reader class for the Alice ESD ///
	4	/// Reads in ESD information and converts it into internal AliFemtoEvent ///
	5	/// Reads in AliESDfriend to create shared hit/quality information ///
	6	/// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl ///
	7	/// Adam Kisiel kisiel@mps.ohio-state.edu ///
	8	/// ///
	9	////////////////////////////////////////////////////////////////////////////////
	10
67427ff7	11	/*
	12	*$Id$
	13	*$Log$
ea77036b	14	*Revision 1.2.2.2 2007/10/04 13:10:52 akisiel
	15	*Add Kink index storageAliFemtoEventReaderESD.cxx AliFemtoTrack.cxx AliFemtoTrack.h
	16	*
	17	*Revision 1.2.2.1 2007/09/30 11:38:59 akisiel
	18	*Adapt the readers to the new AliESDEvent structure
	19	*
	20	*Revision 1.2 2007/05/22 09:01:42 akisiel
	21	*Add the possibiloity to save cut settings in the ROOT file
	22	*
3a74204a	23	*Revision 1.1 2007/05/16 10:22:11 akisiel
	24	*Making the directory structure of AliFemto flat. All files go into one common directory
	25	*
36892005	26	*Revision 1.5 2007/05/03 09:45:20 akisiel
	27	*Fixing Effective C++ warnings
	28	*
0215f606	29	*Revision 1.4 2007/04/27 07:28:34 akisiel
	30	*Remove event number reading due to interface changes
	31	*
a531c6fc	32	*Revision 1.3 2007/04/27 07:25:16 akisiel
	33	*Make revisions needed for compilation from the main AliRoot tree
	34	*
b2f60a91	35	*Revision 1.1.1.1 2007/04/25 15:38:41 panos
	36	*Importing the HBT code dir
	37	*
67427ff7	38	*/
	39
	40	#include "AliFemtoEventReaderESD.h"
	41
	42	#include "TFile.h"
	43	#include "TTree.h"
ea77036b	44	#include "AliESDEvent.h"
67427ff7	45	#include "AliESDtrack.h"
ea77036b	46	#include "AliESDVertex.h"
67427ff7	47
	48	//#include "TSystem.h"
	49
d0e92d9a	50	#include "AliFmPhysicalHelixD.h"
d0e92d9a	51	#include "AliFmThreeVectorF.h"
67427ff7	52
d0e92d9a	53	#include "SystemOfUnits.h"
67427ff7	54
d0e92d9a	55	#include "AliFemtoEvent.h"
ea77036b	56	#include "AliFemtoModelHiddenInfo.h"
67427ff7	57
	58	ClassImp(AliFemtoEventReaderESD)
	59
	60	#if !(ST_NO_NAMESPACES)
	61	using namespace units;
	62	#endif
	63
	64	using namespace std;
	65	//____________________________
	66	//constructor with 0 parameters , look at default settings
	67	AliFemtoEventReaderESD::AliFemtoEventReaderESD():
	68	fInputFile(" "),
	69	fFileName(" "),
	70	fConstrained(true),
ea77036b	71	fReadInner(false),
67427ff7	72	fNumberofEvent(0),
67427ff7	73	fCurEvent(0),
67427ff7	74	fTree(0x0),
67427ff7	75	fEsdFile(0x0),
ea77036b	76	fEvent(0x0)
67427ff7	77	{
d0e92d9a	78	// default constructor
67427ff7	79	}
67427ff7	80
0215f606	81	AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
fcda1d4e	82	AliFemtoEventReader(aReader),
0215f606	83	fInputFile(" "),
	84	fFileName(" "),
	85	fConstrained(true),
ea77036b	86	fReadInner(false),
0215f606	87	fNumberofEvent(0),
0215f606	88	fCurEvent(0),
0215f606	89	fTree(0x0),
0215f606	90	fEsdFile(0x0),
ea77036b	91	fEvent(0x0)
0215f606	92	{
d0e92d9a	93	// copy constructor
0215f606	94	fInputFile = aReader.fInputFile;
	95	fFileName = aReader.fFileName;
	96	fConstrained = aReader.fConstrained;
ea77036b	97	fReadInner = aReader.fReadInner;
0215f606	98	fNumberofEvent = aReader.fNumberofEvent;
0215f606	99	fCurEvent = aReader.fCurEvent;
ea77036b	100	// fTree = aReader.fTree->CloneTree();
d0e92d9a	101	// fEvent = new AliESD(*aReader.fEvent);
ea77036b	102	fEvent = new AliESDEvent();
0215f606	103	fEsdFile = new TFile(aReader.fEsdFile->GetName());
0215f606	104	}
67427ff7	105	//__________________
	106	//Destructor
	107	AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
	108	{
d0e92d9a	109	// destructor
67427ff7	110	//delete fListOfFiles;
	111	delete fTree;
	112	delete fEvent;
	113	delete fEsdFile;
67427ff7	114	}
67427ff7	115
0215f606	116	//__________________
	117	AliFemtoEventReaderESD& AliFemtoEventReaderESD::operator=(const AliFemtoEventReaderESD& aReader)
	118	{
d0e92d9a	119	// assignment operator
0215f606	120	if (this == &aReader)
	121	return *this;
	122
	123	fInputFile = aReader.fInputFile;
	124	fFileName = aReader.fFileName;
	125	fConstrained = aReader.fConstrained;
ea77036b	126	fReadInner = aReader.fReadInner;
0215f606	127	fNumberofEvent = aReader.fNumberofEvent;
0215f606	128	fCurEvent = aReader.fCurEvent;
0215f606	129	if (fTree) delete fTree;
ea77036b	130	// fTree = aReader.fTree->CloneTree();
0215f606	131	if (fEvent) delete fEvent;
ea77036b	132	fEvent = new AliESDEvent();
0215f606	133	if (fEsdFile) delete fEsdFile;
	134	fEsdFile = new TFile(aReader.fEsdFile->GetName());
	135
0215f606	136	return *this;
0215f606	137	}
67427ff7	138	//__________________
	139	AliFemtoString AliFemtoEventReaderESD::Report()
	140	{
d0e92d9a	141	// create reader report
67427ff7	142	AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
	143	return temp;
	144	}
	145
	146	//__________________
67427ff7	147	void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
67427ff7	148	{
d0e92d9a	149	//setting the name of file where names of ESD file are written
d0e92d9a	150	//it takes only this files which have good trees
67427ff7	151	char buffer[256];
	152	fInputFile=string(inputFile);
	153	cout<<"Input File set on "<<fInputFile<<endl;
	154	ifstream infile(inputFile);
ea77036b	155
	156	fTree = new TChain("esdTree");
	157
67427ff7	158	if(infile.good()==true)
	159	{
	160	//checking if all give files have good tree inside
	161	while (infile.eof()==false)
	162	{
	163	infile.getline(buffer,256);
	164	//ifstream test_file(buffer);
	165	TFile *esdFile=TFile::Open(buffer,"READ");
	166	if (esdFile!=0x0)
	167	{
	168	TTree* tree = (TTree*) esdFile->Get("esdTree");
	169	if (tree!=0x0)
	170	{
	171	cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
ea77036b	172	fTree->AddFile(buffer);
67427ff7	173	delete tree;
	174	}
	175	esdFile->Close();
	176	}
	177	delete esdFile;
	178	}
	179	}
	180	}
	181
67427ff7	182	void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
	183	{
	184	fConstrained=constrained;
	185	}
	186
	187	bool AliFemtoEventReaderESD::GetConstrained() const
	188	{
	189	return fConstrained;
	190	}
	191
ea77036b	192	void AliFemtoEventReaderESD::SetReadTPCInner(const bool readinner)
	193	{
	194	fReadInner=readinner;
	195	}
	196
	197	bool AliFemtoEventReaderESD::GetReadTPCInner() const
	198	{
	199	return fReadInner;
	200	}
	201
67427ff7	202	AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
67427ff7	203	{
d0e92d9a	204	// read in a next hbt event from the chain
	205	// convert it to AliFemtoEvent and return
	206	// for further analysis
67427ff7	207	AliFemtoEvent *hbtEvent = 0;
67427ff7	208
	209	if (fCurEvent==fNumberofEvent)//open next file
	210	{
ea77036b	211	if(fNumberofEvent==0)
67427ff7	212	{
ea77036b	213	// delete fEvent;//added 1.04.2007
ea77036b	214	fEvent=new AliESDEvent();
1e1d2996	215	// delete fTree;
67427ff7	216	//fTree=0;
ea77036b	217	// delete fEsdFile;
67427ff7	218
67427ff7	219	//ESD data
ea77036b	220	// fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	221	// fTree = (TTree*) fEsdFile->Get("esdTree");
	222	// fTree->SetBranchAddress("ESD", &fEvent);
	223	fTree->SetBranchStatus("MuonTracks*",0);
	224	fTree->SetBranchStatus("PmdTracks*",0);
	225	fTree->SetBranchStatus("TrdTracks*",0);
	226	fTree->SetBranchStatus("V0s*",0);
	227	fTree->SetBranchStatus("Cascades*",0);
	228	fTree->SetBranchStatus("Kinks*",0);
	229	fTree->SetBranchStatus("CaloClusters*",0);
	230	fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
	231	fTree->SetBranchStatus("ESDfriend*",0);
	232	fEvent->ReadFromTree(fTree);
1e1d2996	233
67427ff7	234	fNumberofEvent=fTree->GetEntries();
	235	cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	236	fCurEvent=0;
	237	//sim data
	238	}
	239	else //no more data to read
	240	{
	241	cout<<"no more files "<<hbtEvent<<endl;
	242	fReaderStatus=1;
	243	return hbtEvent;
	244	}
	245	}
	246	cout<<"starting to read event "<<fCurEvent<<endl;
	247	fTree->GetEvent(fCurEvent);//getting next event
ea77036b	248	cout << "Read event " << fEvent << " from file " << fTree << endl;
d0e92d9a	249	// vector<int> tLabelTable;//to check labels
67427ff7	250
	251	hbtEvent = new AliFemtoEvent;
	252	//setting basic things
a531c6fc	253	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7	254	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	255	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	256	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	257	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	258	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	259	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	260	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	261	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	262	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	263	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	264	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	265
	266	//Vertex
	267	double fV1[3];
	268	fEvent->GetVertex()->GetXYZ(fV1);
	269
	270	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	271	hbtEvent->SetPrimVertPos(vertex);
	272
	273	//starting to reading tracks
	274	int nofTracks=0; //number of reconstructed tracks in event
	275	nofTracks=fEvent->GetNumberOfTracks();
	276	int realnofTracks=0;//number of track which we use ina analysis
ea77036b	277	cout << "Event has " << nofTracks << " tracks " << endl;
67427ff7	278
	279	for (int i=0;i<nofTracks;i++)
	280	{
d0e92d9a	281	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
67427ff7	282
	283	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	284	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
0215f606	285	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
67427ff7	286
	287	trackCopy->SetCharge((short)esdtrack->GetSign());
	288
	289	//in aliroot we have AliPID
	290	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	291	//we use only 5 first
	292	double esdpid[5];
	293	esdtrack->GetESDpid(esdpid);
	294	trackCopy->SetPidProbElectron(esdpid[0]);
	295	trackCopy->SetPidProbMuon(esdpid[1]);
	296	trackCopy->SetPidProbPion(esdpid[2]);
	297	trackCopy->SetPidProbKaon(esdpid[3]);
	298	trackCopy->SetPidProbProton(esdpid[4]);
	299
	300	double pxyz[3];
ea77036b	301	if (fReadInner == true) {
	302
	303	if (esdtrack->GetTPCInnerParam()) {
	304	AliExternalTrackParam param = new AliExternalTrackParam(esdtrack->GetTPCInnerParam());
	305	param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
	306	param->GetPxPyPz(pxyz);//reading noconstarined momentum
	307	delete param;
	308
	309	AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
	310	tInfo->SetPDGPid(211);
	311	tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
	312	tInfo->SetMass(0.13957);
	313	trackCopy->SetHiddenInfo(tInfo);
	314	}
	315	}
67427ff7	316	if (fConstrained==true)
d0e92d9a	317	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
67427ff7	318	else
d0e92d9a	319	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
67427ff7	320	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	321	trackCopy->SetP(v);//setting momentum
	322	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
d0e92d9a	323	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
69c1c8ff	324	if (ktP.Mag() == 0) {
67427ff7	325	delete trackCopy;
	326	continue;
	327	}
	328	const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
	329	//setting helix I do not if it is ok
d0e92d9a	330	AliFmPhysicalHelixD helix(ktP,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
67427ff7	331	trackCopy->SetHelix(helix);
	332
	333	trackCopy->SetTrackId(esdtrack->GetID());
	334	trackCopy->SetFlags(esdtrack->GetStatus());
ea77036b	335	trackCopy->SetLabel(esdtrack->GetLabel());
67427ff7	336
	337	//some stuff which could be useful
	338	float impact[2];
	339	float covimpact[3];
	340	esdtrack->GetImpactParameters(impact,covimpact);
	341	trackCopy->SetImpactD(impact[0]);
	342	trackCopy->SetImpactZ(impact[1]);
	343	trackCopy->SetCdd(covimpact[0]);
	344	trackCopy->SetCdz(covimpact[1]);
	345	trackCopy->SetCzz(covimpact[2]);
	346	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
	347	trackCopy->SetITSncls(esdtrack->GetNcls(0));
	348	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
	349	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
	350	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
	351	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
	352	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
	353
ea77036b	354	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
ea77036b	355	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
67427ff7	356
0b3bd1ac	357	double pvrt[3];
	358	fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
	359
	360	double xtpc[3];
	361	esdtrack->GetInnerXYZ(xtpc);
	362	xtpc[2] -= pvrt[2];
	363	trackCopy->SetNominalTPCEntrancePoint(xtpc);
	364
	365	esdtrack->GetOuterXYZ(xtpc);
	366	xtpc[2] -= pvrt[2];
	367	trackCopy->SetNominalTPCExitPoint(xtpc);
	368
ea77036b	369	int indexes[3];
	370	for (int ik=0; ik<3; ik++) {
	371	indexes[ik] = esdtrack->GetKinkIndex(ik);
	372	}
	373	trackCopy->SetKinkIndexes(indexes);
67427ff7	374	//decision if we want this track
	375	//if we using diffrent labels we want that this label was use for first time
	376	//if we use hidden info we want to have match between sim data and ESD
d0e92d9a	377	if (tGoodMomentum==true)
67427ff7	378	{
	379	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	380	realnofTracks++;//real number of tracks
	381	}
	382	else
	383	{
	384	delete trackCopy;
	385	}
	386
	387	}
	388
	389	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
	390	fCurEvent++;
	391	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
ea77036b	392	// if (fCurEvent== fNumberofEvent)//if end of current file close all
	393	// {
	394	// fTree->Reset();
	395	// delete fTree;
	396	// fEsdFile->Close();
	397	// }
67427ff7	398	return hbtEvent;
	399	}
	400
	401
	402
	403
	404
	405
	406
	407
	408