[u/mrichter/AliRoot.git] / RICH / AliRICHTracker.cxx

#include "AliRICHTracker.h"
#include <AliESD.h>
#include <TVector3.h>
#include <TTree.h>
#include "AliRICH.h"
#include "AliRICHHelix.h"
#include <AliMagF.h>
#include "AliRICHRecon.h"
#include <AliStack.h>
#include <TParticle.h>
#include <TDatabasePDG.h>
#include <TMath.h>
ClassImp(AliRICHTracker)
//__________________________________________________________________________________________________
Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
{
  //invoked by AliRecontruction for RICH
  //if ESD doesn't contain tracks, try to reconstruct with particle from STACK 
  //(this case is just to forsee a standalone RICH simulation
  
  AliDebug(1,"Start pattern recognition");
  if(pESD->GetNumberOfTracks()) RecWithESD(pESD);
  else
    RecWithStack();
  AliDebug(1,"Stop pattern recognition");

  return 0; // error code: 0=no error;
} //pure virtual from AliTracker
//__________________________________________________________________________________________________
void AliRICHTracker::RecWithESD(AliESD *pESD)
{
  //recontruction from ESD
  //
  Int_t iNtracks=pESD->GetNumberOfTracks();
  Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
  AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b));
  Double_t xb[3],pb[3];//tmp storage for track parameters
  TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
  
  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
  
  for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
    AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
//  if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
    pTrack->GetXYZ(xb); 
      pTrack->GetPxPyPz(pb); 
          Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters
    AliDebug(1,Form("Track %i pmod=%f mass=%f stat=%i",iTrackN,pTrack->GetP(),pTrack->GetMass(),status));
    x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]);
    AliRICHHelix helix(x0,p0,pTrack->GetSign(),b);   
    Int_t iChamber=helix.RichIntersect(pRich->P());        
    AliDebug(1,Form("intersection with %i chamber found",iChamber));
    if(!iChamber) continue;//intersection with no chamber found
    
    Double_t distMip=9999; //min distance between clusters and track position on PC 
    Int_t iMipId=0; //index of that min distance cluster 
    for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
      AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
      Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
      if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track 
      
      AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
                                                                    helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
    }////clusters loop for intersected chamber
    
    AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
    
    AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
    Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
    AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
    pTrack->SetRICHsignal(thetaCerenkov);
        
    Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
    CalcProb(thetaCerenkov,p0.Mag(),richPID);
    pTrack->SetRICHpid(richPID);         
    
  }//ESD tracks loop
  AliDebug(1,"Stop.");  
} //RecWithESD
//__________________________________________________________________________________________________
void AliRICHTracker::RecWithStack()
{
  // reconstruction for particles from STACK
  //
  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
  
  pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
  AliStack *pStack =   pRich->GetLoader()->GetRunLoader()->Stack();
  Int_t iNtracks=pStack->GetNtrack();

  Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
  AliDebug(1,Form("Start with simulated %i tracks in %f Tesla field",iNtracks,b));
  TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
  

  for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
    TParticle *pParticle = pStack->Particle(iTrackN);
    p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());
    x0.SetXYZ(pParticle->Vx(),pParticle->Vy(),pParticle->Vz());
    AliRICHHelix helix(x0,p0,Int_t(pParticle->GetPDG()->Charge()/3+0.1),b);   
    Int_t iChamber=helix.RichIntersect(pRich->P());        
    AliDebug(1,Form("intersection with %i chamber found",iChamber));
    if(!iChamber) continue;//intersection with no chamber found
    
    Double_t distMip=9999; //min distance between clusters and track position on PC 
    Int_t iMipId=0; //index of that min distance cluster 
    for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
      AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
      Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
      if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track 
      
      AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
                                                                    helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
    }////clusters loop for intersected chamber
    
    AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
    
    AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
    Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
    AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
//    pTrack->SetRICHsignal(thetaCerenkov);

    Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
    CalcProb(thetaCerenkov,p0.Mag(),richPID);
    
  }//ESD tracks loop
  AliDebug(1,"Stop.");  
} //RecWithStack

