PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutRadialDistanceKK.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 "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.1*Bfield)*0.15/ptv1);// 0.15 for D=1m
 106     if(phistar1 > 2*PI) phistar1 -= 2*PI;
 107     if(phistar1 < 0) phistar1 += 2*PI;
 108     Double_t phistar2 = phi2 - asin(chg2*(0.1*Bfield)*0.15/ptv2);// 0.15 for D=1m
 109     if(phistar2 > 2*PI) phistar2 -= 2*PI;
 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.1*Bfield)*0.24/ptv1);// mine. 0.24 for D=1.6m
 130     if(phistar1 > 2*PI) phistar1 -= 2*PI;
 131     if(phistar1 < 0) phistar1 += 2*PI;
 132     phistar2 = phi2 - asin(chg2*(0.1*Bfield)*0.24/ptv2);// mine. 0.24 for D=1.6m
 133     if(phistar2 > 2*PI) phistar2 -= 2*PI;
 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 }