6d79abb93b7a1532b69a49f1d76d0563ff52a312
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoPairCutRadialDistance.cxx
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__
22 ClassImp(AliFemtoPairCutRadialDistance)
23 #endif
24
25 //__________________
26 AliFemtoPairCutRadialDistance::AliFemtoPairCutRadialDistance():
27   AliFemtoPairCutAntiGamma(),
28   fDPhiStarMin(0),
29   fEtaMin(0),
30   fMinRad(0.8),
31   fMagSign(1)
32 {
33 }
34 //__________________
35 AliFemtoPairCutRadialDistance::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 //__________________
49 AliFemtoPairCutRadialDistance::~AliFemtoPairCutRadialDistance(){
50   /* no-op */
51 }
52 AliFemtoPairCutRadialDistance& 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 //__________________
64 bool 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
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
108   Double_t rad;
109   Bool_t pass5 = kTRUE;
110
111     rad = fMinRad;
112     for (Double_t iter=fMinRad*10; iter<251; iter+=1.0) {
113       Double_t dps = (phi1-phi2+(TMath::ASin(-0.075*chg1*fMagSign*rad/ptv1))-(TMath::ASin(-0.075*chg2*fMagSign*rad/ptv2)));
114         dps = TVector2::Phi_mpi_pi(dps);
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 //__________________
134 AliFemtoString 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
145 TList *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
156 void AliFemtoPairCutRadialDistance::SetPhiStarDifferenceMinimum(double dtpc)
157 {
158   fDPhiStarMin = dtpc;
159 }
160
161 void AliFemtoPairCutRadialDistance::SetEtaDifferenceMinimum(double etpc) 
162 {
163   fEtaMin = etpc;
164 }
165
166
167 void AliFemtoPairCutRadialDistance::SetMinimumRadius(double minrad) 
168 {
169   fMinRad = minrad;
170 }
171
172 void AliFemtoPairCutRadialDistance::SetMagneticFieldSign(int magsign)
173 {
174   if(magsign>1) fMagSign = 1;
175   else if(magsign<1) fMagSign = -1;
176   else fMagSign = magsign;
177 }