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

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoEventReaderAOD - the reader class for the Alice AOD                //
// Reads in AOD information and converts it into internal AliFemtoEvent       //
// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl                          //
//          Adam Kisiel kisiel@mps.ohio-state.edu                             //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoEventReaderAOD.h"

#include "TFile.h"
#include "TTree.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODVertex.h"
#include "AliAODMCHeader.h"

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

#include "SystemOfUnits.h"

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

ClassImp(AliFemtoEventReaderAOD)

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

using namespace std;
//____________________________
//constructor with 0 parameters , look at default settings 
AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
  fNumberofEvent(0),
  fCurEvent(0),
  fEvent(0x0),
  fAllTrue(160),
  fAllFalse(160),
  fFilterBit(0),
  fPWG2AODTracks(0x0),
  fReadMC(0),
  fInputFile(" "),
  fFileName(" "),
  fTree(0x0),
  fAodFile(0x0)
{
  // default constructor
  fAllTrue.ResetAllBits(kTRUE);
  fAllFalse.ResetAllBits(kFALSE);
}

AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
  AliFemtoEventReader(),
  fNumberofEvent(0),
  fCurEvent(0),
  fEvent(0x0),
  fAllTrue(160),
  fAllFalse(160),
  fFilterBit(0),
  fPWG2AODTracks(0x0),
  fReadMC(0),
  fInputFile(" "),
  fFileName(" "),
  fTree(0x0),
  fAodFile(0x0)
{
  // copy constructor
  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fEvent = new AliAODEvent();
  fAodFile = new TFile(aReader.fAodFile->GetName());
  fAllTrue.ResetAllBits(kTRUE);
  fAllFalse.ResetAllBits(kFALSE);
  fFilterBit = aReader.fFilterBit;
  fPWG2AODTracks = aReader.fPWG2AODTracks;
}
//__________________
//Destructor
AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
{
  // destructor
  delete fTree;
  delete fEvent;
  delete fAodFile;
  if (fPWG2AODTracks) {
    fPWG2AODTracks->Delete();
    delete fPWG2AODTracks;
  }
}

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

  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  if (fTree) delete fTree;
  if (fEvent) delete fEvent;
  fEvent = new AliAODEvent();
  if (fAodFile) delete fAodFile;
  fAodFile = new TFile(aReader.fAodFile->GetName());
  fAllTrue.ResetAllBits(kTRUE);
  fAllFalse.ResetAllBits(kFALSE);
  fFilterBit = aReader.fFilterBit;
  fPWG2AODTracks = aReader.fPWG2AODTracks;

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

//__________________
void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
{
  //setting the name of file where names of AOD file are written 
  //it takes only this files which have good trees
  char buffer[256];
  fInputFile=string(inputFile);
  ifstream infile(inputFile);

  fTree = new TChain("aodTree");

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

AliFemtoEvent* AliFemtoEventReaderAOD::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)	
	{
	  fEvent=new AliAODEvent();
	  fEvent->ReadFromTree(fTree);

	  // Check for the existence of the additional information
	  fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");

	  if (fPWG2AODTracks) {
	    cout << "Found additional PWG2 specific information in the AOD!" << endl;
	    cout << "Reading only tracks with the additional information" << endl;
	  }

	  fNumberofEvent=fTree->GetEntries();
	  //	  cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	  fCurEvent=0;
	}
      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;
	
  hbtEvent = new AliFemtoEvent;

  CopyAODtoFemtoEvent(hbtEvent);

  fCurEvent++;	

  return hbtEvent; 
}

