[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
	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"
	17	#include "AliAODMCHeader.h"
	18
	19	#include "AliFmPhysicalHelixD.h"
	20	#include "AliFmThreeVectorF.h"
	21
	22	#include "SystemOfUnits.h"
	23
	24	#include "AliFemtoEvent.h"
	25	#include "AliFemtoModelHiddenInfo.h"
	26	#include "AliFemtoModelGlobalHiddenInfo.h"
	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():
	38	fNumberofEvent(0),
	39	fCurEvent(0),
	40	fEvent(0x0),
	41	fAllTrue(160),
	42	fAllFalse(160),
	43	fFilterBit(0),
	44	fPWG2AODTracks(0x0),
	45	fReadMC(0),
	46	fInputFile(" "),
	47	fFileName(" "),
	48	fTree(0x0),
	49	fAodFile(0x0)
	50	{
	51	// default constructor
	52	fAllTrue.ResetAllBits(kTRUE);
	53	fAllFalse.ResetAllBits(kFALSE);
	54	}
	55
	56	AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
	57	AliFemtoEventReader(),
	58	fNumberofEvent(0),
	59	fCurEvent(0),
	60	fEvent(0x0),
	61	fAllTrue(160),
	62	fAllFalse(160),
	63	fFilterBit(0),
	64	fPWG2AODTracks(0x0),
	65	fReadMC(0),
	66	fInputFile(" "),
	67	fFileName(" "),
	68	fTree(0x0),
	69	fAodFile(0x0)
	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);
	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	{
	151	// cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
	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();
	185	// cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	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
	198	// cout << "Read event " << fEvent << " from file " << fTree << endl;
	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));
	223	tEvent->SetZDCParticipants(0);
	224	tEvent->SetTriggerMask(fEvent->GetTriggerMask());
	225	tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	226
	227	// Attempt to access MC header
	228	AliAODMCHeader *mcH;
	229	TClonesArray *mcP=0;
	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
	242	tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
	243
	244	Int_t *motherids=0;
	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
	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
	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
	296	if (mcP) {
	297	// Fill the hidden information with the simulated data
	298	// Int_t pLabel = aodtrack->GetLabel();
	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
	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
	389	if (ktP.Mag() == 0) {
	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
	399	// if (!aodtrack->TestFilterBit(fFilterBit))
	400	// continue;
	401
	402	CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
	403
	404	if (mcP) {
	405	// Fill the hidden information with the simulated data
	406	// Int_t pLabel = aodtrack->GetLabel();
	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
	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
	510	if (ktP.Mag() == 0) {
	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
	523	if (mcP) delete [] motherids;
	524
	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);
	584	tFemtoTrack->SetTPCnclsF(1);
	585	tFemtoTrack->SetTPCsignalN(1);
	586	tFemtoTrack->SetTPCsignalS(1);
	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
	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	}
	656
	657
	658
	659
	660