[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoEventReaderESDKine.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoEventReaderESDKine - 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.1  2007/05/25 12:42:54  akisiel
 *Adding a reader for the Kine information
 *
 *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 "AliFemtoEventReaderESDKine.h"

#include "TFile.h"
#include "TChain.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliStack.h"
//#include "AliAODParticle.h"
#include "TParticle.h"

//#include "TSystem.h"

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

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"

#include "TMath.h"
#include "TParticle.h"
#include "AliFemtoModelHiddenInfo.h"
#include "AliFemtoModelGlobalHiddenInfo.h"
#include "AliGenHijingEventHeader.h"

ClassImp(AliFemtoEventReaderESDKine)

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

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

AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
  AliFemtoEventReader(),
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurRLEvent(0),
  fTree(0x0),
  fEvent(0x0),
  fRunLoader(0x0)
{
  // copy constructor
  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fEvent = new AliESDEvent();
}
//__________________
//Destructor
AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
{
  // destructor
  //delete fListOfFiles;
  delete fTree;
  delete fEvent;
  if (fRunLoader) delete fRunLoader;
}

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

  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurRLEvent = aReader.fCurRLEvent;
  if (fTree) delete fTree;
  //  fTree = aReader.fTree->CloneTree();
  if (fEvent) delete fEvent;
  fEvent = new AliESDEvent();
  if (fRunLoader) delete fRunLoader;
  fRunLoader = new AliRunLoader(*aReader.fRunLoader);

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

//__________________
void AliFemtoEventReaderESDKine::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 AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}

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

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

  if (fCurEvent==fNumberofEvent)//open next file  
    {
      if (fNumberofEvent == 0) {
	fEvent=new AliESDEvent();
		
	  //ESD data
// 	  fEsdFile=TFile::Open(fFileName.c_str(),"READ");
// 	  fTree = (TTree*) fEsdFile->Get("esdTree");			

	  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);

// 	  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();

	if (fNumberofEvent == 0) {
	  cout<<"no event in input "<<endl;
	  fReaderStatus=1;
	  return hbtEvent; 
	}

	cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
	fCurEvent=0;
	// simulation data reading setup
	
      }
      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
  //  vector<int> tLabelTable;//to check labels
  
  cout << "fFileName is " << fFileName.Data() << endl;
  cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
  if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
    fFileName = fTree->GetCurrentFile()->GetName();
    tGAliceFilename = fFileName;
    tGAliceFilename.ReplaceAll("AliESDs","galice");
    cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
    if (fRunLoader) delete fRunLoader;
    fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
    if (fRunLoader==0x0)
      {
	cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
	exit(0);
      }
    if(fRunLoader->LoadHeader())
      {
	cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
	exit(0);
      }
    fRunLoader->LoadKinematics();
    fCurRLEvent = 0;
  }

  fRunLoader->GetEvent(fCurRLEvent);
  AliStack* tStack = 0x0;
  tStack = fRunLoader->Stack();
	
  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];
  double fVCov[6];
  if (fUseTPCOnly) {
    fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
    fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
    if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
      fVCov[4] = -1001.0;
  }
  else {
    fEvent->GetPrimaryVertex()->GetXYZ(fV1);
    fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
    if (!fEvent->GetPrimaryVertex()->GetStatus())
      fVCov[4] = -1001.0;
  }

  AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
  hbtEvent->SetPrimVertPos(vertex);
  hbtEvent->SetPrimVertCov(fVCov);

  AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
	
  Double_t tReactionPlane = 0;
  if (hdh)
    {
      tReactionPlane = hdh->ReactionPlaneAngle();
    }
  //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

  Int_t *motherids;
  motherids = new Int_t[fStack->GetNtrack()];
  for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;

  // Read in mother ids
  TParticle *motherpart;
  for (int ip=0; ip<fStack->GetNtrack(); ip++) {
    motherpart = fStack->Particle(ip);
    if (motherpart->GetDaughter(0) > 0)
      motherids[motherpart->GetDaughter(0)] = ip;
    if (motherpart->GetDaughter(1) > 0)
      motherids[motherpart->GetDaughter(1)] = ip;

//     if (motherpart->GetPdgCode() == 211) {
//       cout << "Mother " << ip << " has daughters " 
// 	   << motherpart->GetDaughter(0) << " " 
// 	   << motherpart->GetDaughter(1) << " " 
// 	   << motherpart->Vx() << " " 
// 	   << motherpart->Vy() << " " 
// 	   << motherpart->Vz() << " " 
// 	   << 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];
      double rxyz[3];
      double impact[2];
      double covimpact[3];
      
      if (fUseTPCOnly) {
	if (!esdtrack->GetTPCInnerParam()) {
	  delete trackCopy;
	  continue;
	}

	AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
	param->GetXYZ(rxyz);
	param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
	param->GetPxPyPz(pxyz);//reading noconstarined momentum

	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	if (v.mag() < 0.0001) {
	  //	cout << "Found 0 momentum ???? " <<endl;
	  delete trackCopy;
	  continue;
	}
	trackCopy->SetP(v);//setting momentum
	trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));

	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	//setting helix I do not if it is ok
	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
	trackCopy->SetHelix(helix);

	//some stuff which could be useful 
	trackCopy->SetImpactD(impact[0]);
	trackCopy->SetImpactZ(impact[1]);
	trackCopy->SetCdd(covimpact[0]);
	trackCopy->SetCdz(covimpact[1]);
	trackCopy->SetCzz(covimpact[2]);
	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));	

	delete param;
      }
      else {
	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]);
	if (v.mag() < 0.0001) {
	  //	cout << "Found 0 momentum ???? " <<endl;
	  delete trackCopy;
	  continue;
	}
	trackCopy->SetP(v);//setting momentum
	trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	//setting helix I do not if it is ok
	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
	trackCopy->SetHelix(helix);
	
	//some stuff which could be useful 
	float imp[2];
	float cim[3];
	esdtrack->GetImpactParameters(imp,cim);

	impact[0] = imp[0];
	impact[1] = imp[1];
	covimpact[0] = cim[0];
	covimpact[1] = cim[1];
	covimpact[2] = cim[2];

	trackCopy->SetImpactD(impact[0]);
	trackCopy->SetImpactZ(impact[1]);
	trackCopy->SetCdd(covimpact[0]);
	trackCopy->SetCdz(covimpact[1]);
	trackCopy->SetCzz(covimpact[2]);
	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
      }
	    	
      trackCopy->SetTrackId(esdtrack->GetID());
      trackCopy->SetFlags(esdtrack->GetStatus());
      trackCopy->SetLabel(esdtrack->GetLabel());
		
      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 xtpc[3];
      esdtrack->GetInnerXYZ(xtpc);
      xtpc[2] -= fV1[2];
      trackCopy->SetNominalTPCEntrancePoint(xtpc);

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

      int indexes[3];
      for (int ik=0; ik<3; ik++) {
	indexes[ik] = esdtrack->GetKinkIndex(ik);
      }
      trackCopy->SetKinkIndexes(indexes);

      // Fill the hidden information with the simulated data
      TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));

      // Check the mother information

      // Using the new way of storing the freeze-out information
      // Final state particle is stored twice on the stack
      // one copy (mother) is stored with original freeze-out information
      //   and is not tracked
      // the other one (daughter) is stored with primary vertex position
      //   and is tracked
	
      // Freeze-out coordinates
      double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
      fpx = tPart->Vx() - fV1[0];
      fpy = tPart->Vy() - fV1[1];
      fpz = tPart->Vz() - fV1[2];
      fpt = tPart->T();

      AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
      tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);

      if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 	TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 	// Check if this is the same particle stored twice on the stack
  	if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
	  // It is the same particle
	  // Read in the original freeze-out information
	  // and convert it from to [fm]
 	  fpx = mother->Vx()*1e13;
 	  fpy = mother->Vy()*1e13;
 	  fpz = mother->Vz()*1e13;
 	  fpt = mother->T()*1e13*3e10;
	  
	}
      }

      tInfo->SetPDGPid(tPart->GetPdgCode());
      tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
      Double_t mass2 = (tPart->Energy() *tPart->Energy() -
			tPart->Px()*tPart->Px() -
			tPart->Py()*tPart->Py() -
			tPart->Pz()*tPart->Pz());
      if (mass2>0.0)
	tInfo->SetMass(TMath::Sqrt(mass2));
      else 
	tInfo->SetMass(0.0);

      tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
      trackCopy->SetHiddenInfo(tInfo);
      
      //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
	  //	  delete trackCopy;
	}
      else
	{
	  delete  trackCopy;
	}
		
    }

  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event	
  fCurEvent++;	
  fCurRLEvent++;
  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
  if (fCurEvent== fNumberofEvent)//if end of current file close all
    {   
      fTree->Reset(); 
      delete fTree;
    }
  return hbtEvent; 
}
//____________________________________________________________________
Float_t AliFemtoEventReaderESDKine::GetSigmaToVertex(const AliESDtrack* esdTrack)
{
  // Calculates the number of sigma to the vertex.

  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];

  b[0] = impact[0];
  b[1] = impact[1];
  bCov[0] = covar[0];
  bCov[1] = covar[1];
  bCov[2] = covar[2];

  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);

  // -----------------------------------
  // How to get to a n-sigma cut?
  //
  // The accumulated statistics from 0 to d is
  //
  // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
  // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
  //
  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
  // Can this be expressed in a different way?

  if (bRes[0] == 0 || bRes[1] ==0)
    return -1;

  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));

  // stupid rounding problem screws up everything:
  // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
  if (TMath::Exp(-d * d / 2) < 1e-10)
    return 1000;

  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  return d;
}
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoEventReaderESDKine - 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
	11	/*
	12	*$Id$
	13	*$Log$
	14	*Revision 1.1 2007/05/25 12:42:54 akisiel
	15	*Adding a reader for the Kine information
	16	*
	17	*Revision 1.2 2007/05/22 09:01:42 akisiel
	18	*Add the possibiloity to save cut settings in the ROOT file
	19	*
	20	*Revision 1.1 2007/05/16 10:22:11 akisiel
	21	*Making the directory structure of AliFemto flat. All files go into one common directory
	22	*
	23	*Revision 1.5 2007/05/03 09:45:20 akisiel
	24	*Fixing Effective C++ warnings
	25	*
	26	*Revision 1.4 2007/04/27 07:28:34 akisiel
	27	*Remove event number reading due to interface changes
	28	*
	29	*Revision 1.3 2007/04/27 07:25:16 akisiel
	30	*Make revisions needed for compilation from the main AliRoot tree
	31	*
	32	*Revision 1.1.1.1 2007/04/25 15:38:41 panos
	33	*Importing the HBT code dir
	34	*
	35	*/
	36
	37	#include "AliFemtoEventReaderESDKine.h"
	38
	39	#include "TFile.h"
	40	#include "TChain.h"
	41	#include "AliESDEvent.h"
	42	#include "AliESDtrack.h"
	43	#include "AliESDVertex.h"
	44	#include "AliStack.h"
	45	//#include "AliAODParticle.h"
	46	#include "TParticle.h"
	47
	48	//#include "TSystem.h"
	49
	50	#include "AliFmPhysicalHelixD.h"
	51	#include "AliFmThreeVectorF.h"
	52
	53	#include "SystemOfUnits.h"
	54
	55	#include "AliFemtoEvent.h"
	56
	57	#include "TMath.h"
	58	#include "TParticle.h"
	59	#include "AliFemtoModelHiddenInfo.h"
	60	#include "AliFemtoModelGlobalHiddenInfo.h"
	61	#include "AliGenHijingEventHeader.h"
	62
	63	ClassImp(AliFemtoEventReaderESDKine)
	64