void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
{
  // A function that reads in the AOD event
  // and transfers the neccessary information into
  // the internal AliFemtoEvent

  // setting global event characteristics
  tEvent->SetRunNumber(fEvent->GetRunNumber());
  tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
  tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
  tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
  tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
  tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
  tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
  tEvent->SetZDCParticipants(0);
  tEvent->SetTriggerMask(fEvent->GetTriggerMask());
  tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());

  // Attempt to access MC header
  AliAODMCHeader *mcH;
  TClonesArray *mcP=0;
  if (fReadMC) {
    mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
    if (!mcH) {
      cout << "AOD MC information requested, but no header found!" << endl;
    }

    mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
    if (!mcP) {
      cout << "AOD MC information requested, but no particle array found!" << endl;
    }
  }

  tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);

  Int_t *motherids=0;
  if (mcP) {
    motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
    for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;

    // Read in mother ids
    AliAODMCParticle *motherpart;
    for (int ip=0; ip<mcP->GetEntries(); ip++) {
      motherpart = (AliAODMCParticle *) mcP->At(ip);
      if (motherpart->GetDaughter(0) > 0)
	motherids[motherpart->GetDaughter(0)] = ip;
      if (motherpart->GetDaughter(1) > 0)
	motherids[motherpart->GetDaughter(1)] = ip;
    }
  }

  // Primary Vertex position
  double fV1[3];
  fEvent->GetPrimaryVertex()->GetPosition(fV1);

  AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
  tEvent->SetPrimVertPos(vertex);
	
  //starting to reading tracks
  int nofTracks=0;  //number of reconstructed tracks in event

  // Check to see whether the additional info exists
  if (fPWG2AODTracks)
    nofTracks=fPWG2AODTracks->GetEntries();
  else
    nofTracks=fEvent->GetNumberOfTracks();

  int realnofTracks=0;   // number of track which we use in a analysis

  for (int i=0;i<nofTracks;i++)
    {
      AliFemtoTrack* trackCopy = new AliFemtoTrack();	

      if (fPWG2AODTracks) {
	// Read tracks from the additional pwg2 specific AOD part
	// if they exist
	// Note that in that case all the AOD tracks without the 
	// additional information will be ignored !
	AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);

	// Getting the AOD track through the ref of the additional info
	AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();	
	if (!aodtrack->TestFilterBit(fFilterBit))
	  continue;
	
	CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
	
	if (mcP) {
	  // Fill the hidden information with the simulated data
	  //	  Int_t pLabel = aodtrack->GetLabel();
	  AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->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->Xv() - fV1[0];
	  fpy = tPart->Yv() - fV1[1];
	  fpz = tPart->Zv() - fV1[2];
	  fpt = tPart->T();

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

	  fpx *= 1e13;
	  fpy *= 1e13;
	  fpz *= 1e13;
	  fpt *= 1e13;
	  
	  //      cout << "Looking for mother ids " << endl;
	  if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
	    //	cout << "Got mother id" << endl;
	    AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->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]
	      
	      // EPOS style 
	      // 	  fpx = mother->Xv()*1e13*0.197327;
	      // 	  fpy = mother->Yv()*1e13*0.197327;
	      // 	  fpz = mother->Zv()*1e13*0.197327;
	      // 	  fpt = mother->T() *1e13*0.197327*0.5;
	      
	      
	      // Therminator style 
	      fpx = mother->Xv()*1e13;
	      fpy = mother->Yv()*1e13;
	      fpz = mother->Zv()*1e13;
	      fpt = mother->T() *1e13*3e10;
	      
	    }
	  }
	  
//       if (fRotateToEventPlane) {
// 	double tPhi = TMath::ATan2(fpy, fpx);
// 	double tRad = TMath::Hypot(fpx, fpy);
	
// 	fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
// 	fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
//       }

	  tInfo->SetPDGPid(tPart->GetPdgCode());

// 	  if (fRotateToEventPlane) {
// 	    double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
// 	    double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
	    
// 	    tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
// 				   tRad*TMath::Sin(tPhi - tReactionPlane),
// 				   tPart->Pz());
// 	  }
//       else
	  tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
	  Double_t mass2 = (tPart->E() *tPart->E() -
			    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);

	}

	double pxyz[3];
	aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	// Check the sanity of the tracks - reject zero momentum tracks
	if (ktP.Mag() == 0) {
	  delete trackCopy;
	  continue;
	}
      }
      else {
	// No additional information exists
	// Read in the normal AliAODTracks 
	const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
	
// 	if (!aodtrack->TestFilterBit(fFilterBit))
// 	  continue;
	
	CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
	
	if (mcP) {
	  // Fill the hidden information with the simulated data
	  //	  Int_t pLabel = aodtrack->GetLabel();
	  AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
	  
	  AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
	  double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
	  if (!tPart) {
	    fpx = fV1[0];
	    fpy = fV1[1];
	    fpz = fV1[2];
	    tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
	    tInfo->SetPDGPid(0);
	    tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
	    tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
	    tInfo->SetMass(0);
	  }
	  else {
	    // 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
	    fpx = tPart->Xv() - fV1[0];
	    fpy = tPart->Yv() - fV1[1];
	    fpz = tPart->Zv() - fV1[2];
	    //	  fpt = tPart->T();
	    
	    tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
	    
	    fpx *= 1e13;
	    fpy *= 1e13;
	    fpz *= 1e13;
	    //	  fpt *= 1e13;
	    
	    //      cout << "Looking for mother ids " << endl;
	    if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
	      //	cout << "Got mother id" << endl;
	      AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
	      // Check if this is the same particle stored twice on the stack
	      if (mother) {
		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]
		  
		  // EPOS style 
		  // 	  fpx = mother->Xv()*1e13*0.197327;
		  // 	  fpy = mother->Yv()*1e13*0.197327;
		  // 	  fpz = mother->Zv()*1e13*0.197327;
		  // 	  fpt = mother->T() *1e13*0.197327*0.5;
		  
		  
		  // Therminator style 
		  fpx = mother->Xv()*1e13;
		  fpy = mother->Yv()*1e13;
		  fpz = mother->Zv()*1e13;
		  //	      fpt = mother->T() *1e13*3e10;
		  
		}
	      }
	    }
	    
	    //       if (fRotateToEventPlane) {
	    // 	double tPhi = TMath::ATan2(fpy, fpx);
	    // 	double tRad = TMath::Hypot(fpx, fpy);
	    
	    // 	fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
	    // 	fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
	    //       }
	    
	    tInfo->SetPDGPid(tPart->GetPdgCode());
	    
	    // 	  if (fRotateToEventPlane) {
	    // 	    double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
	    // 	    double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
	    
	    // 	    tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
	    // 				   tRad*TMath::Sin(tPhi - tReactionPlane),
	    // 				   tPart->Pz());
	    // 	  }
	    //       else
	    tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
	    Double_t mass2 = (tPart->E() *tPart->E() -
			      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);
	}

	double pxyz[3];
	aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	// Check the sanity of the tracks - reject zero momentum tracks
	if (ktP.Mag() == 0) {
	  delete trackCopy;
	  continue;
	}
      }


      tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
      realnofTracks++;//real number of tracks		
    }

  tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event	

  if (mcP) delete [] motherids;

  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
}