Int_t AliRICHTracker::LoadClusters(TTree *pTree)
{
// Load clusters for RICH
  AliDebug(1,"Start.");  pTree->GetEntry(0);  AliDebug(1,"Stop."); return 0;
}

//__________________________________________________________________________________________________
void AliRICHTracker::CalcProb(Double_t thetaCer,Double_t pmod, Double_t *richPID)
{
// 
  Double_t height[5];Double_t totalHeight=0;
  Int_t code[5]={kElectron,kMuonPlus,kPiPlus,kKPlus,kProton};
  TDatabasePDG *db = new TDatabasePDG();
  for(Int_t iPart=0;iPart<4;iPart++){
    Double_t mass = db->GetParticle(code[iPart])->Mass();
    Double_t refIndex=AliRICHParam::IndOfRefC6F14(6.755);
    Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(refIndex*pmod);
    if(cosThetaTh>=1) {break;}
    Double_t thetaTh = TMath::ACos(cosThetaTh);
    Double_t sinThetaThNorm = TMath::Sin(thetaTh)/TMath::Sqrt(1-1/(refIndex*refIndex));
    Double_t sigmaThetaTh = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
    height[iPart] = TMath::Gaus(thetaCer,thetaTh,sigmaThetaTh);
    totalHeight +=height[iPart];
  }
  for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;    
}//CalcProb
Commit	Line	Data
998b831f	1	#include "AliRICHTracker.h"
	2	#include <AliESD.h>
	3	#include <TVector3.h>
	4	#include <TTree.h>
	5	#include "AliRICH.h"
	6	#include "AliRICHHelix.h"
	7	#include <AliMagF.h>
	8	#include "AliRICHRecon.h"
	9	#include <AliStack.h>
	10	#include <TParticle.h>
30c60010	11	#include <TDatabasePDG.h>
30c60010	12	#include <TMath.h>
998b831f	13	ClassImp(AliRICHTracker)
30c60010	14	//__________________________________________________________________________________________________
998b831f	15	Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
998b831f	16	{
30c60010	17	//invoked by AliRecontruction for RICH
	18	//if ESD doesn't contain tracks, try to reconstruct with particle from STACK
	19	//(this case is just to forsee a standalone RICH simulation
	20
	21	AliDebug(1,"Start pattern recognition");
	22	if(pESD->GetNumberOfTracks()) RecWithESD(pESD);
	23	else
	24	RecWithStack();
	25	AliDebug(1,"Stop pattern recognition");
	26
	27	return 0; // error code: 0=no error;
	28	} //pure virtual from AliTracker
	29	//__________________________________________________________________________________________________
	30	void AliRICHTracker::RecWithESD(AliESD *pESD)
	31	{
	32	//recontruction from ESD
	33	//
998b831f	34	Int_t iNtracks=pESD->GetNumberOfTracks();
	35	Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
	36	AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b));
	37	Double_t xb[3],pb[3];//tmp storage for track parameters
	38	TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
	39
	40	AliRICH pRich=((AliRICH)gAlice->GetDetector("RICH"));
	41
998b831f	42	for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
	43	AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
	44	// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
	45	pTrack->GetXYZ(xb);
	46	pTrack->GetPxPyPz(pb);
	47	Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters
	48	AliDebug(1,Form("Track %i pmod=%f mass=%f stat=%i",iTrackN,pTrack->GetP(),pTrack->GetMass(),status));
998b831f	49	x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]);
	50	AliRICHHelix helix(x0,p0,pTrack->GetSign(),b);
	51	Int_t iChamber=helix.RichIntersect(pRich->P());
	52	AliDebug(1,Form("intersection with %i chamber found",iChamber));
	53	if(!iChamber) continue;//intersection with no chamber found
	54
	55	Double_t distMip=9999; //min distance between clusters and track position on PC
	56	Int_t iMipId=0; //index of that min distance cluster
	57	for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
	58	AliRICHcluster pClus=(AliRICHcluster)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
	59	Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
	60	if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track
	61
	62	AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
	63	helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
	64	}////clusters loop for intersected chamber
	65
	66	AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
	67
	68	AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
	69	Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
	70	AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
	71	pTrack->SetRICHsignal(thetaCerenkov);
	72