65	#if !(ST_NO_NAMESPACES)
66	using namespace units;
67	#endif
68
69	using namespace std;
70	//____________________________
71	//constructor with 0 parameters , look at default settings
72	AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
73	fInputFile(" "),
74	fFileName(" "),
75	fConstrained(true),
76	fNumberofEvent(0),
77	fCurEvent(0),
78	fCurRLEvent(0),
79	fTree(0x0),
80	fEvent(0x0),
81	fRunLoader(0x0)
82	{
83	}
84
85	AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
86	AliFemtoEventReader(),
87	fInputFile(" "),
88	fFileName(" "),
89	fConstrained(true),
90	fNumberofEvent(0),
91	fCurEvent(0),
92	fCurRLEvent(0),
93	fTree(0x0),
94	fEvent(0x0),
95	fRunLoader(0x0)
96	{
97	// copy constructor
98	fInputFile = aReader.fInputFile;
99	fFileName = aReader.fFileName;
100	fConstrained = aReader.fConstrained;
101	fNumberofEvent = aReader.fNumberofEvent;
102	fCurEvent = aReader.fCurEvent;
103	fEvent = new AliESDEvent();
104	}
105	//__________________
106	//Destructor
107	AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
108	{
109	// destructor
110	//delete fListOfFiles;
111	delete fTree;
112	delete fEvent;
113	if (fRunLoader) delete fRunLoader;
114	}
115
116	//__________________
117	AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
118	{
119	// assignment operator
120	if (this == &aReader)
121	return *this;
122
123	fInputFile = aReader.fInputFile;
124	fFileName = aReader.fFileName;
125	fConstrained = aReader.fConstrained;
126	fNumberofEvent = aReader.fNumberofEvent;
127	fCurEvent = aReader.fCurEvent;
128	fCurRLEvent = aReader.fCurRLEvent;
129	if (fTree) delete fTree;
130	// fTree = aReader.fTree->CloneTree();
131	if (fEvent) delete fEvent;
132	fEvent = new AliESDEvent();
133	if (fRunLoader) delete fRunLoader;
134	fRunLoader = new AliRunLoader(*aReader.fRunLoader);
135
136	return *this;
137	}
138	//__________________
139	AliFemtoString AliFemtoEventReaderESDKine::Report()
140	{
141	// create reader report
142	AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
143	return temp;
144	}
145
146	//__________________
147	void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
148	{
149	//setting the name of file where names of ESD file are written
150	//it takes only this files which have good trees
151	char buffer[256];
152	fInputFile=string(inputFile);
153	cout<<"Input File set on "<<fInputFile<<endl;
154	ifstream infile(inputFile);
155
156	fTree = new TChain("esdTree");
157
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;
172	fTree->AddFile(buffer);
173	delete tree;
174	}
175	esdFile->Close();
176	}
177	delete esdFile;
178	}
179	}
180	}
181
182	void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
183	{
184	fConstrained=constrained;
185	}
186
187	bool AliFemtoEventReaderESDKine::GetConstrained() const
188	{
189	return fConstrained;
190	}
191
192	AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
193	{
194	// read in a next hbt event from the chain
195	// convert it to AliFemtoEvent and return
196	// for further analysis
197	AliFemtoEvent *hbtEvent = 0;
198	TString tGAliceFilename;
199
200	if (fCurEvent==fNumberofEvent)//open next file
201	{
202	if (fNumberofEvent == 0) {
203	fEvent=new AliESDEvent();
204
205	//ESD data
206	// fEsdFile=TFile::Open(fFileName.c_str(),"READ");
207	// fTree = (TTree*) fEsdFile->Get("esdTree");
208
209	fTree->SetBranchStatus("MuonTracks*",0);
210	fTree->SetBranchStatus("PmdTracks*",0);
211	fTree->SetBranchStatus("TrdTracks*",0);
212	fTree->SetBranchStatus("V0s*",0);
213	fTree->SetBranchStatus("Cascades*",0);
214	fTree->SetBranchStatus("Kinks*",0);
215	fTree->SetBranchStatus("CaloClusters*",0);
216	fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
217	fTree->SetBranchStatus("ESDfriend*",0);
218	fEvent->ReadFromTree(fTree);
219
220	// chain->SetBranchStatus("*",0);
221	// chain->SetBranchStatus("fUniqueID",1);
222	// chain->SetBranchStatus("fTracks",1);
223	// chain->SetBranchStatus("fTracks.*",1);
224	// chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
225	// fTree->SetBranchStatus("fTracks.fCalibContainer",0);
226
227
228	fNumberofEvent=fTree->GetEntries();
229
230	if (fNumberofEvent == 0) {
231	cout<<"no event in input "<<endl;
232	fReaderStatus=1;
233	return hbtEvent;
234	}
235
236	cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
237	fCurEvent=0;
238	// simulation data reading setup
239
240	}
241	else //no more data to read
242	{
243	cout<<"no more files "<<hbtEvent<<endl;
244	fReaderStatus=1;
245	return hbtEvent;
246	}
247	}
248	cout<<"starting to read event "<<fCurEvent<<endl;
249	fTree->GetEvent(fCurEvent);//getting next event
250	// vector<int> tLabelTable;//to check labels
251
252	cout << "fFileName is " << fFileName.Data() << endl;
253	cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
254	if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
255	fFileName = fTree->GetCurrentFile()->GetName();
256	tGAliceFilename = fFileName;
257	tGAliceFilename.ReplaceAll("AliESDs","galice");
258	cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
259	if (fRunLoader) delete fRunLoader;
260	fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
261	if (fRunLoader==0x0)
262	{
263	cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
264	exit(0);
265	}
266	if(fRunLoader->LoadHeader())
267	{
268	cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
269	exit(0);
270	}
271	fRunLoader->LoadKinematics();
272	fCurRLEvent = 0;
273	}
274
275	fRunLoader->GetEvent(fCurRLEvent);
276	AliStack* tStack = 0x0;
277	tStack = fRunLoader->Stack();
278
279	hbtEvent = new AliFemtoEvent;
280	//setting basic things
281	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
282	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
283	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
284	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
285	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
286	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
287	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
288	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
289	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
290	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
291	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
292	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
293
294	//Vertex
295	double fV1[3];
296	double fVCov[6];
297	if (fUseTPCOnly) {
298	fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
299	fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
300	if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
301	fVCov[4] = -1001.0;
302	}
303	else {
304	fEvent->GetPrimaryVertex()->GetXYZ(fV1);
305	fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
306	if (!fEvent->GetPrimaryVertex()->GetStatus())
307	fVCov[4] = -1001.0;
308	}
309
310	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
311	hbtEvent->SetPrimVertPos(vertex);
312	hbtEvent->SetPrimVertCov(fVCov);
313
314	AliGenHijingEventHeader hdh = dynamic_cast<AliGenHijingEventHeader > (fGenHeader);
315
316	Double_t tReactionPlane = 0;
317	if (hdh)
318	{
319	tReactionPlane = hdh->ReactionPlaneAngle();
320	}
321	//starting to reading tracks
322	int nofTracks=0; //number of reconstructed tracks in event
323	nofTracks=fEvent->GetNumberOfTracks();
324	int realnofTracks=0;//number of track which we use ina analysis
325
326	Int_t *motherids;
327	motherids = new Int_t[fStack->GetNtrack()];
328	for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
329
330	// Read in mother ids
331	TParticle *motherpart;
332	for (int ip=0; ip<fStack->GetNtrack(); ip++) {
333	motherpart = fStack->Particle(ip);
334	if (motherpart->GetDaughter(0) > 0)
335	motherids[motherpart->GetDaughter(0)] = ip;
336	if (motherpart->GetDaughter(1) > 0)
337	motherids[motherpart->GetDaughter(1)] = ip;
338
339	// if (motherpart->GetPdgCode() == 211) {
340	// cout << "Mother " << ip << " has daughters "
341	// << motherpart->GetDaughter(0) << " "
342	// << motherpart->GetDaughter(1) << " "
343	// << motherpart->Vx() << " "
344	// << motherpart->Vy() << " "
345	// << motherpart->Vz() << " "
346	// << endl;
347
348	// }
349	}
350
351	for (int i=0;i<nofTracks;i++)
352	{
353	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
354
355	AliFemtoTrack* trackCopy = new AliFemtoTrack();
356	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
357	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
358
359	trackCopy->SetCharge((short)esdtrack->GetSign());
360
361	//in aliroot we have AliPID
362	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
363	//we use only 5 first
364	double esdpid[5];
365	esdtrack->GetESDpid(esdpid);
366	trackCopy->SetPidProbElectron(esdpid[0]);
367	trackCopy->SetPidProbMuon(esdpid[1]);
368	trackCopy->SetPidProbPion(esdpid[2]);
369	trackCopy->SetPidProbKaon(esdpid[3]);
370	trackCopy->SetPidProbProton(esdpid[4]);
371
372	double pxyz[3];
373	double rxyz[3];
374	double impact[2];
375	double covimpact[3];
376
377	if (fUseTPCOnly) {
378	if (!esdtrack->GetTPCInnerParam()) {
379	delete trackCopy;
380	continue;
381	}
382
383	AliExternalTrackParam param = new AliExternalTrackParam(esdtrack->GetTPCInnerParam());
384	param->GetXYZ(rxyz);
385	param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
386	param->GetPxPyPz(pxyz);//reading noconstarined momentum
387
388	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
389	if (v.mag() < 0.0001) {
390	// cout << "Found 0 momentum ???? " <<endl;
391	delete trackCopy;
392	continue;
393	}
394	trackCopy->SetP(v);//setting momentum
395	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
396
397	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
398	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
399	//setting helix I do not if it is ok
400	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
401	trackCopy->SetHelix(helix);
402
403	//some stuff which could be useful
404	trackCopy->SetImpactD(impact[0]);
405	trackCopy->SetImpactZ(impact[1]);
406	trackCopy->SetCdd(covimpact[0]);
407	trackCopy->SetCdz(covimpact[1]);
408	trackCopy->SetCzz(covimpact[2]);
409	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
410
411	delete param;
412	}
413	else {
414	if (fConstrained==true)
415	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
416	else
417	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
418
419	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
420	if (v.mag() < 0.0001) {
421	// cout << "Found 0 momentum ???? " <<endl;
422	delete trackCopy;
423	continue;
424	}
425	trackCopy->SetP(v);//setting momentum
426	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
427	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
428	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
429	//setting helix I do not if it is ok
430	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
431	trackCopy->SetHelix(helix);
432
433	//some stuff which could be useful
434	float imp[2];
435	float cim[3];
436	esdtrack->GetImpactParameters(imp,cim);
437
438	impact[0] = imp[0];
439	impact[1] = imp[1];
440	covimpact[0] = cim[0];
441	covimpact[1] = cim[1];
442	covimpact[2] = cim[2];
443
444	trackCopy->SetImpactD(impact[0]);
445	trackCopy->SetImpactZ(impact[1]);
446	trackCopy->SetCdd(covimpact[0]);
447	trackCopy->SetCdz(covimpact[1]);
448	trackCopy->SetCzz(covimpact[2]);
449	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
450	}
451
452	trackCopy->SetTrackId(esdtrack->GetID());
453	trackCopy->SetFlags(esdtrack->GetStatus());
454	trackCopy->SetLabel(esdtrack->GetLabel());
455
456	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
457	trackCopy->SetITSncls(esdtrack->GetNcls(0));
458	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
459	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
460	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
461	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
462	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
463
464
465	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
466	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
467
468	double xtpc[3];
469	esdtrack->GetInnerXYZ(xtpc);
470	xtpc[2] -= fV1[2];
471	trackCopy->SetNominalTPCEntrancePoint(xtpc);
472
473	esdtrack->GetOuterXYZ(xtpc);
474	xtpc[2] -= fV1[2];
475	trackCopy->SetNominalTPCExitPoint(xtpc);
476
477	int indexes[3];
478	for (int ik=0; ik<3; ik++) {
479	indexes[ik] = esdtrack->GetKinkIndex(ik);
480	}
481	trackCopy->SetKinkIndexes(indexes);
482
483	// Fill the hidden information with the simulated data
484	TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
485
486	// Check the mother information
487
488	// Using the new way of storing the freeze-out information
489	// Final state particle is stored twice on the stack
490	// one copy (mother) is stored with original freeze-out information
491	// and is not tracked
492	// the other one (daughter) is stored with primary vertex position
493	// and is tracked
494
495	// Freeze-out coordinates
496	double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
497	fpx = tPart->Vx() - fV1[0];
498	fpy = tPart->Vy() - fV1[1];
499	fpz = tPart->Vz() - fV1[2];
500	fpt = tPart->T();
501
502	AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
503	tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
504
505	if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
506	TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
507	// Check if this is the same particle stored twice on the stack
508	if ((mother->GetPdgCode() == tPart->GetPdgCode() \|\| (mother->Px() == tPart->Px()))) {
509	// It is the same particle
510	// Read in the original freeze-out information
511	// and convert it from to [fm]
512	fpx = mother->Vx()*1e13;
513	fpy = mother->Vy()*1e13;
514	fpz = mother->Vz()*1e13;
515	fpt = mother->T()1e133e10;
516
517	}
518	}
519
520	tInfo->SetPDGPid(tPart->GetPdgCode());
521	tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
522	Double_t mass2 = (tPart->Energy() *tPart->Energy() -
523	tPart->Px()*tPart->Px() -
524	tPart->Py()*tPart->Py() -
525	tPart->Pz()*tPart->Pz());
526	if (mass2>0.0)
527	tInfo->SetMass(TMath::Sqrt(mass2));
528	else
529	tInfo->SetMass(0.0);
530
531	tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
532	trackCopy->SetHiddenInfo(tInfo);
533
534	//decision if we want this track
535	//if we using diffrent labels we want that this label was use for first time
536	//if we use hidden info we want to have match between sim data and ESD
537	if (tGoodMomentum==true)
538	{
539	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
540	realnofTracks++;//real number of tracks
541	// delete trackCopy;
542	}
543	else
544	{
545	delete trackCopy;
546	}
547
548	}
549
550	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
551	fCurEvent++;
552	fCurRLEvent++;
553	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
554	if (fCurEvent== fNumberofEvent)//if end of current file close all
555	{
556	fTree->Reset();
557	delete fTree;
558	}
559	return hbtEvent;
560	}
561	//____________________________________________________________________
562	Float_t AliFemtoEventReaderESDKine::GetSigmaToVertex(const AliESDtrack* esdTrack)
563	{
564	// Calculates the number of sigma to the vertex.
565
566	Float_t b[2];
567	Float_t bRes[2];
568	Float_t bCov[3];
569
570	b[0] = impact[0];
571	b[1] = impact[1];
572	bCov[0] = covar[0];
573	bCov[1] = covar[1];
574	bCov[2] = covar[2];
575
576	bRes[0] = TMath::Sqrt(bCov[0]);
577	bRes[1] = TMath::Sqrt(bCov[2]);
578
579	// -----------------------------------
580	// How to get to a n-sigma cut?
581	//
582	// The accumulated statistics from 0 to d is
583	//
584	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
585	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
586	//
587	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
588	// Can this be expressed in a different way?
589
590	if (bRes[0] == 0 \|\| bRes[1] ==0)
591	return -1;
592
593	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
594
595	// stupid rounding problem screws up everything:
596	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
597	if (TMath::Exp(-d * d / 2) < 1e-10)
598	return 1000;
599
600	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
601	return d;
602	}