[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Reader / 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.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 "AliESD.h"
#include "AliESDtrack.h"

//#include "TSystem.h"

#include "Infrastructure/AliFmPhysicalHelixD.h"
#include "Infrastructure/AliFmThreeVectorF.h"

#include "Base/SystemOfUnits.h"

#include "Infrastructure/AliFemtoEvent.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),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fListOfFiles(0x0),
  fTree(0x0),
  fEvent(0x0),
  fEsdFile(0x0),
  fEventFriend(0),
  fSharedList(0x0),
  fClusterPerPadrow(0x0)
{
  fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad] = new list<Int_t>();
  }
  fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad] = new list<Int_t>();
  }
}

AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fListOfFiles(0x0),
  fTree(0x0),
  fEvent(0x0),
  fEsdFile(0x0),
  fEventFriend(0),
  fSharedList(0x0),
  fClusterPerPadrow(0x0)
{
  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fTree = aReader.fTree->CloneTree();
  fEvent = new AliESD(*aReader.fEvent);
  fEsdFile = new TFile(aReader.fEsdFile->GetName());
  fEventFriend = aReader.fEventFriend;
  fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad] = new list<Int_t>();
    list<Int_t>::iterator iter;
    for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
      fClusterPerPadrow[tPad]->push_back(*iter);
    }
  }
  fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad] = new list<Int_t>();
    list<Int_t>::iterator iter;
    for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
      fSharedList[tPad]->push_back(*iter);
    }
  }
  for (unsigned int veciter = 0; veciter<aReader.fListOfFiles.size(); veciter++)
    fListOfFiles.push_back(aReader.fListOfFiles[veciter]);
}
//__________________
//Destructor
AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
{
  //delete fListOfFiles;
  delete fTree;
  delete fEvent;
  delete fEsdFile;

  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad]->clear();
    delete fClusterPerPadrow[tPad];
  }
  delete [] fClusterPerPadrow;
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad]->clear();
    delete fSharedList[tPad];
  }
  delete [] fSharedList;
}

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

  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  if (fTree) delete fTree;
  fTree = aReader.fTree->CloneTree();
  if (fEvent) delete fEvent;
  fEvent = new AliESD(*aReader.fEvent);
  if (fEsdFile) delete fEsdFile;
  fEsdFile = new TFile(aReader.fEsdFile->GetName());

  fEventFriend = aReader.fEventFriend;
  
  if (fClusterPerPadrow) {
    for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
      fClusterPerPadrow[tPad]->clear();
      delete fClusterPerPadrow[tPad];
    }
    delete [] fClusterPerPadrow;
  }
  
  if (fSharedList) {
    for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
      fSharedList[tPad]->clear();
      delete fSharedList[tPad];
    }
    delete [] fSharedList;
  }

  fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad] = new list<Int_t>();
    list<Int_t>::iterator iter;
    for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
      fClusterPerPadrow[tPad]->push_back(*iter);
    }
  }
  fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad] = new list<Int_t>();
    list<Int_t>::iterator iter;
    for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
      fSharedList[tPad]->push_back(*iter);
    }
  }
  for (unsigned int veciter = 0; veciter<aReader.fListOfFiles.size(); veciter++)
    fListOfFiles.push_back(aReader.fListOfFiles[veciter]);
  
  return *this;
}
//__________________
AliFemtoString AliFemtoEventReaderESD::Report()
{
  AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
  return temp;
}

//__________________
//setting the name of file where names of ESD file are written 
//it takes only this files which have good trees
void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
{
  char buffer[256];
  fInputFile=string(inputFile);
  cout<<"Input File set on "<<fInputFile<<endl;
  ifstream infile(inputFile);
  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;
		  fListOfFiles.push_back(string(buffer));
		  delete tree;
		}
	      esdFile->Close();	
	    }
	  delete esdFile;
	}
    }
}

//setting the next file to read	
bool AliFemtoEventReaderESD::GetNextFile()
{ 	
  if (fCurFile>=fListOfFiles.size())
    return false;
  fFileName=fListOfFiles.at(fCurFile);	
  cout<<"FileName set on "<<fFileName<<" "<<fCurFile<<endl;

  fCurFile++;
  return true;
}
void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}

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

AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
{
  AliFemtoEvent *hbtEvent = 0;
  string tFriendFileName;

  if (fCurEvent==fNumberofEvent)//open next file  
    {
      cout<<"next file"<<endl;
      if(GetNextFile())	
	{
	  delete fEventFriend;
	  fEventFriend = 0;
	  delete fEvent;//added 1.04.2007
	  fEvent=new AliESD();
	  //	  delete fTree;
	  //fTree=0;
	  delete fEsdFile;
		
	  //ESD data
	  fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	  fTree = (TTree*) fEsdFile->Get("esdTree");			
	  fTree->SetBranchAddress("ESD", &fEvent);			

	  // Attach the friend tree with additional information
 	  tFriendFileName = fFileName;
 	  tFriendFileName.insert(tFriendFileName.find("s.root"),"friend");
 	  cout << "Reading friend " << tFriendFileName.c_str() << endl;;
  	  fTree->AddFriend("esdFriendTree",tFriendFileName.c_str());
  	  fTree->SetBranchAddress("ESDfriend",&fEventFriend);

// 	  chain->SetBranchStatus("*",0);
// 	  chain->SetBranchStatus("fUniqueID",1);
// 	  chain->SetBranchStatus("fTracks",1);
// 	  chain->SetBranchStatus("fTracks.*",1);
// 	  chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
	  fTree->SetBranchStatus("fTracks.fCalibContainer",0);


	  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
  fEvent->SetESDfriend(fEventFriend);
  vector<int> label_table;//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

  // Clear the shared cluster list
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad]->clear();
  }
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad]->clear();
  }


  for (int i=0;i<nofTracks;i++) {
    const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track

    list<Int_t>::iterator tClustIter;

    Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
    Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
    for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
      if (tTrackIndices[tNcl] >= 0) {
	tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
	  if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
	  fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
	}
	else {
	  fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
	}
      }
    }
      
  }

  for (int i=0;i<nofTracks;i++)
    {
      bool  good_momentum=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 (fConstrained==true)		    
	good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
      else
	good_momentum=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 p(pxyz[0],pxyz[1],pxyz[2]);
      if (p.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(p,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()); 

      // Fill cluster per padrow information
      Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
      Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
      for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	if (tTrackIndices[tNcl] > 0)
	  trackCopy->SetTPCcluster(tNcl, 1);
	else
	  trackCopy->SetTPCcluster(tNcl, 0);
      }
      
      // Fill shared cluster information
      list<Int_t>::iterator tClustIter;

      for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	if (tTrackIndices[tNcl] > 0) {
	  tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
	  if (tClustIter != fSharedList[tNcl]->end()) {
	    trackCopy->SetTPCshared(tNcl, 1);
	    cout << "Event next" <<  endl;
	    cout << "Track: " << i << endl;
	    cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
	  }
	  else {
	    trackCopy->SetTPCshared(tNcl, 0);
	  }
	}
      }

      //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 (good_momentum==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$
0215f606	14	*Revision 1.4 2007/04/27 07:28:34 akisiel
	15	*Remove event number reading due to interface changes
	16	*
a531c6fc	17	*Revision 1.3 2007/04/27 07:25:16 akisiel
	18	*Make revisions needed for compilation from the main AliRoot tree
	19	*
b2f60a91	20	*Revision 1.1.1.1 2007/04/25 15:38:41 panos
	21	*Importing the HBT code dir
	22	*
67427ff7	23	*/
	24
	25	#include "AliFemtoEventReaderESD.h"
	26
	27	#include "TFile.h"
	28	#include "TTree.h"
	29	#include "AliESD.h"
	30	#include "AliESDtrack.h"
	31
	32	//#include "TSystem.h"
	33
b2f60a91	34	#include "Infrastructure/AliFmPhysicalHelixD.h"
b2f60a91	35	#include "Infrastructure/AliFmThreeVectorF.h"
67427ff7	36
b2f60a91	37	#include "Base/SystemOfUnits.h"
67427ff7	38
b2f60a91	39	#include "Infrastructure/AliFemtoEvent.h"
67427ff7	40
	41	ClassImp(AliFemtoEventReaderESD)
	42
	43	#if !(ST_NO_NAMESPACES)
	44	using namespace units;
	45	#endif
	46
	47	using namespace std;
	48	//____________________________
	49	//constructor with 0 parameters , look at default settings
	50	AliFemtoEventReaderESD::AliFemtoEventReaderESD():
	51	fInputFile(" "),
	52	fFileName(" "),
	53	fConstrained(true),
	54	fNumberofEvent(0),
	55	fCurEvent(0),
	56	fCurFile(0),
0215f606	57	fListOfFiles(0x0),
67427ff7	58	fTree(0x0),
	59	fEvent(0x0),
	60	fEsdFile(0x0),
