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

/////////////////////////////////////////////////////////////////////////////////
//                                                                             //
// AliFemtoPairCutRadialDistance - 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                                                  //
//                                                                             //
/////////////////////////////////////////////////////////////////////////////////
/********************************************************************************
 *
 * Author: Johanna Gramling, University of Heidelberg, jgramlin@cern.ch
 *         Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch
 *         Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch
 *
 ********************************************************************************/

#include "AliFemtoPairCutRadialDistanceKK.h"
#include <string>
#include <cstdio>

#ifdef __ROOT__
ClassImp(AliFemtoPairCutRadialDistanceKK)
#endif

//__________________
AliFemtoPairCutRadialDistanceKK::AliFemtoPairCutRadialDistanceKK():
  AliFemtoPairCutAntiGamma(),
  fDPhiStarMin(0),
  fEtaMin(0),
  fMinRad(0.8),
  fMagSign(1)
{
}
//__________________
AliFemtoPairCutRadialDistanceKK::AliFemtoPairCutRadialDistanceKK(const AliFemtoPairCutRadialDistanceKK& c) : 
  AliFemtoPairCutAntiGamma(c),
  fDPhiStarMin(0), 
  fEtaMin(0),
  fMinRad(0.8),
  fMagSign(1)
{ 
  fDPhiStarMin = c.fDPhiStarMin;
  fEtaMin = c.fEtaMin;
  fMinRad = c.fMinRad;
  fMagSign = c.fMagSign;
}

