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