30c60010	73	Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
	74	CalcProb(thetaCerenkov,p0.Mag(),richPID);
	75	pTrack->SetRICHpid(richPID);
998b831f	76
	77	}//ESD tracks loop
	78	AliDebug(1,"Stop.");
30c60010	79	} //RecWithESD
	80	//__________________________________________________________________________________________________
	81	void AliRICHTracker::RecWithStack()
	82	{
	83	// reconstruction for particles from STACK
	84	//
	85	AliRICH pRich=((AliRICH)gAlice->GetDetector("RICH"));
	86
	87	pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
	88	AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
	89	Int_t iNtracks=pStack->GetNtrack();
	90
	91	Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
	92	AliDebug(1,Form("Start with simulated %i tracks in %f Tesla field",iNtracks,b));
	93	TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
	94
	95
	96	for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
	97	TParticle *pParticle = pStack->Particle(iTrackN);
	98	p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());
	99	x0.SetXYZ(pParticle->Vx(),pParticle->Vy(),pParticle->Vz());
	100	AliRICHHelix helix(x0,p0,Int_t(pParticle->GetPDG()->Charge()/3+0.1),b);
	101	Int_t iChamber=helix.RichIntersect(pRich->P());
	102	AliDebug(1,Form("intersection with %i chamber found",iChamber));
	103	if(!iChamber) continue;//intersection with no chamber found
	104
	105	Double_t distMip=9999; //min distance between clusters and track position on PC
	106	Int_t iMipId=0; //index of that min distance cluster
	107	for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
	108	AliRICHcluster pClus=(AliRICHcluster)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
	109	Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
	110	if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track
	111
	112	AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
	113	helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
	114	}////clusters loop for intersected chamber
	115
	116	AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
	117
	118	AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
	119	Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
	120	AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
	121	// pTrack->SetRICHsignal(thetaCerenkov);
	122
	123	Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p)
	124	CalcProb(thetaCerenkov,p0.Mag(),richPID);
	125
	126	}//ESD tracks loop
	127	AliDebug(1,"Stop.");
	128	} //RecWithStack
998b831f	129
	130	Int_t AliRICHTracker::LoadClusters(TTree *pTree)
	131	{
	132	// Load clusters for RICH
	133	AliDebug(1,"Start."); pTree->GetEntry(0); AliDebug(1,"Stop."); return 0;
	134	}
30c60010	135
	136	//__________________________________________________________________________________________________
	137	void AliRICHTracker::CalcProb(Double_t thetaCer,Double_t pmod, Double_t *richPID)
	138	{
	139	//
	140	Double_t height[5];Double_t totalHeight=0;
	141	Int_t code[5]={kElectron,kMuonPlus,kPiPlus,kKPlus,kProton};
	142	TDatabasePDG *db = new TDatabasePDG();
	143	for(Int_t iPart=0;iPart<4;iPart++){
	144	Double_t mass = db->GetParticle(code[iPart])->Mass();
	145	Double_t refIndex=AliRICHParam::IndOfRefC6F14(6.755);
	146	Double_t cosThetaTh = TMath::Sqrt(massmass+pmodpmod)/(refIndex*pmod);
	147	if(cosThetaTh>=1) {break;}
	148	Double_t thetaTh = TMath::ACos(cosThetaTh);
	149	Double_t sinThetaThNorm = TMath::Sin(thetaTh)/TMath::Sqrt(1-1/(refIndex*refIndex));
	150	Double_t sigmaThetaTh = (0.014(1/sinThetaThNorm-1) + 0.0043)1.25;
	151	height[iPart] = TMath::Gaus(thetaCer,thetaTh,sigmaThetaTh);
	152	totalHeight +=height[iPart];
	153	}
	154	for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;
	155	}//CalcProb