0215f606	61	fEventFriend(0),
	62	fSharedList(0x0),
	63	fClusterPerPadrow(0x0)
67427ff7	64	{
	65	fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	66	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	67	fClusterPerPadrow[tPad] = new list<Int_t>();
	68	}
	69	fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	70	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	71	fSharedList[tPad] = new list<Int_t>();
	72	}
	73	}
	74
0215f606	75	AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
	76	fInputFile(" "),
	77	fFileName(" "),
	78	fConstrained(true),
	79	fNumberofEvent(0),
	80	fCurEvent(0),
	81	fCurFile(0),
	82	fListOfFiles(0x0),
	83	fTree(0x0),
	84	fEvent(0x0),
	85	fEsdFile(0x0),
	86	fEventFriend(0),
	87	fSharedList(0x0),
	88	fClusterPerPadrow(0x0)
	89	{
	90	fInputFile = aReader.fInputFile;
	91	fFileName = aReader.fFileName;
	92	fConstrained = aReader.fConstrained;
	93	fNumberofEvent = aReader.fNumberofEvent;
	94	fCurEvent = aReader.fCurEvent;
	95	fCurFile = aReader.fCurFile;
	96	fTree = aReader.fTree->CloneTree();
	97	fEvent = new AliESD(*aReader.fEvent);
	98	fEsdFile = new TFile(aReader.fEsdFile->GetName());
	99	fEventFriend = aReader.fEventFriend;
	100	fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	101	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	102	fClusterPerPadrow[tPad] = new list<Int_t>();
	103	list<Int_t>::iterator iter;
	104	for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
	105	fClusterPerPadrow[tPad]->push_back(*iter);
	106	}
	107	}
	108	fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	109	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	110	fSharedList[tPad] = new list<Int_t>();
	111	list<Int_t>::iterator iter;
	112	for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
	113	fSharedList[tPad]->push_back(*iter);
	114	}
	115	}
	116	for (unsigned int veciter = 0; veciter<aReader.fListOfFiles.size(); veciter++)
	117	fListOfFiles.push_back(aReader.fListOfFiles[veciter]);
	118	}
67427ff7	119	//__________________
	120	//Destructor
	121	AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
	122	{
	123	//delete fListOfFiles;
	124	delete fTree;
	125	delete fEvent;
	126	delete fEsdFile;
	127
	128	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	129	fClusterPerPadrow[tPad]->clear();
	130	delete fClusterPerPadrow[tPad];
	131	}
	132	delete [] fClusterPerPadrow;
	133	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	134	fSharedList[tPad]->clear();
	135	delete fSharedList[tPad];
	136	}
	137	delete [] fSharedList;
	138	}
	139
0215f606	140	//__________________
	141	AliFemtoEventReaderESD& AliFemtoEventReaderESD::operator=(const AliFemtoEventReaderESD& aReader)
	142	{
	143	if (this == &aReader)
	144	return *this;
	145
	146	fInputFile = aReader.fInputFile;
	147	fFileName = aReader.fFileName;
	148	fConstrained = aReader.fConstrained;
	149	fNumberofEvent = aReader.fNumberofEvent;
	150	fCurEvent = aReader.fCurEvent;
	151	fCurFile = aReader.fCurFile;
	152	if (fTree) delete fTree;
	153	fTree = aReader.fTree->CloneTree();
	154	if (fEvent) delete fEvent;
	155	fEvent = new AliESD(*aReader.fEvent);
	156	if (fEsdFile) delete fEsdFile;
	157	fEsdFile = new TFile(aReader.fEsdFile->GetName());
	158
	159	fEventFriend = aReader.fEventFriend;
	160
	161	if (fClusterPerPadrow) {
	162	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	163	fClusterPerPadrow[tPad]->clear();
	164	delete fClusterPerPadrow[tPad];
	165	}
	166	delete [] fClusterPerPadrow;
	167	}
	168
	169	if (fSharedList) {
	170	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	171	fSharedList[tPad]->clear();
	172	delete fSharedList[tPad];
	173	}
	174	delete [] fSharedList;
	175	}
	176
	177	fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	178	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	179	fClusterPerPadrow[tPad] = new list<Int_t>();
	180	list<Int_t>::iterator iter;
	181	for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
	182	fClusterPerPadrow[tPad]->push_back(*iter);
	183	}
	184	}
	185	fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	186	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	187	fSharedList[tPad] = new list<Int_t>();
	188	list<Int_t>::iterator iter;
	189	for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
	190	fSharedList[tPad]->push_back(*iter);
	191	}
	192	}
	193	for (unsigned int veciter = 0; veciter<aReader.fListOfFiles.size(); veciter++)
	194	fListOfFiles.push_back(aReader.fListOfFiles[veciter]);
	195
	196	return *this;
	197	}