void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack, 
						 AliFemtoTrack *tFemtoTrack, 
						 AliPWG2AODTrack *tPWG2AODTrack)
{
  // Copy the track information from the AOD into the internal AliFemtoTrack
  // If it exists, use the additional information from the PWG2 AOD

  // Primary Vertex position
  double fV1[3];
  fEvent->GetPrimaryVertex()->GetPosition(fV1);

  tFemtoTrack->SetCharge(tAodTrack->Charge());
  
  //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

  // AOD pid has 10 components
  double aodpid[10];
  tAodTrack->GetPID(aodpid);
  tFemtoTrack->SetPidProbElectron(aodpid[0]);
  tFemtoTrack->SetPidProbMuon(aodpid[1]);
  tFemtoTrack->SetPidProbPion(aodpid[2]);
  tFemtoTrack->SetPidProbKaon(aodpid[3]);
  tFemtoTrack->SetPidProbProton(aodpid[4]);
						
  double pxyz[3];
  tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
  AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
  tFemtoTrack->SetP(v);//setting momentum
  tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
  const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
  //setting track helix 
  const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
  AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge())); 
  tFemtoTrack->SetHelix(helix);
	    	
  // Flags
  tFemtoTrack->SetTrackId(tAodTrack->GetID());
  tFemtoTrack->SetFlags(1);
  tFemtoTrack->SetLabel(tAodTrack->GetLabel());
		
  // Track quality information 
  float covmat[6];
  tAodTrack->GetCovMatrix(covmat);
  tFemtoTrack->SetImpactD(covmat[0]);
  tFemtoTrack->SetImpactZ(covmat[2]);
  tFemtoTrack->SetCdd(covmat[3]);
  tFemtoTrack->SetCdz(covmat[4]);
  tFemtoTrack->SetCzz(covmat[5]);
  // This information is only available in the ESD
  // We put in fake values or reasonable estimates
  tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());    
  tFemtoTrack->SetITSncls(1);     
  tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());       
  tFemtoTrack->SetTPCncls(1);       
  tFemtoTrack->SetTPCnclsF(1);      
  tFemtoTrack->SetTPCsignalN(1); 
  tFemtoTrack->SetTPCsignalS(1); 

  if (tPWG2AODTrack) {
    // Copy the PWG2 specific information if it exists
    tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
    tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
    
    double xtpc[3] = {0,0,0};
    tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
    tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
    tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
    tFemtoTrack->SetNominalTPCExitPoint(xtpc);
  }
  else {
    // If not use dummy values
    tFemtoTrack->SetTPCClusterMap(fAllTrue);
    tFemtoTrack->SetTPCSharedMap(fAllFalse);
    
    double xtpc[3] = {0,0,0};
    tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
    tFemtoTrack->SetNominalTPCExitPoint(xtpc);
  }

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

void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
{
  fFilterBit = (1 << (ibit));
}

void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
{
  fReadMC = a;
}

AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
{
  if (aLabel < 0) return 0;
  AliAODMCParticle *aodP;
  Int_t posstack = 0;
  if (aLabel > mcP->GetEntries())
    posstack = mcP->GetEntries();
  else
    posstack = aLabel;

  aodP = (AliAODMCParticle *) mcP->At(posstack);
  if (aodP->GetLabel() > posstack) {
    do {
      aodP = (AliAODMCParticle *) mcP->At(posstack);
      if (aodP->GetLabel() == aLabel) return aodP;
      posstack--;
    }
    while (posstack > 0);
  }
  else {
    do {
      aodP = (AliAODMCParticle *) mcP->At(posstack);
      if (aodP->GetLabel() == aLabel) return aodP;
      posstack++;
    }
    while (posstack < mcP->GetEntries());
  }
  
  return 0;
}
Commit	Line	Data
de568db1	1	////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoEventReaderAOD - the reader class for the Alice AOD //
	4	// Reads in AOD information and converts it into internal AliFemtoEvent //
	5	// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl //
	6	// Adam Kisiel kisiel@mps.ohio-state.edu //
	7	// //
	8	////////////////////////////////////////////////////////////////////////////////
	9
	10	#include "AliFemtoEventReaderAOD.h"
	11
	12	#include "TFile.h"
	13	#include "TTree.h"
	14	#include "AliAODEvent.h"
	15	#include "AliAODTrack.h"
	16	#include "AliAODVertex.h"
683877d2	17	#include "AliAODMCHeader.h"
de568db1	18
	19	#include "AliFmPhysicalHelixD.h"
	20	#include "AliFmThreeVectorF.h"
	21
	22	#include "SystemOfUnits.h"
	23
	24	#include "AliFemtoEvent.h"
	25	#include "AliFemtoModelHiddenInfo.h"
683877d2	26	#include "AliFemtoModelGlobalHiddenInfo.h"
de568db1	27
	28	ClassImp(AliFemtoEventReaderAOD)
	29
	30	#if !(ST_NO_NAMESPACES)
	31	using namespace units;
	32	#endif
	33
	34	using namespace std;
	35	//____________________________
	36	//constructor with 0 parameters , look at default settings
	37	AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
de568db1	38	fNumberofEvent(0),
de568db1	39	fCurEvent(0),
de568db1	40	fEvent(0x0),
	41	fAllTrue(160),
	42	fAllFalse(160),
	43	fFilterBit(0),
fcda1d4e	44	fPWG2AODTracks(0x0),
683877d2	45	fReadMC(0),
fcda1d4e	46	fInputFile(" "),
	47	fFileName(" "),
	48	fTree(0x0),
	49	fAodFile(0x0)
de568db1	50	{
	51	// default constructor
	52	fAllTrue.ResetAllBits(kTRUE);
	53	fAllFalse.ResetAllBits(kFALSE);
	54	}
	55
	56	AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
fcda1d4e	57	AliFemtoEventReader(),
de568db1	58	fNumberofEvent(0),
de568db1	59	fCurEvent(0),
de568db1	60	fEvent(0x0),
	61	fAllTrue(160),
	62	fAllFalse(160),
	63	fFilterBit(0),
fcda1d4e	64	fPWG2AODTracks(0x0),
683877d2	65	fReadMC(0),
fcda1d4e	66	fInputFile(" "),
	67	fFileName(" "),
	68	fTree(0x0),
	69	fAodFile(0x0)
