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   fMaxRad(2.5),
  32   fMagSign(1),
  33   fPhistarmin(kTRUE)
  34 {
  35 }
  36 //__________________
  37 AliFemtoPairCutRadialDistance::AliFemtoPairCutRadialDistance(const AliFemtoPairCutRadialDistance& c) :
  38   AliFemtoPairCutAntiGamma(c),
  39   fDPhiStarMin(0),
  40   fEtaMin(0),
  41   fMinRad(0.8),
  42   fMaxRad(2.5),
  43   fMagSign(1),
  44   fPhistarmin(kTRUE)
  45 {
  46   fDPhiStarMin = c.fDPhiStarMin;
  47   fEtaMin = c.fEtaMin;
  48   fMinRad = c.fMinRad;
  49   fMaxRad = c.fMaxRad;
  50   fMagSign = c.fMagSign;
  51   fPhistarmin = c.fPhistarmin;
  52 }
  53
  54 //__________________
  55 AliFemtoPairCutRadialDistance::~AliFemtoPairCutRadialDistance(){
  56   /* no-op */
  57 }
  58 AliFemtoPairCutRadialDistance& AliFemtoPairCutRadialDistance::operator=(const AliFemtoPairCutRadialDistance& c)
  59 {
  60   if (this != &c) {
  61     fDPhiStarMin = c.fDPhiStarMin;
  62     fEtaMin = c.fEtaMin;
  63     fMinRad = c.fMinRad;
  64     fMaxRad = c.fMaxRad;
  65     fMagSign = c.fMagSign;
  66     fPhistarmin = c.fPhistarmin;
  67
  68   }
  69
  70   return *this;
  71 }
  72 //__________________
  73 bool AliFemtoPairCutRadialDistance::Pass(const AliFemtoPair* pair){
  74   // Accept pairs based on their TPC entrance separation and
  75   // quality and sharity
  76   //  bool temp = true;
  77
  78 //    double pih = 3.14159265358979312;
  79 //    double pit = 6.28318530717958623;
  80
  81   double phi1 = pair->Track1()->Track()->P().Phi();
  82   double phi2 = pair->Track2()->Track()->P().Phi();
  83   double chg1 = pair->Track1()->Track()->Charge();
  84   double chg2 = pair->Track2()->Track()->Charge();
  85   double ptv1 = pair->Track1()->Track()->Pt();
  86   double ptv2 = pair->Track2()->Track()->Pt();
  87   double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
  88   double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
  89
  90
  91   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  92   Double_t magsign = 0.0;
  93
  94   if (!aodH) {
  95     // AliWarning("Could not get AODInputHandler");
  96     return false;
  97   }
  98   else {
  99     AliAODEvent *fAOD;
 100     fAOD = aodH->GetEvent();
 101     //cout<<fAOD<<endl;
 102     magsign = fAOD->GetMagneticField();
 103   }
 104
 105   if (magsign > 1)
 106     fMagSign = 1;
 107   else if ( magsign < 1)
 108     fMagSign = -1;
 109   else
 110     fMagSign = magsign;
 111
 112   //cout << "mag sign = " << magsign << endl;
 113
 114   Double_t rad;
 115   Bool_t pass5 = kTRUE;
 116
 117   rad = fMinRad;
 118
 119   if (fPhistarmin) {
 120     for (rad = fMinRad; rad < fMaxRad; rad += 0.01) {
 121       Double_t dps = (phi2-phi1+(TMath::ASin(-0.075*chg2*fMagSign*rad/ptv2))-(TMath::ASin(-0.075*chg1*fMagSign*rad/ptv1)));
 122       dps = TVector2::Phi_mpi_pi(dps);
 123       Double_t etad = eta2 - eta1;
 124       if (fabs(etad)<fEtaMin && fabs(dps)<fDPhiStarMin) {
 125         // cout << "5% cut is not passed - returning" << endl;
 126         pass5 = kFALSE;
 127         break;
 128       }
 129     }
 130   }
 131   else {
 132
 133     double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
 134     double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
 135
 136     if (fabs(afsi0b) >=1.) return kTRUE;
 137     if (fabs(afsi1b) >=1.) return kTRUE;
 138
 139     Double_t dps =  phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
 140     dps = TVector2::Phi_mpi_pi(dps);
 141
 142     Double_t etad = eta2 - eta1;
 143     if (fabs(etad)<fEtaMin && fabs(dps)<fDPhiStarMin) {
 144       pass5 = kFALSE;
 145     }
 146
 147   }
 148
 149   if (pass5) {
 150     pass5 = AliFemtoPairCutAntiGamma::Pass(pair);
 151   }
 152   else {
 153     fNPairsFailed++;
 154   }
 155
 156   return pass5;
 157 }
 158 //__________________
 159 AliFemtoString AliFemtoPairCutRadialDistance::Report(){
 160   // Prepare a report from the execution
 161   string stemp = "AliFemtoRadialDistance Pair Cut - remove shared and split pairs and pairs with small separation at the specified radius\n";  char ctemp[100];
 162   snprintf(ctemp , 100, "Accept pair with separation more that %f",fDPhiStarMin);
 163   stemp += ctemp;
 164   snprintf(ctemp , 100, "Number of pairs which passed:\t%ld  Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
 165   stemp += ctemp;
 166   AliFemtoString returnThis = stemp;
 167   return returnThis;}
 168 //__________________
 169
 170 TList *AliFemtoPairCutRadialDistance::ListSettings()
 171 {
 172   // return a list of settings in a writable form
 173   TList *tListSetttings =  AliFemtoPairCutAntiGamma::ListSettings();
 174   char buf[200];
 175   snprintf(buf, 200, "AliFemtoPairCutRadialDistance.phistarsepmin=%f", fDPhiStarMin);
 176   tListSetttings->AddLast(new TObjString(buf));
 177
 178   return tListSetttings;
 179 }
 180
 181 void AliFemtoPairCutRadialDistance::SetPhiStarDifferenceMinimum(double dtpc)
 182 {
 183   fDPhiStarMin = dtpc;
 184 }
 185
 186 void AliFemtoPairCutRadialDistance::SetEtaDifferenceMinimum(double etpc)
 187 {
 188   fEtaMin = etpc;
 189 }
 190
 191
 192 void AliFemtoPairCutRadialDistance::SetMinimumRadius(double minrad)
 193 {
 194   fMinRad = minrad;
 195 }
 196
 197 void AliFemtoPairCutRadialDistance::SetMaximumRadius(double maxrad)
 198 {
 199   fMaxRad = maxrad;
 200 }
 201
 202 void AliFemtoPairCutRadialDistance::SetMagneticFieldSign(int magsign)
 203 {
 204   if(magsign>1) fMagSign = 1;
 205   else if(magsign<1) fMagSign = -1;
 206   else fMagSign = magsign;
 207 }
 208
 209 void AliFemtoPairCutRadialDistance::SetPhiStarMin(Bool_t phistarmin)
 210 {
 211   fPhistarmin = phistarmin;
 212 }