67427ff7	198	//__________________
	199	AliFemtoString AliFemtoEventReaderESD::Report()
	200	{
	201	AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
	202	return temp;
	203	}
	204
	205	//__________________
	206	//setting the name of file where names of ESD file are written
	207	//it takes only this files which have good trees
	208	void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
	209	{
	210	char buffer[256];
	211	fInputFile=string(inputFile);
	212	cout<<"Input File set on "<<fInputFile<<endl;
	213	ifstream infile(inputFile);
	214	if(infile.good()==true)
	215	{
	216	//checking if all give files have good tree inside
	217	while (infile.eof()==false)
	218	{
	219	infile.getline(buffer,256);
	220	//ifstream test_file(buffer);
	221	TFile *esdFile=TFile::Open(buffer,"READ");
	222	if (esdFile!=0x0)
	223	{
	224	TTree* tree = (TTree*) esdFile->Get("esdTree");
	225	if (tree!=0x0)
	226	{
	227	cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
	228	fListOfFiles.push_back(string(buffer));
	229	delete tree;
	230	}
	231	esdFile->Close();
	232	}
	233	delete esdFile;
	234	}
	235	}
	236	}
	237
	238	//setting the next file to read
	239	bool AliFemtoEventReaderESD::GetNextFile()
	240	{
	241	if (fCurFile>=fListOfFiles.size())
	242	return false;
	243	fFileName=fListOfFiles.at(fCurFile);
	244	cout<<"FileName set on "<<fFileName<<" "<<fCurFile<<endl;
	245
	246	fCurFile++;
	247	return true;
	248	}
	249	void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
	250	{
	251	fConstrained=constrained;
	252	}
	253
	254	bool AliFemtoEventReaderESD::GetConstrained() const
	255	{
	256	return fConstrained;
	257	}
	258
	259	AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
	260	{
	261	AliFemtoEvent *hbtEvent = 0;
262	string tFriendFileName;
263
264	if (fCurEvent==fNumberofEvent)//open next file
265	{
266	cout<<"next file"<<endl;
267	if(GetNextFile())
268	{
269	delete fEventFriend;
270	fEventFriend = 0;
271	delete fEvent;//added 1.04.2007
272	fEvent=new AliESD();
1e1d2996	273	// delete fTree;
67427ff7	274	//fTree=0;
	275	delete fEsdFile;
	276
	277	//ESD data
	278	fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	279	fTree = (TTree*) fEsdFile->Get("esdTree");
	280	fTree->SetBranchAddress("ESD", &fEvent);
	281
	282	// Attach the friend tree with additional information
	283	tFriendFileName = fFileName;
	284	tFriendFileName.insert(tFriendFileName.find("s.root"),"friend");
	285	cout << "Reading friend " << tFriendFileName.c_str() << endl;;
	286	fTree->AddFriend("esdFriendTree",tFriendFileName.c_str());
	287	fTree->SetBranchAddress("ESDfriend",&fEventFriend);
	288
1e1d2996	289	// chain->SetBranchStatus("*",0);
	290	// chain->SetBranchStatus("fUniqueID",1);
	291	// chain->SetBranchStatus("fTracks",1);
	292	// chain->SetBranchStatus("fTracks.*",1);
	293	// chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
	294	fTree->SetBranchStatus("fTracks.fCalibContainer",0);
	295
	296
67427ff7	297	fNumberofEvent=fTree->GetEntries();
	298	cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	299	fCurEvent=0;
	300	//sim data
	301	}
	302	else //no more data to read
	303	{
	304	cout<<"no more files "<<hbtEvent<<endl;
	305	fReaderStatus=1;
	306	return hbtEvent;
	307	}
	308	}
	309	cout<<"starting to read event "<<fCurEvent<<endl;
	310	fTree->GetEvent(fCurEvent);//getting next event
	311	fEvent->SetESDfriend(fEventFriend);
	312	vector<int> label_table;//to check labels
	313
	314	hbtEvent = new AliFemtoEvent;
	315	//setting basic things
a531c6fc	316	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7	317	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	318	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	319	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	320	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	321	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	322	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	323	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	324	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	325	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	326	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	327	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	328
	329	//Vertex
	330	double fV1[3];
	331	fEvent->GetVertex()->GetXYZ(fV1);
	332
	333	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	334	hbtEvent->SetPrimVertPos(vertex);
	335
	336	//starting to reading tracks
	337	int nofTracks=0; //number of reconstructed tracks in event
	338	nofTracks=fEvent->GetNumberOfTracks();
	339	int realnofTracks=0;//number of track which we use ina analysis
	340
	341	// Clear the shared cluster list
	342	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	343	fClusterPerPadrow[tPad]->clear();
	344	}
	345	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	346	fSharedList[tPad]->clear();
	347	}
	348
	349
	350	for (int i=0;i<nofTracks;i++) {
	351	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	352
	353	list<Int_t>::iterator tClustIter;
	354
	355	Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
	356	Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
	357	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	358	if (tTrackIndices[tNcl] >= 0) {
	359	tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
0215f606	360	if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
67427ff7	361	fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
	362	}
	363	else {
	364	fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
	365	}
	366	}
	367	}
	368
	369	}
	370
	371	for (int i=0;i<nofTracks;i++)
	372	{
	373	bool good_momentum=true; //flaga to chcek if we can read momentum of this track
	374
	375	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	376	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
0215f606	377	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
67427ff7	378
	379	trackCopy->SetCharge((short)esdtrack->GetSign());
	380
	381	//in aliroot we have AliPID
	382	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	383	//we use only 5 first
	384	double esdpid[5];
	385	esdtrack->GetESDpid(esdpid);
	386	trackCopy->SetPidProbElectron(esdpid[0]);
	387	trackCopy->SetPidProbMuon(esdpid[1]);
	388	trackCopy->SetPidProbPion(esdpid[2]);
	389	trackCopy->SetPidProbKaon(esdpid[3]);
	390	trackCopy->SetPidProbProton(esdpid[4]);
	391
	392	double pxyz[3];
	393	if (fConstrained==true)
	394	good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
	395	else
	396	good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
	397	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	398	trackCopy->SetP(v);//setting momentum
	399	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	400	const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
	401	if (p.mag() == 0) {
	402	delete trackCopy;
	403	continue;
	404	}
	405	const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
	406	//setting helix I do not if it is ok
	407	AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
	408	trackCopy->SetHelix(helix);
	409
	410	trackCopy->SetTrackId(esdtrack->GetID());
	411	trackCopy->SetFlags(esdtrack->GetStatus());
	412	//trackCopy->SetLabel(esdtrack->GetLabel());
	413
	414	//some stuff which could be useful
	415	float impact[2];
	416	float covimpact[3];
	417	esdtrack->GetImpactParameters(impact,covimpact);
	418	trackCopy->SetImpactD(impact[0]);
	419	trackCopy->SetImpactZ(impact[1]);
	420	trackCopy->SetCdd(covimpact[0]);
	421	trackCopy->SetCdz(covimpact[1]);
	422	trackCopy->SetCzz(covimpact[2]);
	423	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
	424	trackCopy->SetITSncls(esdtrack->GetNcls(0));
	425	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
	426	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
	427	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
	428	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
	429	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
	430
	431	// Fill cluster per padrow information
	432	Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
	433	Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
	434	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	435	if (tTrackIndices[tNcl] > 0)
	436	trackCopy->SetTPCcluster(tNcl, 1);
	437	else
	438	trackCopy->SetTPCcluster(tNcl, 0);
	439	}
	440
	441	// Fill shared cluster information
442	list<Int_t>::iterator tClustIter;
443
444	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
445	if (tTrackIndices[tNcl] > 0) {
446	tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
447	if (tClustIter != fSharedList[tNcl]->end()) {
448	trackCopy->SetTPCshared(tNcl, 1);
449	cout << "Event next" << endl;
450	cout << "Track: " << i << endl;
451	cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
452	}
453	else {
454	trackCopy->SetTPCshared(tNcl, 0);
455	}
456	}
457	}
458
459	//decision if we want this track
460	//if we using diffrent labels we want that this label was use for first time
461	//if we use hidden info we want to have match between sim data and ESD
462	if (good_momentum==true)
463	{
464	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
465	realnofTracks++;//real number of tracks
466	}
467	else
468	{
469	delete trackCopy;
470	}
471
472	}
473
474	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
475	fCurEvent++;
476	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
477	if (fCurEvent== fNumberofEvent)//if end of current file close all
478	{
479	fTree->Reset();
480	delete fTree;
481	fEsdFile->Close();
482	}
483	return hbtEvent;
484	}
485
486
487
488
489
490
491
492
493