de568db1	70	{
	71	// copy constructor
	72	fInputFile = aReader.fInputFile;
	73	fFileName = aReader.fFileName;
	74	fNumberofEvent = aReader.fNumberofEvent;
	75	fCurEvent = aReader.fCurEvent;
	76	fEvent = new AliAODEvent();
	77	fAodFile = new TFile(aReader.fAodFile->GetName());
	78	fAllTrue.ResetAllBits(kTRUE);
	79	fAllFalse.ResetAllBits(kFALSE);
	80	fFilterBit = aReader.fFilterBit;
	81	fPWG2AODTracks = aReader.fPWG2AODTracks;
	82	}
	83	//__________________
	84	//Destructor
	85	AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
	86	{
	87	// destructor
	88	delete fTree;
	89	delete fEvent;
	90	delete fAodFile;
	91	if (fPWG2AODTracks) {
	92	fPWG2AODTracks->Delete();
	93	delete fPWG2AODTracks;
	94	}
	95	}
	96
	97	//__________________
	98	AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
	99	{
	100	// assignment operator
	101	if (this == &aReader)
	102	return *this;
	103
	104	fInputFile = aReader.fInputFile;
	105	fFileName = aReader.fFileName;
	106	fNumberofEvent = aReader.fNumberofEvent;
	107	fCurEvent = aReader.fCurEvent;
	108	if (fTree) delete fTree;
	109	if (fEvent) delete fEvent;
	110	fEvent = new AliAODEvent();
	111	if (fAodFile) delete fAodFile;
	112	fAodFile = new TFile(aReader.fAodFile->GetName());
	113	fAllTrue.ResetAllBits(kTRUE);
	114	fAllFalse.ResetAllBits(kFALSE);
	115	fFilterBit = aReader.fFilterBit;
	116	fPWG2AODTracks = aReader.fPWG2AODTracks;
	117
	118	return *this;
	119	}
	120	//__________________
	121	AliFemtoString AliFemtoEventReaderAOD::Report()
	122	{
	123	// create reader report
	124	AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
	125	return temp;
	126	}
	127
	128	//__________________
	129	void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
	130	{
	131	//setting the name of file where names of AOD file are written
	132	//it takes only this files which have good trees
	133	char buffer[256];
134	fInputFile=string(inputFile);
de568db1	135	ifstream infile(inputFile);
	136
	137	fTree = new TChain("aodTree");
	138
	139	if(infile.good()==true)
	140	{
	141	//checking if all give files have good tree inside
	142	while (infile.eof()==false)
	143	{
	144	infile.getline(buffer,256);
	145	TFile *aodFile=TFile::Open(buffer,"READ");
	146	if (aodFile!=0x0)
	147	{
	148	TTree* tree = (TTree*) aodFile->Get("aodTree");
	149	if (tree!=0x0)
	150	{
683877d2	151	// cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
de568db1	152	fTree->AddFile(buffer);
	153	delete tree;
	154	}
	155	aodFile->Close();
	156	}
	157	delete aodFile;
	158	}
	159	}
	160	}
	161
	162	AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
	163	{
	164	// read in a next hbt event from the chain
	165	// convert it to AliFemtoEvent and return
	166	// for further analysis
	167	AliFemtoEvent *hbtEvent = 0;
	168
	169	if (fCurEvent==fNumberofEvent)//open next file
	170	{
	171	if(fNumberofEvent==0)
	172	{
	173	fEvent=new AliAODEvent();
	174	fEvent->ReadFromTree(fTree);
	175
	176	// Check for the existence of the additional information
	177	fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
	178
	179	if (fPWG2AODTracks) {
	180	cout << "Found additional PWG2 specific information in the AOD!" << endl;
	181	cout << "Reading only tracks with the additional information" << endl;
	182	}
	183
	184	fNumberofEvent=fTree->GetEntries();
683877d2	185	// cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
de568db1	186	fCurEvent=0;
	187	}
	188	else //no more data to read
	189	{
	190	cout<<"no more files "<<hbtEvent<<endl;
	191	fReaderStatus=1;
	192	return hbtEvent;
	193	}
	194	}
	195
	196	cout<<"starting to read event "<<fCurEvent<<endl;
	197	fTree->GetEvent(fCurEvent);//getting next event
683877d2	198	// cout << "Read event " << fEvent << " from file " << fTree << endl;
de568db1	199
	200	hbtEvent = new AliFemtoEvent;
	201
	202	CopyAODtoFemtoEvent(hbtEvent);
	203
	204	fCurEvent++;
	205
	206	return hbtEvent;
	207	}
	208
	209	void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
	210	{
	211	// A function that reads in the AOD event
	212	// and transfers the neccessary information into
	213	// the internal AliFemtoEvent
	214
	215	// setting global event characteristics
	216	tEvent->SetRunNumber(fEvent->GetRunNumber());
	217	tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	218	tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	219	tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	220	tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	221	tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	222	tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
b976850d	223	tEvent->SetZDCParticipants(0);
de568db1	224	tEvent->SetTriggerMask(fEvent->GetTriggerMask());
de568db1	225	tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
683877d2	226
	227	// Attempt to access MC header
	228	AliAODMCHeader *mcH;
d9ae2120	229	TClonesArray *mcP=0;
683877d2	230	if (fReadMC) {
	231	mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
	232	if (!mcH) {
	233	cout << "AOD MC information requested, but no header found!" << endl;
	234	}
	235
	236	mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
	237	if (!mcP) {
	238	cout << "AOD MC information requested, but no particle array found!" << endl;
	239	}
	240	}
	241
fee52126	242	tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
fee52126	243
d9ae2120	244	Int_t *motherids=0;
683877d2	245	if (mcP) {
	246	motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
	247	for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
	248
	249	// Read in mother ids
	250	AliAODMCParticle *motherpart;
	251	for (int ip=0; ip<mcP->GetEntries(); ip++) {
	252	motherpart = (AliAODMCParticle *) mcP->At(ip);
	253	if (motherpart->GetDaughter(0) > 0)
	254	motherids[motherpart->GetDaughter(0)] = ip;
	255	if (motherpart->GetDaughter(1) > 0)
	256	motherids[motherpart->GetDaughter(1)] = ip;
	257	}
	258	}
	259
de568db1	260	// Primary Vertex position
	261	double fV1[3];
	262	fEvent->GetPrimaryVertex()->GetPosition(fV1);
	263
	264	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	265	tEvent->SetPrimVertPos(vertex);
	266
	267	//starting to reading tracks
	268	int nofTracks=0; //number of reconstructed tracks in event
	269
	270	// Check to see whether the additional info exists
	271	if (fPWG2AODTracks)
	272	nofTracks=fPWG2AODTracks->GetEntries();
	273	else
	274	nofTracks=fEvent->GetNumberOfTracks();
	275
	276	int realnofTracks=0; // number of track which we use in a analysis
de568db1	277
	278	for (int i=0;i<nofTracks;i++)
	279	{
	280	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	281
	282	if (fPWG2AODTracks) {
	283	// Read tracks from the additional pwg2 specific AOD part
	284	// if they exist
	285	// Note that in that case all the AOD tracks without the
	286	// additional information will be ignored !
	287	AliPWG2AODTrack pwg2aodtrack = (AliPWG2AODTrack ) fPWG2AODTracks->At(i);
	288
	289	// Getting the AOD track through the ref of the additional info
	290	AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
	291	if (!aodtrack->TestFilterBit(fFilterBit))
	292	continue;
	293
	294	CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
	295
683877d2	296	if (mcP) {
683877d2	297	// Fill the hidden information with the simulated data
d9ae2120	298	// Int_t pLabel = aodtrack->GetLabel();
683877d2	299	AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
	300
	301	// Check the mother information
	302
	303	// Using the new way of storing the freeze-out information
	304	// Final state particle is stored twice on the stack
	305	// one copy (mother) is stored with original freeze-out information
	306	// and is not tracked
	307	// the other one (daughter) is stored with primary vertex position
	308	// and is tracked
	309
	310	// Freeze-out coordinates
	311	double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
	312	fpx = tPart->Xv() - fV1[0];
	313	fpy = tPart->Yv() - fV1[1];
	314	fpz = tPart->Zv() - fV1[2];
	315	fpt = tPart->T();
	316
	317	AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
	318	tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
	319
	320	fpx *= 1e13;
	321	fpy *= 1e13;
	322	fpz *= 1e13;
	323	fpt *= 1e13;
	324
	325	// cout << "Looking for mother ids " << endl;
	326	if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
	327	// cout << "Got mother id" << endl;
	328	AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
	329	// Check if this is the same particle stored twice on the stack
	330	if ((mother->GetPdgCode() == tPart->GetPdgCode() \|\| (mother->Px() == tPart->Px()))) {
	331	// It is the same particle
	332	// Read in the original freeze-out information
	333	// and convert it from to [fm]
	334
	335	// EPOS style
	336	// fpx = mother->Xv()1e130.197327;
	337	// fpy = mother->Yv()1e130.197327;
	338	// fpz = mother->Zv()1e130.197327;
	339	// fpt = mother->T() 1e130.197327*0.5;
	340
	341
	342	// Therminator style
	343	fpx = mother->Xv()*1e13;
	344	fpy = mother->Yv()*1e13;
	345	fpz = mother->Zv()*1e13;
	346	fpt = mother->T() 1e133e10;
	347
	348	}
	349	}
	350
	351	// if (fRotateToEventPlane) {
	352	// double tPhi = TMath::ATan2(fpy, fpx);
	353	// double tRad = TMath::Hypot(fpx, fpy);
	354
	355	// fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
	356	// fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
	357	// }
	358
	359	tInfo->SetPDGPid(tPart->GetPdgCode());
	360
	361	// if (fRotateToEventPlane) {
	362	// double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
363	// double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
364
365	// tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
366	// tRad*TMath::Sin(tPhi - tReactionPlane),
367	// tPart->Pz());
368	// }
369	// else
370	tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
371	Double_t mass2 = (tPart->E() *tPart->E() -
372	tPart->Px()*tPart->Px() -
373	tPart->Py()*tPart->Py() -
374	tPart->Pz()*tPart->Pz());
375	if (mass2>0.0)
376	tInfo->SetMass(TMath::Sqrt(mass2));
377	else
378	tInfo->SetMass(0.0);
379
380	tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
381	trackCopy->SetHiddenInfo(tInfo);
382
383	}
384
de568db1	385	double pxyz[3];
	386	aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
	387	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	388	// Check the sanity of the tracks - reject zero momentum tracks
69c1c8ff	389	if (ktP.Mag() == 0) {
de568db1	390	delete trackCopy;
	391	continue;
	392	}
	393	}
	394	else {
	395	// No additional information exists
	396	// Read in the normal AliAODTracks
	397	const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
	398
683877d2	399	// if (!aodtrack->TestFilterBit(fFilterBit))
683877d2	400	// continue;
de568db1	401
	402	CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
	403
683877d2	404	if (mcP) {
683877d2	405	// Fill the hidden information with the simulated data
d9ae2120	406	// Int_t pLabel = aodtrack->GetLabel();
683877d2	407	AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
	408
	409	AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
	410	double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
	411	if (!tPart) {
	412	fpx = fV1[0];
	413	fpy = fV1[1];
	414	fpz = fV1[2];
	415	tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
	416	tInfo->SetPDGPid(0);
	417	tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
	418	tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
	419	tInfo->SetMass(0);
	420	}
	421	else {
	422	// Check the mother information
	423
	424	// Using the new way of storing the freeze-out information
	425	// Final state particle is stored twice on the stack
	426	// one copy (mother) is stored with original freeze-out information
	427	// and is not tracked
	428	// the other one (daughter) is stored with primary vertex position
	429	// and is tracked
	430
	431	// Freeze-out coordinates
	432	fpx = tPart->Xv() - fV1[0];
	433	fpy = tPart->Yv() - fV1[1];
	434	fpz = tPart->Zv() - fV1[2];
	435	// fpt = tPart->T();
	436
	437	tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
	438
	439	fpx *= 1e13;
	440	fpy *= 1e13;
	441	fpz *= 1e13;
	442	// fpt *= 1e13;
	443
	444	// cout << "Looking for mother ids " << endl;
	445	if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
	446	// cout << "Got mother id" << endl;
	447	AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
	448	// Check if this is the same particle stored twice on the stack
	449	if (mother) {
	450	if ((mother->GetPdgCode() == tPart->GetPdgCode() \|\| (mother->Px() == tPart->Px()))) {
	451	// It is the same particle
	452	// Read in the original freeze-out information
	453	// and convert it from to [fm]
	454
	455	// EPOS style
	456	// fpx = mother->Xv()1e130.197327;
	457	// fpy = mother->Yv()1e130.197327;
	458	// fpz = mother->Zv()1e130.197327;
	459	// fpt = mother->T() 1e130.197327*0.5;
	460
	461
	462	// Therminator style
	463	fpx = mother->Xv()*1e13;
	464	fpy = mother->Yv()*1e13;
	465	fpz = mother->Zv()*1e13;
	466	// fpt = mother->T() 1e133e10;
	467
	468	}
	469	}
	470	}
471
472	// if (fRotateToEventPlane) {
473	// double tPhi = TMath::ATan2(fpy, fpx);
474	// double tRad = TMath::Hypot(fpx, fpy);
475
476	// fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
477	// fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
478	// }
479
480	tInfo->SetPDGPid(tPart->GetPdgCode());
481
482	// if (fRotateToEventPlane) {
483	// double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
484	// double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
485
486	// tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
487	// tRad*TMath::Sin(tPhi - tReactionPlane),
488	// tPart->Pz());
489	// }
490	// else
491	tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
492	Double_t mass2 = (tPart->E() *tPart->E() -
493	tPart->Px()*tPart->Px() -
494	tPart->Py()*tPart->Py() -
495	tPart->Pz()*tPart->Pz());
496	if (mass2>0.0)
497	tInfo->SetMass(TMath::Sqrt(mass2));
498	else
499	tInfo->SetMass(0.0);
500
501	tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
502	}
503	trackCopy->SetHiddenInfo(tInfo);
504	}
505
de568db1	506	double pxyz[3];
	507	aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
	508	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	509	// Check the sanity of the tracks - reject zero momentum tracks
69c1c8ff	510	if (ktP.Mag() == 0) {
de568db1	511	delete trackCopy;
	512	continue;
	513	}
	514	}
	515
	516
	517	tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	518	realnofTracks++;//real number of tracks
	519	}
	520
	521	tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
	522
