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

/////////////////////////////////////////////////////////////////////////////////////
//                                                                                 //
// AliFemtoEventReaderKinematicsChain - the reader class for the Alice ESD and     //
// the model Kinematics information tailored for the Task framework and the        //
// Reads in AliESDfriend to create shared hit/quality information                  //
// Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
//          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
//  	       				                                           //
/////////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoEventReaderKinematicsChain.h"

#include "TFile.h"
#include "TTree.h"
#include "TList.h"

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

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"

#include "TParticle.h"
#include "AliStack.h"
#include "TParticlePDG.h"
#include "AliFemtoModelHiddenInfo.h"
#include "AliFemtoModelGlobalHiddenInfo.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenCocktailEventHeader.h"

#include "AliVertexerTracks.h"

ClassImp(AliFemtoEventReaderKinematicsChain)

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

using namespace std;
//____________________________
AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain():
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fStack(0x0),
  fGenHeader(0x0),
  fRotateToEventPlane(0)
{
  //constructor with 0 parameters , look at default settings 
}

//__________________
AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
  AliFemtoEventReader(aReader),
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fStack(0x0),
  fGenHeader(0x0),
  fRotateToEventPlane(0)
{
  // Copy constructor
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fStack = aReader.fStack;
  fRotateToEventPlane = aReader.fRotateToEventPlane;
}
//__________________
AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
{
  //Destructor
  //delete fEvent;
}

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

  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fStack = aReader.fStack;
  fGenHeader = aReader.fGenHeader;
  fRotateToEventPlane = aReader.fRotateToEventPlane;
  return *this;
}
//__________________
// Simple report
AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
{
  AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
  return temp;
}

//__________________
void AliFemtoEventReaderKinematicsChain::SetConstrained(const bool constrained)
{
  // Select whether to read constrained or not constrained momentum
  fConstrained=constrained;
}
//__________________
bool AliFemtoEventReaderKinematicsChain::GetConstrained() const
{
  // Check whether we read constrained or not constrained momentum
  return fConstrained;
}
//__________________
AliFemtoEvent* AliFemtoEventReaderKinematicsChain::ReturnHbtEvent()
{
  // Get the event, read all the relevant information from the stack
  // and fill the AliFemtoEvent class
  // Returns a valid AliFemtoEvent
  AliFemtoEvent *hbtEvent = 0;
  string tFriendFileName;

  cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
	
  hbtEvent = new AliFemtoEvent;
  //setting basic things
  //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
  hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
  hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
  hbtEvent->SetZDCN1Energy(0);
  hbtEvent->SetZDCP1Energy(0);
  hbtEvent->SetZDCN2Energy(0);
  hbtEvent->SetZDCP2Energy(0);
  hbtEvent->SetZDCEMEnergy(0);
  hbtEvent->SetZDCParticipants(0);
  hbtEvent->SetTriggerMask(0);
  hbtEvent->SetTriggerCluster(0);
	
  //Vertex
  double fV1[3] = {0.0,0.0,0.0};
  double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};


  AliFmThreeVectorF vertex(0,0,0);
    
  
  hbtEvent->SetPrimVertPos(vertex);
  hbtEvent->SetPrimVertCov(fVCov);

  Double_t tReactionPlane = 0;

  AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);	
  if (!hdh) {
    AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
    if (cdh) {
      TList *tGenHeaders = cdh->GetHeaders();
      for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
	hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);	
	if (hdh) break;
      }
    }
  }

  if (hdh)
    {
      tReactionPlane = hdh->ReactionPlaneAngle();
      cout << "Got reaction plane " << tReactionPlane << endl;
    }

  hbtEvent->SetReactionPlaneAngle(tReactionPlane);

  //starting to reading tracks
  int nofTracks=0;  //number of all tracks in MC event
  nofTracks=fStack->GetNtrack(); 
  int realnofTracks=0;//number of track which we use in analysis


  int tNormMult = 0;
  for (int i=0;i<nofTracks;i++)
    {

      //take only primaries
      if(!fStack->IsPhysicalPrimary(i)) {continue;}
	  	  
      AliFemtoTrack* trackCopy = new AliFemtoTrack();	
	
      	  //getting next track
      TParticle *kinetrack= fStack->Particle(i);

      //setting multiplicity
        realnofTracks++;//real number of tracks (only primary particles)

      //setting normalized multiplicity
      if (kinetrack->Eta() < 0.9)
	if(kinetrack->GetPDG()->Charge()/3!=0)
	  tNormMult++;
	  
	  
	  //charge
      trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));


      //in aliroot we have AliPID 
      //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
      //we use only 5 first
      double kinepid[5];
      for(int pid_iter=0;pid_iter<5;pid_iter++)
	  kinepid[pid_iter]=0;

      int pdgcode = kinetrack->GetPdgCode();
      //proton
      if(pdgcode==2212 || pdgcode==-2212)
        kinepid[4]=1000;
      //kaon
      if(pdgcode==321 || pdgcode==-321 )
        kinepid[3]=1000;
      //pion
      if( pdgcode==211 || pdgcode==-211)
        kinepid[2]=1000;
      //electron
      if(pdgcode==11 || pdgcode==-11)
        kinepid[0]=1000;
      //muon
      if(pdgcode==13 || pdgcode==-13)
        kinepid[1]=1000;

      trackCopy->SetPidProbElectron(kinepid[0]);
      trackCopy->SetPidProbMuon(kinepid[1]);
      trackCopy->SetPidProbPion(kinepid[2]);
      trackCopy->SetPidProbKaon(kinepid[3]);
      trackCopy->SetPidProbProton(kinepid[4]);
					
					
	//Momentum
      double pxyz[3];
      double rxyz[3];
     
	pxyz[0]=kinetrack->Px();
	pxyz[1]=kinetrack->Py();
	pxyz[2]=kinetrack->Pz();

	rxyz[0]=kinetrack->Vx();
	rxyz[1]=kinetrack->Vy();
	rxyz[2]=kinetrack->Vz();

	if (fRotateToEventPlane) {
	  double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
	  double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
	  
	  pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
	  pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
	}

	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	if (v.Mag() < 0.0001) {
	  //	  cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
	  delete trackCopy;
	  continue;
	}

	trackCopy->SetP(v);//setting momentum
	trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);

	//label
    trackCopy->SetLabel(i);


	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	
		
    }
  
  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
  hbtEvent->SetNormalizedMult(tNormMult);
  fCurEvent++;	

  return hbtEvent; 
}

