Bug fix in HBT code: AliFemtoPairCutRadialDistance class
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoPairCutRadialDistance.cxx
CommitLineData
76ce4b5b 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 "AliFemtoPairCutRadialDistance.h"
18#include <string>
19#include <cstdio>
20
21#ifdef __ROOT__
22ClassImp(AliFemtoPairCutRadialDistance)
23#endif
24
25//__________________
26AliFemtoPairCutRadialDistance::AliFemtoPairCutRadialDistance():
27 AliFemtoPairCutAntiGamma(),
28 fDPhiStarMin(0),
29 fEtaMin(0),
30 fMinRad(0.8),
31 fMagSign(1)
32{
33}
34//__________________
35AliFemtoPairCutRadialDistance::AliFemtoPairCutRadialDistance(const AliFemtoPairCutRadialDistance& 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//__________________
49AliFemtoPairCutRadialDistance::~AliFemtoPairCutRadialDistance(){
50 /* no-op */
51}
52AliFemtoPairCutRadialDistance& AliFemtoPairCutRadialDistance::operator=(const AliFemtoPairCutRadialDistance& 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//__________________
64bool AliFemtoPairCutRadialDistance::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 pih = 3.14159265358979312;
70// double pit = 6.28318530717958623;
71
72
73 double phi1 = pair->Track1()->Track()->P().Phi();
74 double phi2 = pair->Track2()->Track()->P().Phi();
75 double chg1 = pair->Track1()->Track()->Charge();
76 double chg2 = pair->Track2()->Track()->Charge();
77 double ptv1 = pair->Track1()->Track()->Pt();
78 double ptv2 = pair->Track2()->Track()->Pt();
79 double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
80 double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
81
82
5e2038a9 83 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
84
85
86 AliAODEvent *fAOD;
87
88 if (!aodH) {
89 // AliWarning("Could not get AODInputHandler");
90 return false;
91 }
92 else {
93
94 fAOD = aodH->GetEvent();
95 }
96
97 Int_t magsign = fAOD->GetMagneticField();
98 if (magsign > 1)
99 fMagSign = 1;
100 else if ( magsign < 1)
101 fMagSign = -1;
102 else
103 fMagSign = magsign;
104
105
106 // cout << "mag sign = " << fMagSign << endl;
107
76ce4b5b 108 Double_t rad;
109 Bool_t pass5 = kTRUE;
110
111 rad = fMinRad;
1c95f38e 112 for (Double_t iter=fMinRad*100; iter<251; iter+=1.0) {
76ce4b5b 113 Double_t dps = (phi1-phi2+(TMath::ASin(-0.075*chg1*fMagSign*rad/ptv1))-(TMath::ASin(-0.075*chg2*fMagSign*rad/ptv2)));
5e2038a9 114 dps = TVector2::Phi_mpi_pi(dps);
76ce4b5b 115 double etad = eta2 - eta1;
116 if (fabs(etad)<fEtaMin && fabs(dps)<fDPhiStarMin) {
117 // cout << "5% cut is not passed - returning" << endl;
118 pass5 = kFALSE;
119 break;
120 }
121 rad+=0.01;
122 }
123
124
125 if (pass5) {
126 pass5 = AliFemtoPairCutAntiGamma::Pass(pair);
127 }
128 else
129 fNPairsFailed++;
130
131 return pass5;
132}
133//__________________
134AliFemtoString AliFemtoPairCutRadialDistance::Report(){
135 // Prepare a report from the execution
136 string stemp = "AliFemtoRadialDistance Pair Cut - remove shared and split pairs and pairs with small separation at the specified radius\n"; char ctemp[100];
137 snprintf(ctemp , 100, "Accept pair with separation more that %f",fDPhiStarMin);
138 stemp += ctemp;
139 snprintf(ctemp , 100, "Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
140 stemp += ctemp;
141 AliFemtoString returnThis = stemp;
142 return returnThis;}
143//__________________
144
145TList *AliFemtoPairCutRadialDistance::ListSettings()
146{
147 // return a list of settings in a writable form
148 TList *tListSetttings = AliFemtoPairCutAntiGamma::ListSettings();
149 char buf[200];
150 snprintf(buf, 200, "AliFemtoPairCutRadialDistance.phistarsepmin=%f", fDPhiStarMin);
151 tListSetttings->AddLast(new TObjString(buf));
152
153 return tListSetttings;
154}
155
156void AliFemtoPairCutRadialDistance::SetPhiStarDifferenceMinimum(double dtpc)
157{
158 fDPhiStarMin = dtpc;
159}
160
161void AliFemtoPairCutRadialDistance::SetEtaDifferenceMinimum(double etpc)
162{
163 fEtaMin = etpc;
164}
165
166
167void AliFemtoPairCutRadialDistance::SetMinimumRadius(double minrad)
168{
169 fMinRad = minrad;
170}
171
172void AliFemtoPairCutRadialDistance::SetMagneticFieldSign(int magsign)
173{
174 if(magsign>1) fMagSign = 1;
175 else if(magsign<1) fMagSign = -1;
176 else fMagSign = magsign;
177}