4ce766eb	523	if (mcP) delete [] motherids;
683877d2	524
de568db1	525	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
	526	}
	527
	528	void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack,
	529	AliFemtoTrack *tFemtoTrack,
	530	AliPWG2AODTrack *tPWG2AODTrack)
	531	{
	532	// Copy the track information from the AOD into the internal AliFemtoTrack
	533	// If it exists, use the additional information from the PWG2 AOD
	534
	535	// Primary Vertex position
	536	double fV1[3];
	537	fEvent->GetPrimaryVertex()->GetPosition(fV1);
	538
	539	tFemtoTrack->SetCharge(tAodTrack->Charge());
	540
	541	//in aliroot we have AliPID
	542	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	543	//we use only 5 first
	544
	545	// AOD pid has 10 components
	546	double aodpid[10];
	547	tAodTrack->GetPID(aodpid);
	548	tFemtoTrack->SetPidProbElectron(aodpid[0]);
	549	tFemtoTrack->SetPidProbMuon(aodpid[1]);
	550	tFemtoTrack->SetPidProbPion(aodpid[2]);
	551	tFemtoTrack->SetPidProbKaon(aodpid[3]);
	552	tFemtoTrack->SetPidProbProton(aodpid[4]);
	553
	554	double pxyz[3];
	555	tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
	556	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	557	tFemtoTrack->SetP(v);//setting momentum
	558	tFemtoTrack->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	559	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	560	//setting track helix
	561	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	562	AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
	563	tFemtoTrack->SetHelix(helix);
	564
	565	// Flags
	566	tFemtoTrack->SetTrackId(tAodTrack->GetID());
	567	tFemtoTrack->SetFlags(1);
	568	tFemtoTrack->SetLabel(tAodTrack->GetLabel());
	569
	570	// Track quality information
	571	float covmat[6];
	572	tAodTrack->GetCovMatrix(covmat);
	573	tFemtoTrack->SetImpactD(covmat[0]);
	574	tFemtoTrack->SetImpactZ(covmat[2]);
	575	tFemtoTrack->SetCdd(covmat[3]);
	576	tFemtoTrack->SetCdz(covmat[4]);
	577	tFemtoTrack->SetCzz(covmat[5]);
	578	// This information is only available in the ESD
	579	// We put in fake values or reasonable estimates
	580	tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
	581	tFemtoTrack->SetITSncls(1);
	582	tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
	583	tFemtoTrack->SetTPCncls(1);
683877d2	584	tFemtoTrack->SetTPCnclsF(1);
	585	tFemtoTrack->SetTPCsignalN(1);
	586	tFemtoTrack->SetTPCsignalS(1);
de568db1	587
	588	if (tPWG2AODTrack) {
	589	// Copy the PWG2 specific information if it exists
	590	tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
	591	tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
	592
	593	double xtpc[3] = {0,0,0};
	594	tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
	595	tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
	596	tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
	597	tFemtoTrack->SetNominalTPCExitPoint(xtpc);
	598	}
	599	else {
	600	// If not use dummy values
	601	tFemtoTrack->SetTPCClusterMap(fAllTrue);
	602	tFemtoTrack->SetTPCSharedMap(fAllFalse);
	603
	604	double xtpc[3] = {0,0,0};
	605	tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
	606	tFemtoTrack->SetNominalTPCExitPoint(xtpc);
	607	}
	608
	609	int indexes[3];
	610	for (int ik=0; ik<3; ik++) {
	611	indexes[ik] = 0;
	612	}
	613	tFemtoTrack->SetKinkIndexes(indexes);
	614	}
	615
	616	void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
	617	{
	618	fFilterBit = (1 << (ibit));
	619	}
	620
683877d2	621	void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
	622	{
	623	fReadMC = a;
	624	}
	625
	626	AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
	627	{
	628	if (aLabel < 0) return 0;
	629	AliAODMCParticle *aodP;
	630	Int_t posstack = 0;
	631	if (aLabel > mcP->GetEntries())
	632	posstack = mcP->GetEntries();
	633	else
	634	posstack = aLabel;
	635
	636	aodP = (AliAODMCParticle *) mcP->At(posstack);
	637	if (aodP->GetLabel() > posstack) {
	638	do {
	639	aodP = (AliAODMCParticle *) mcP->At(posstack);
	640	if (aodP->GetLabel() == aLabel) return aodP;
	641	posstack--;
	642	}
	643	while (posstack > 0);
	644	}
	645	else {
	646	do {
	647	aodP = (AliAODMCParticle *) mcP->At(posstack);
	648	if (aodP->GetLabel() == aLabel) return aodP;
	649	posstack++;
	650	}
	651	while (posstack < mcP->GetEntries());
	652	}
	653
	654	return 0;
	655	}
de568db1	656
	657
	658
	659
	660