//___________________
void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
{
  // The chain loads the stack for us
  // You must provide the address where it can be found
  fStack = aStack;
}
//___________________
void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
{
  // The chain loads the generator event header for us
  // You must provide the address where it can be found
  fGenHeader = aGenHeader;
}

//__________________
void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
{
  fRotateToEventPlane=dorotate;
}

Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double *impact, double *covar)
{
  // Calculates the number of sigma to the vertex.

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

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

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

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

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

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

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

  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  return d;
}
Commit	Line	Data
859f1d8d	1	/////////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoEventReaderKinematicsChain - the reader class for the Alice ESD and //
	4	// the model Kinematics information tailored for the Task framework and the //
	5	// Reads in AliESDfriend to create shared hit/quality information //
	6	// Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch //
	7	// Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
	8	// //
	9	/////////////////////////////////////////////////////////////////////////////////////
	10
	11	#include "AliFemtoEventReaderKinematicsChain.h"
	12
	13	#include "TFile.h"
	14	#include "TTree.h"
	15	#include "TList.h"
	16
	17	#include "AliFmPhysicalHelixD.h"
	18	#include "AliFmThreeVectorF.h"
	19
	20	#include "SystemOfUnits.h"
	21
	22	#include "AliFemtoEvent.h"
	23
	24	#include "TParticle.h"
	25	#include "AliStack.h"
	26	#include "TParticlePDG.h"
	27	#include "AliFemtoModelHiddenInfo.h"
	28	#include "AliFemtoModelGlobalHiddenInfo.h"
	29	#include "AliGenHijingEventHeader.h"
	30	#include "AliGenCocktailEventHeader.h"
	31
	32	#include "AliVertexerTracks.h"
	33
	34	ClassImp(AliFemtoEventReaderKinematicsChain)
	35
	36	#if !(ST_NO_NAMESPACES)
	37	using namespace units;
	38	#endif
	39
	40	using namespace std;
	41	//____________________________
	42	AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain():
	43	fFileName(" "),
	44	fConstrained(true),
	45	fNumberofEvent(0),
	46	fCurEvent(0),
	47	fCurFile(0),
	48	fStack(0x0),
	49	fGenHeader(0x0),
	50	fRotateToEventPlane(0)
	51	{
	52	//constructor with 0 parameters , look at default settings
	53	}
	54
	55	//__________________
	56	AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
	57	AliFemtoEventReader(aReader),
	58	fFileName(" "),
	59	fConstrained(true),
	60	fNumberofEvent(0),
	61	fCurEvent(0),
	62	fCurFile(0),
	63	fStack(0x0),
	64	fGenHeader(0x0),