//__________________
AliFemtoPairCutRadialDistanceKK::~AliFemtoPairCutRadialDistanceKK(){
  /* no-op */
}
AliFemtoPairCutRadialDistanceKK& AliFemtoPairCutRadialDistanceKK::operator=(const AliFemtoPairCutRadialDistanceKK& c)
{
  if (this != &c) {
    fDPhiStarMin = c.fDPhiStarMin;
    fEtaMin = c.fEtaMin;
    fMinRad = c.fMinRad;
    fMagSign = c.fMagSign;
  }

  return *this;
}
//__________________
bool AliFemtoPairCutRadialDistanceKK::Pass(const AliFemtoPair* pair){
  // Accept pairs based on their TPC entrance separation and
  // quality and sharity
  //  bool temp = true;
  
    double PI = 3.14159265358979312;
//    double pit = 6.28318530717958623;


 Bool_t pass5 = kTRUE;
  
  double phi1 = pair->Track1()->Track()->P().Phi();
  double phi2 = pair->Track2()->Track()->P().Phi();
  double chg1 = pair->Track1()->Track()->Charge();
  double chg2 = pair->Track2()->Track()->Charge();
  double ptv1 = pair->Track1()->Track()->Pt();
  double ptv2 = pair->Track2()->Track()->Pt();
  double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
  double eta2 = pair->Track2()->Track()->P().PseudoRapidity();


    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());


    AliAODEvent *fAOD;

    if (!aodH) {
        // AliWarning("Could not get AODInputHandler");
        return false;
    }
    else {

        fAOD = aodH->GetEvent();
    }
    
     Double_t Bfield = fAOD->GetMagneticField();
  

    if(fabs(eta1-eta2) < fEtaMin) {
  
     // propagate through B field to r=1m
    Double_t phistar1 = phi1 - asin(chg1*(0.1*Bfield)*0.15/ptv1);// 0.15 for D=1m
    if(phistar1 > 2*PI) phistar1 -= 2*PI;
    if(phistar1 < 0) phistar1 += 2*PI;
    Double_t phistar2 = phi2 - asin(chg2*(0.1*Bfield)*0.15/ptv2);// 0.15 for D=1m 
    if(phistar2 > 2*PI) phistar2 -= 2*PI;
    if(phistar2 < 0) phistar2 += 2*PI;

    Double_t deltaphi = phistar1 - phistar2;
    if(deltaphi > PI) deltaphi -= 2*PI;
    if(deltaphi < -PI) deltaphi += 2*PI;
    deltaphi = fabs(deltaphi);

   
     if(deltaphi < fDPhiStarMin){ 
    //  cout<<"---1--- phi1 = "<<phi1<<" phi2 = "<<phi2<<endl;
    //  cout<<" eta1 = "<<eta1<<" eta2 = "<<eta2<<endl;
    //  cout<<" ptv1 = "<<ptv1<<" ptv2 = "<<ptv2<<endl;
    //   cout<<"----1-- "<< deltaphi<<" deltaeta"<< fabs(eta1-eta2)<<endl;
       pass5 = kFALSE;
//	break;
     }
     else
     {
    // propagate through B field to r=1.6m
    phistar1 = phi1 - asin(chg1*(0.1*Bfield)*0.24/ptv1);// mine. 0.24 for D=1.6m
    if(phistar1 > 2*PI) phistar1 -= 2*PI;
    if(phistar1 < 0) phistar1 += 2*PI;
    phistar2 = phi2 - asin(chg2*(0.1*Bfield)*0.24/ptv2);// mine. 0.24 for D=1.6m 
    if(phistar2 > 2*PI) phistar2 -= 2*PI;
    if(phistar2 < 0) phistar2 += 2*PI;

    deltaphi = phistar1 - phistar2;
    if(deltaphi > PI) deltaphi -= 2*PI;
    if(deltaphi < -PI) deltaphi += 2*PI;
    deltaphi = fabs(deltaphi);

    if(deltaphi < fDPhiStarMin){ 	
//      cout<<"---2--- phi1 = "<<phi1<<" phi2 = "<<phi2<<endl;
//      cout<<" eta1 = "<<eta1<<" eta2 = "<<eta2<<endl;
//      cout<<" ptv1 = "<<ptv1<<" ptv2 = "<<ptv2<<endl;
//      cout<<"----2-- "<< deltaphi<<" deltaeta"<< fabs(eta1-eta2)<<endl;
        pass5 = kFALSE;
//	break;
      }
      else
      {
       pass5 = kTRUE;
      }
    
    }  

 }

  if (pass5) {
    pass5 = AliFemtoPairCutAntiGamma::Pass(pair);
  }
  else
    fNPairsFailed++;    

  return pass5;
}
//__________________
AliFemtoString AliFemtoPairCutRadialDistanceKK::Report(){
  // Prepare a report from the execution
  string stemp = "AliFemtoRadialDistance Pair Cut - remove shared and split pairs and pairs with small separation at the specified radius\n";  char ctemp[100];
  snprintf(ctemp , 100, "Accept pair with separation more that %f",fDPhiStarMin);
  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 *AliFemtoPairCutRadialDistanceKK::ListSettings()
{
  // return a list of settings in a writable form
  TList *tListSetttings =  AliFemtoPairCutAntiGamma::ListSettings();
  char buf[200];
  snprintf(buf, 200, "AliFemtoPairCutRadialDistanceKK.phistarsepmin=%f", fDPhiStarMin);
  tListSetttings->AddLast(new TObjString(buf));

  return tListSetttings;
}

void AliFemtoPairCutRadialDistanceKK::SetPhiStarDifferenceMinimum(double dtpc)
{
  fDPhiStarMin = dtpc;
}

void AliFemtoPairCutRadialDistanceKK::SetEtaDifferenceMinimum(double etpc) 
{
  fEtaMin = etpc;
}


void AliFemtoPairCutRadialDistanceKK::SetMinimumRadius(double minrad) 
{
  fMinRad = minrad;
}

void AliFemtoPairCutRadialDistanceKK::SetMagneticFieldSign(int magsign)
{
  if(magsign>1) fMagSign = 1;
  else if(magsign<1) fMagSign = -1;
  else fMagSign = magsign;
}
Commit	Line	Data
78153016	1	/////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoPairCutRadialDistance - 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	*
	11	* Author: Johanna Gramling, University of Heidelberg, jgramlin@cern.ch
	12	* Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch
	13	* Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch
	14	*
	15	********************************************************************************/
	16
	17	#include "AliFemtoPairCutRadialDistanceKK.h"
	18	#include <string>
	19	#include <cstdio>
	20
	21	#ifdef __ROOT__
	22	ClassImp(AliFemtoPairCutRadialDistanceKK)
	23	#endif
	24
	25	//__________________
	26	AliFemtoPairCutRadialDistanceKK::AliFemtoPairCutRadialDistanceKK():
	27	AliFemtoPairCutAntiGamma(),
	28	fDPhiStarMin(0),
	29	fEtaMin(0),
	30	fMinRad(0.8),
	31	fMagSign(1)
	32	{
	33	}
	34	//__________________
	35	AliFemtoPairCutRadialDistanceKK::AliFemtoPairCutRadialDistanceKK(const AliFemtoPairCutRadialDistanceKK& c) :
	36	AliFemtoPairCutAntiGamma(c),
	37	fDPhiStarMin(0),
	38	fEtaMin(0),
	39	fMinRad(0.8),
	40	fMagSign(1)
	41	{
	42	fDPhiStarMin = c.fDPhiStarMin;
	43	fEtaMin = c.fEtaMin;
	44	fMinRad = c.fMinRad;
	45	fMagSign = c.fMagSign;
	46	}
	47
	48	//__________________
	49	AliFemtoPairCutRadialDistanceKK::~AliFemtoPairCutRadialDistanceKK(){
	50	/* no-op */
	51	}
	52	AliFemtoPairCutRadialDistanceKK& AliFemtoPairCutRadialDistanceKK::operator=(const AliFemtoPairCutRadialDistanceKK& c)
	53	{
	54	if (this != &c) {
	55	fDPhiStarMin = c.fDPhiStarMin;
	56	fEtaMin = c.fEtaMin;
	57	fMinRad = c.fMinRad;
	58	fMagSign = c.fMagSign;
	59	}
	60
	61	return *this;
	62	}
	63	//__________________
	64	bool AliFemtoPairCutRadialDistanceKK::Pass(const AliFemtoPair* pair){
65	// Accept pairs based on their TPC entrance separation and
66	// quality and sharity
67	// bool temp = true;
68
69	double PI = 3.14159265358979312;
70	// double pit = 6.28318530717958623;
71
72
73	Bool_t pass5 = kTRUE;
74
75	double phi1 = pair->Track1()->Track()->P().Phi();
76	double phi2 = pair->Track2()->Track()->P().Phi();
77	double chg1 = pair->Track1()->Track()->Charge();
78	double chg2 = pair->Track2()->Track()->Charge();
79	double ptv1 = pair->Track1()->Track()->Pt();
80	double ptv2 = pair->Track2()->Track()->Pt();
81	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
82	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
83
84
85	AliAODInputHandler aodH = dynamic_cast<AliAODInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
86
87
88	AliAODEvent *fAOD;
89
90	if (!aodH) {
91	// AliWarning("Could not get AODInputHandler");
92	return false;
93	}
94	else {
95
96	fAOD = aodH->GetEvent();
97	}
98
99	Double_t Bfield = fAOD->GetMagneticField();
100
101
102	if(fabs(eta1-eta2) < fEtaMin) {
103
104	// propagate through B field to r=1m
105	Double_t phistar1 = phi1 - asin(chg1(0.1Bfield)*0.15/ptv1);// 0.15 for D=1m
106	if(phistar1 > 2PI) phistar1 -= 2PI;
107	if(phistar1 < 0) phistar1 += 2*PI;
108	Double_t phistar2 = phi2 - asin(chg2(0.1Bfield)*0.15/ptv2);// 0.15 for D=1m
109	if(phistar2 > 2PI) phistar2 -= 2PI;
110	if(phistar2 < 0) phistar2 += 2*PI;
111
112	Double_t deltaphi = phistar1 - phistar2;
113	if(deltaphi > PI) deltaphi -= 2*PI;
114	if(deltaphi < -PI) deltaphi += 2*PI;
115	deltaphi = fabs(deltaphi);
116
117
118	if(deltaphi < fDPhiStarMin){
119	// cout<<"---1--- phi1 = "<<phi1<<" phi2 = "<<phi2<<endl;
120	// cout<<" eta1 = "<<eta1<<" eta2 = "<<eta2<<endl;
121	// cout<<" ptv1 = "<<ptv1<<" ptv2 = "<<ptv2<<endl;
122	// cout<<"----1-- "<< deltaphi<<" deltaeta"<< fabs(eta1-eta2)<<endl;
123	pass5 = kFALSE;
124	// break;
125	}
126	else
127	{
128	// propagate through B field to r=1.6m
129	phistar1 = phi1 - asin(chg1(0.1Bfield)*0.24/ptv1);// mine. 0.24 for D=1.6m
130	if(phistar1 > 2PI) phistar1 -= 2PI;
131	if(phistar1 < 0) phistar1 += 2*PI;
132	phistar2 = phi2 - asin(chg2(0.1Bfield)*0.24/ptv2);// mine. 0.24 for D=1.6m
133	if(phistar2 > 2PI) phistar2 -= 2PI;
134	if(phistar2 < 0) phistar2 += 2*PI;
135
136	deltaphi = phistar1 - phistar2;
137	if(deltaphi > PI) deltaphi -= 2*PI;
138	if(deltaphi < -PI) deltaphi += 2*PI;
139	deltaphi = fabs(deltaphi);
140
141	if(deltaphi < fDPhiStarMin){
142	// cout<<"---2--- phi1 = "<<phi1<<" phi2 = "<<phi2<<endl;
143	// cout<<" eta1 = "<<eta1<<" eta2 = "<<eta2<<endl;
144	// cout<<" ptv1 = "<<ptv1<<" ptv2 = "<<ptv2<<endl;
145	// cout<<"----2-- "<< deltaphi<<" deltaeta"<< fabs(eta1-eta2)<<endl;
146	pass5 = kFALSE;
147	// break;
148	}
149	else
150	{
151	pass5 = kTRUE;
152	}
153
154	}
155
156	}
157
158	if (pass5) {
159	pass5 = AliFemtoPairCutAntiGamma::Pass(pair);
160	}
161	else
162	fNPairsFailed++;
163
164	return pass5;
165	}
166	//__________________
167	AliFemtoString AliFemtoPairCutRadialDistanceKK::Report(){
168	// Prepare a report from the execution
169	string stemp = "AliFemtoRadialDistance Pair Cut - remove shared and split pairs and pairs with small separation at the specified radius\n"; char ctemp[100];
170	snprintf(ctemp , 100, "Accept pair with separation more that %f",fDPhiStarMin);
171	stemp += ctemp;
172	snprintf(ctemp , 100, "Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
173	stemp += ctemp;
174	AliFemtoString returnThis = stemp;
175	return returnThis;}
176	//__________________
177
178	TList *AliFemtoPairCutRadialDistanceKK::ListSettings()
179	{
180	// return a list of settings in a writable form
181	TList *tListSetttings = AliFemtoPairCutAntiGamma::ListSettings();
182	char buf[200];
183	snprintf(buf, 200, "AliFemtoPairCutRadialDistanceKK.phistarsepmin=%f", fDPhiStarMin);
184	tListSetttings->AddLast(new TObjString(buf));
185
186	return tListSetttings;
187	}
188
189	void AliFemtoPairCutRadialDistanceKK::SetPhiStarDifferenceMinimum(double dtpc)
190	{
191	fDPhiStarMin = dtpc;
192	}
193
194	void AliFemtoPairCutRadialDistanceKK::SetEtaDifferenceMinimum(double etpc)
195	{
196	fEtaMin = etpc;
197	}
198
199
200	void AliFemtoPairCutRadialDistanceKK::SetMinimumRadius(double minrad)
201	{
202	fMinRad = minrad;
203	}
204
205	void AliFemtoPairCutRadialDistanceKK::SetMagneticFieldSign(int magsign)
206	{
207	if(magsign>1) fMagSign = 1;
208	else if(magsign<1) fMagSign = -1;
209	else fMagSign = magsign;
210	}