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

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// AliFemtoPairCutResonances - a pair cut which checks     //
// for some pair qualities that attempt to identify slit/doubly            //
// reconstructed tracks and also selects pairs based on their separation   //
// at the entrance to the TPC                                              //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "AliFemtoPairCutResonances.h"
#include <string>
#include <cstdio>
#include <TMath.h>

#ifdef __ROOT__
ClassImp(AliFemtoPairCutResonances)
#endif

//__________________
AliFemtoPairCutResonances::AliFemtoPairCutResonances():
   AliFemtoShareQualityPairCut(),
   fMaxEEMinv(0.0),
   fMaxDTheta(0.0),
   fDataType(kAOD),
   fSwitchPassFail(0)
{
}
//__________________
AliFemtoPairCutResonances::AliFemtoPairCutResonances(const AliFemtoPairCutResonances& c) : 
  AliFemtoShareQualityPairCut(c),
  fMaxEEMinv(0.0),
  fMaxDTheta(0.0),
  fDataType(kAOD),
  fSwitchPassFail(0)
{
  fMaxEEMinv = c.fMaxEEMinv;
  fMaxDTheta = c.fMaxDTheta;
  fDataType = c.fDataType;
  fSwitchPassFail=c.fSwitchPassFail;
}

AliFemtoPairCutResonances& AliFemtoPairCutResonances::operator=(const AliFemtoPairCutResonances& c)
{
  if (this != &c) {
    fMaxEEMinv = c.fMaxEEMinv;
    fMaxDTheta = c.fMaxDTheta;
    fDataType = c.fDataType;
    fSwitchPassFail=c.fSwitchPassFail;
  }

  return *this;

}
//__________________
AliFemtoPairCutResonances::~AliFemtoPairCutResonances(){
}
//__________________
bool AliFemtoPairCutResonances::Pass(const AliFemtoPair* pair){
  // Accept pairs based on their TPC entrance separation and
  // quality and sharity
  bool temp = true;

  if(fDataType==kKine)
    return true;

  double me = 0.000511;
  double mPi = 0.13957018;
  double mp = 0.938272046;

  double mK0min = 0.00049;
  double mK0max = 0.00051;
  //double mK0 = 0.000497614;
  double mRhomin = 0.000765;
  double mRhomax = 0.000785;
  //double mRho = 0.00077526;
  double mLmin = 1.095;
  double mLmax = 1.135;
  //double mL = 1.115683;

  if ((pair->Track1()->Track()->Charge() * pair->Track2()->Track()->Charge()) < 0.0) {
    // double theta1 = pair->Track1()->Track()->P().Theta();
    // double theta2 = pair->Track2()->Track()->P().Theta();
    // double dtheta = TMath::Abs(theta1 - theta2);

    // check on ee pairs (gamma)
    double e1 = TMath::Sqrt(me*me + pair->Track1()->Track()->P().Mag2());
    double e2 = TMath::Sqrt(me*me + pair->Track2()->Track()->P().Mag2());
    double minvGamma = 2*me*me + 2*(e1*e2 -
			       pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
			       pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
			       pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
    if ( minvGamma < fMaxEEMinv )
       temp = false;
    //check on resonances
    double pi1 =  TMath::Sqrt(mPi*mPi + pair->Track1()->Track()->P().Mag2());
    double pi2 =  TMath::Sqrt(mPi*mPi + pair->Track2()->Track()->P().Mag2());
    double p1 =  TMath::Sqrt(mp*mp + pair->Track1()->Track()->P().Mag2());
    double p2 =  TMath::Sqrt(mp*mp + pair->Track2()->Track()->P().Mag2());
    //check on K0 and Rho
    double minv2pi = 2*mPi*mPi + 2*(pi1*pi2 -
			       pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
			       pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
			       pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
    if ( ((minv2pi>mK0min && minv2pi<mK0max) || (minv2pi>mRhomin && minv2pi<mRhomax)) )
       temp = false;
    //check on L0
    double minvpPi = 2*mp*mPi + 2*(p1*pi2 -
			       pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
			       pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
			       pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
    double minvPip = 2*mPi*mp + 2*(pi1*p2 -
			       pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
			       pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
			       pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
    if( ((minvpPi>mLmin) && (minvpPi<mLmax)) || ((minvPip>mLmin) && (minvPip<mLmax)) )
       temp = false;
  }
  if (fSwitchPassFail) // choose only resonances
  {
     if (!temp) {
        temp = AliFemtoShareQualityPairCut::Pass(pair);
        if (temp) {fNPairsPassed++;}
        else fNPairsFailed++;
        return temp;
     }
     else
     {
        fNPairsFailed++;
        return false;
     }
  }
  else // cut resonances
  {
     if (temp) {
        temp = AliFemtoShareQualityPairCut::Pass(pair);
        if (temp) {fNPairsPassed++;}
        else fNPairsFailed++;
        return temp;
     }
     else
     {
        fNPairsFailed++;
        return false;
     }
  }
}
//__________________
AliFemtoString AliFemtoPairCutResonances::Report(){
  // Prepare a report from the execution
  string stemp = "AliFemtoPairCutResonances Pair Cut - remove pairs possibly coming from Gamma conversions\n";  
  char ctemp[100];
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of pairs which passed:\t%ld  Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
  stemp += ctemp;
  AliFemtoString returnThis = stemp;
  return returnThis;}
//__________________

TList *AliFemtoPairCutResonances::ListSettings()
{
  // return a list of settings in a writable form
  TList *tListSetttings =  AliFemtoShareQualityPairCut::ListSettings();
  char buf[200];
  snprintf(buf, 200, "AliFemtoPairCutResonances.maxeeminv=%f", fMaxEEMinv);
  snprintf(buf, 200, "AliFemtoPairCutResonances.maxdtheta=%f", fMaxDTheta);
  tListSetttings->AddLast(new TObjString(buf));

  return tListSetttings;
}

void AliFemtoPairCutResonances::SetMaxEEMinv(Double_t maxeeminv)
{
  fMaxEEMinv = maxeeminv;
}

void AliFemtoPairCutResonances::SetMaxThetaDiff(Double_t maxdtheta)
{
  fMaxDTheta = maxdtheta;
}

void AliFemtoPairCutResonances::SetDataType(AliFemtoDataType type)
{
  fDataType = type;
}

void AliFemtoPairCutResonances::SetChooseResonances(bool onlyResonances)
{
  fSwitchPassFail = onlyResonances;
}
Commit	Line	Data
2e04885f	1	/////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoPairCutResonances - a pair cut which checks //
	4	// for some pair qualities that attempt to identify slit/doubly //
	5	// reconstructed tracks and also selects pairs based on their separation //
	6	// at the entrance to the TPC //
	7	// //
	8	/////////////////////////////////////////////////////////////////////////////
	9
	10	#include "AliFemtoPairCutResonances.h"
	11	#include <string>
	12	#include <cstdio>
	13	#include <TMath.h>
	14
	15	#ifdef __ROOT__
	16	ClassImp(AliFemtoPairCutResonances)
	17	#endif
	18
	19	//__________________
	20	AliFemtoPairCutResonances::AliFemtoPairCutResonances():
	21	AliFemtoShareQualityPairCut(),
	22	fMaxEEMinv(0.0),
	23	fMaxDTheta(0.0),
	24	fDataType(kAOD),
	25	fSwitchPassFail(0)
	26	{
	27	}
	28	//__________________
	29	AliFemtoPairCutResonances::AliFemtoPairCutResonances(const AliFemtoPairCutResonances& c) :
	30	AliFemtoShareQualityPairCut(c),
	31	fMaxEEMinv(0.0),
	32	fMaxDTheta(0.0),
	33	fDataType(kAOD),
	34	fSwitchPassFail(0)
	35	{
	36	fMaxEEMinv = c.fMaxEEMinv;
	37	fMaxDTheta = c.fMaxDTheta;
	38	fDataType = c.fDataType;
	39	fSwitchPassFail=c.fSwitchPassFail;
	40	}
	41
	42	AliFemtoPairCutResonances& AliFemtoPairCutResonances::operator=(const AliFemtoPairCutResonances& c)
	43	{
	44	if (this != &c) {
	45	fMaxEEMinv = c.fMaxEEMinv;
	46	fMaxDTheta = c.fMaxDTheta;
	47	fDataType = c.fDataType;
	48	fSwitchPassFail=c.fSwitchPassFail;
	49	}
	50
	51	return *this;
	52
	53	}
	54	//__________________
	55	AliFemtoPairCutResonances::~AliFemtoPairCutResonances(){
	56	}
	57	//__________________
	58	bool AliFemtoPairCutResonances::Pass(const AliFemtoPair* pair){
	59	// Accept pairs based on their TPC entrance separation and
	60	// quality and sharity
	61	bool temp = true;
	62
	63	if(fDataType==kKine)
	64	return true;
65
66	double me = 0.000511;
67	double mPi = 0.13957018;
68	double mp = 0.938272046;
69
70	double mK0min = 0.00049;
71	double mK0max = 0.00051;
0cf72019	72	//double mK0 = 0.000497614;
2e04885f	73	double mRhomin = 0.000765;
2e04885f	74	double mRhomax = 0.000785;
0cf72019	75	//double mRho = 0.00077526;
2e04885f	76	double mLmin = 1.095;
2e04885f	77	double mLmax = 1.135;
0cf72019	78	//double mL = 1.115683;
2e04885f	79
2e04885f	80	if ((pair->Track1()->Track()->Charge() * pair->Track2()->Track()->Charge()) < 0.0) {
0cf72019	81	// double theta1 = pair->Track1()->Track()->P().Theta();
	82	// double theta2 = pair->Track2()->Track()->P().Theta();
	83	// double dtheta = TMath::Abs(theta1 - theta2);
2e04885f	84
	85	// check on ee pairs (gamma)
	86	double e1 = TMath::Sqrt(me*me + pair->Track1()->Track()->P().Mag2());
	87	double e2 = TMath::Sqrt(me*me + pair->Track2()->Track()->P().Mag2());
	88	double minvGamma = 2meme + 2(e1e2 -
	89	pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
	90	pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
	91	pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
	92	if ( minvGamma < fMaxEEMinv )
	93	temp = false;
	94	//check on resonances
	95	double pi1 = TMath::Sqrt(mPi*mPi + pair->Track1()->Track()->P().Mag2());
	96	double pi2 = TMath::Sqrt(mPi*mPi + pair->Track2()->Track()->P().Mag2());
	97	double p1 = TMath::Sqrt(mp*mp + pair->Track1()->Track()->P().Mag2());
	98	double p2 = TMath::Sqrt(mp*mp + pair->Track2()->Track()->P().Mag2());
	99	//check on K0 and Rho
	100	double minv2pi = 2mPimPi + 2(pi1pi2 -
	101	pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
	102	pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
	103	pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
	104	if ( ((minv2pi>mK0min && minv2pi<mK0max) \|\| (minv2pi>mRhomin && minv2pi<mRhomax)) )
	105	temp = false;
	106	//check on L0
	107	double minvpPi = 2mpmPi + 2(p1pi2 -
	108	pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
	109	pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
	110	pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
	111	double minvPip = 2mPimp + 2(pi1p2 -
	112	pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
	113	pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
	114	pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
	115	if( ((minvpPi>mLmin) && (minvpPi<mLmax)) \|\| ((minvPip>mLmin) && (minvPip<mLmax)) )
	116	temp = false;
	117	}
	118	if (fSwitchPassFail) // choose only resonances
	119	{
	120	if (!temp) {
	121	temp = AliFemtoShareQualityPairCut::Pass(pair);
	122	if (temp) {fNPairsPassed++;}
	123	else fNPairsFailed++;
	124	return temp;
	125	}
	126	else
	127	{
	128	fNPairsFailed++;
	129	return false;
	130	}
	131	}
	132	else // cut resonances
	133	{
	134	if (temp) {
	135	temp = AliFemtoShareQualityPairCut::Pass(pair);
	136	if (temp) {fNPairsPassed++;}
	137	else fNPairsFailed++;
	138	return temp;
	139	}
	140	else
	141	{
	142	fNPairsFailed++;
	143	return false;
	144	}
	145	}
	146	}
	147	//__________________
148	AliFemtoString AliFemtoPairCutResonances::Report(){
149	// Prepare a report from the execution
150	string stemp = "AliFemtoPairCutResonances Pair Cut - remove pairs possibly coming from Gamma conversions\n";
151	char ctemp[100];
152	stemp += ctemp;
153	snprintf(ctemp , 100, "Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
154	stemp += ctemp;
155	AliFemtoString returnThis = stemp;
156	return returnThis;}
157	//__________________
158
159	TList *AliFemtoPairCutResonances::ListSettings()
160	{
161	// return a list of settings in a writable form
162	TList *tListSetttings = AliFemtoShareQualityPairCut::ListSettings();
163	char buf[200];
164	snprintf(buf, 200, "AliFemtoPairCutResonances.maxeeminv=%f", fMaxEEMinv);
165	snprintf(buf, 200, "AliFemtoPairCutResonances.maxdtheta=%f", fMaxDTheta);
166	tListSetttings->AddLast(new TObjString(buf));
167
168	return tListSetttings;
169	}
170
171	void AliFemtoPairCutResonances::SetMaxEEMinv(Double_t maxeeminv)
172	{
173	fMaxEEMinv = maxeeminv;
174	}
175
176	void AliFemtoPairCutResonances::SetMaxThetaDiff(Double_t maxdtheta)
177	{
178	fMaxDTheta = maxdtheta;
179	}
180
181	void AliFemtoPairCutResonances::SetDataType(AliFemtoDataType type)
182	{
183	fDataType = type;
184	}
185
186	void AliFemtoPairCutResonances::SetChooseResonances(bool onlyResonances)
187	{
188	fSwitchPassFail = onlyResonances;
189	}