65	fRotateToEventPlane(0)
66	{
67	// Copy constructor
68	fConstrained = aReader.fConstrained;
69	fNumberofEvent = aReader.fNumberofEvent;
70	fCurEvent = aReader.fCurEvent;
71	fCurFile = aReader.fCurFile;
72	fStack = aReader.fStack;
73	fRotateToEventPlane = aReader.fRotateToEventPlane;
74	}
75	//__________________
76	AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
77	{
78	//Destructor
79	//delete fEvent;
80	}
81
82	//__________________
83	AliFemtoEventReaderKinematicsChain& AliFemtoEventReaderKinematicsChain::operator=(const AliFemtoEventReaderKinematicsChain& aReader)
84	{
85	// Assignment operator
86	if (this == &aReader)
87	return *this;
88
89	fConstrained = aReader.fConstrained;
90	fNumberofEvent = aReader.fNumberofEvent;
91	fCurEvent = aReader.fCurEvent;
92	fCurFile = aReader.fCurFile;
93	fStack = aReader.fStack;
94	fGenHeader = aReader.fGenHeader;
95	fRotateToEventPlane = aReader.fRotateToEventPlane;
96	return *this;
97	}
98	//__________________
99	// Simple report
100	AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
101	{
102	AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
103	return temp;
104	}
105
106	//__________________
107	void AliFemtoEventReaderKinematicsChain::SetConstrained(const bool constrained)
108	{
109	// Select whether to read constrained or not constrained momentum
110	fConstrained=constrained;
111	}
112	//__________________
113	bool AliFemtoEventReaderKinematicsChain::GetConstrained() const
114	{
115	// Check whether we read constrained or not constrained momentum
116	return fConstrained;
117	}
118	//__________________
119	AliFemtoEvent* AliFemtoEventReaderKinematicsChain::ReturnHbtEvent()
120	{
121	// Get the event, read all the relevant information from the stack
122	// and fill the AliFemtoEvent class
123	// Returns a valid AliFemtoEvent
124	AliFemtoEvent *hbtEvent = 0;
125	string tFriendFileName;
126
127	cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
128
129	hbtEvent = new AliFemtoEvent;
130	//setting basic things
131	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
132	hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
133	hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
134	hbtEvent->SetZDCN1Energy(0);
135	hbtEvent->SetZDCP1Energy(0);
136	hbtEvent->SetZDCN2Energy(0);
137	hbtEvent->SetZDCP2Energy(0);
138	hbtEvent->SetZDCEMEnergy(0);
139	hbtEvent->SetZDCParticipants(0);
140	hbtEvent->SetTriggerMask(0);
141	hbtEvent->SetTriggerCluster(0);
142
143	//Vertex
b84aaa8b	144	double fV1[3] = {0.0,0.0,0.0};
b84aaa8b	145	double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
859f1d8d	146
	147
	148	AliFmThreeVectorF vertex(0,0,0);
	149
	150
	151	hbtEvent->SetPrimVertPos(vertex);
	152	hbtEvent->SetPrimVertCov(fVCov);
	153
	154	Double_t tReactionPlane = 0;
	155
	156	AliGenHijingEventHeader hdh = dynamic_cast<AliGenHijingEventHeader > (fGenHeader);
	157	if (!hdh) {
	158	AliGenCocktailEventHeader cdh = dynamic_cast<AliGenCocktailEventHeader > (fGenHeader);
	159	if (cdh) {
	160	TList *tGenHeaders = cdh->GetHeaders();
	161	for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
	162	hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
	163	if (hdh) break;
	164	}
	165	}
	166	}
	167
	168	if (hdh)
	169	{
	170	tReactionPlane = hdh->ReactionPlaneAngle();
	171	cout << "Got reaction plane " << tReactionPlane << endl;
	172	}
	173
	174	hbtEvent->SetReactionPlaneAngle(tReactionPlane);
	175
	176	//starting to reading tracks
	177	int nofTracks=0; //number of all tracks in MC event
	178	nofTracks=fStack->GetNtrack();
	179	int realnofTracks=0;//number of track which we use in analysis
	180
	181
	182	int tNormMult = 0;
	183	for (int i=0;i<nofTracks;i++)
	184	{
	185
	186	//take only primaries
	187	if(!fStack->IsPhysicalPrimary(i)) {continue;}
	188
	189	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	190
	191	//getting next track
	192	TParticle *kinetrack= fStack->Particle(i);
	193
	194	//setting multiplicity
	195	realnofTracks++;//real number of tracks (only primary particles)
	196
	197	//setting normalized multiplicity
	198	if (kinetrack->Eta() < 0.9)
	199	if(kinetrack->GetPDG()->Charge()/3!=0)
	200	tNormMult++;
	201
	202
	203	//charge
	204	trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
	205
	206
	207	//in aliroot we have AliPID
	208	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	209	//we use only 5 first
210	double kinepid[5];
211	for(int pid_iter=0;pid_iter<5;pid_iter++)
212	kinepid[pid_iter]=0;
213
214	int pdgcode = kinetrack->GetPdgCode();
215	//proton
216	if(pdgcode==2212 \|\| pdgcode==-2212)
217	kinepid[4]=1000;
218	//kaon
219	if(pdgcode==321 \|\| pdgcode==-321 )
220	kinepid[3]=1000;
221	//pion
222	if( pdgcode==211 \|\| pdgcode==-211)
223	kinepid[2]=1000;
224	//electron
225	if(pdgcode==11 \|\| pdgcode==-11)
226	kinepid[0]=1000;
227	//muon
228	if(pdgcode==13 \|\| pdgcode==-13)
229	kinepid[1]=1000;
230
231	trackCopy->SetPidProbElectron(kinepid[0]);
232	trackCopy->SetPidProbMuon(kinepid[1]);
233	trackCopy->SetPidProbPion(kinepid[2]);
234	trackCopy->SetPidProbKaon(kinepid[3]);
235	trackCopy->SetPidProbProton(kinepid[4]);
236
237
238	//Momentum
239	double pxyz[3];
240	double rxyz[3];
241
242	pxyz[0]=kinetrack->Px();
243	pxyz[1]=kinetrack->Py();
244	pxyz[2]=kinetrack->Pz();
245
246	rxyz[0]=kinetrack->Vx();
247	rxyz[1]=kinetrack->Vy();
248	rxyz[2]=kinetrack->Vz();
249
250	if (fRotateToEventPlane) {
251	double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
252	double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
253
254	pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
255	pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
256	}
257
258	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
259	if (v.Mag() < 0.0001) {
260	// cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
261	delete trackCopy;
262	continue;
263	}
264
265	trackCopy->SetP(v);//setting momentum
266	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
267	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
268	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
269
270	//label
271	trackCopy->SetLabel(i);
272
273
274	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
275
276
277	}
278
279	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
280	hbtEvent->SetNormalizedMult(tNormMult);
281	fCurEvent++;
282
283	return hbtEvent;
284	}
285
286	//___________________
287	void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
288	{
289	// The chain loads the stack for us
290	// You must provide the address where it can be found
291	fStack = aStack;
292	}
293	//___________________
294	void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
295	{
296	// The chain loads the generator event header for us
297	// You must provide the address where it can be found
298	fGenHeader = aGenHeader;
299	}
300
301	//__________________
302	void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
303	{
304	fRotateToEventPlane=dorotate;
305	}
306
307	Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double impact, double covar)
308	{
309	// Calculates the number of sigma to the vertex.
310
311	Float_t b[2];
312	Float_t bRes[2];
313	Float_t bCov[3];
314
315	b[0] = impact[0];
316	b[1] = impact[1];
317	bCov[0] = covar[0];
318	bCov[1] = covar[1];
319	bCov[2] = covar[2];
320
321	bRes[0] = TMath::Sqrt(bCov[0]);
322	bRes[1] = TMath::Sqrt(bCov[2]);
323
324	// -----------------------------------
325	// How to get to a n-sigma cut?
326	//
327	// The accumulated statistics from 0 to d is
328	//
329	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
330	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
331	//
332	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
333	// Can this be expressed in a different way?
334
335	if (bRes[0] == 0 \|\| bRes[1] ==0)
336	return -1;
337
338	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
339
340	// stupid rounding problem screws up everything:
341	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
342	if (TMath::Exp(-d * d / 2) < 1e-10)
343	return 1000;
344
345	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
346	return d;
347	}