PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChain.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoEventReaderESDChain - the reader class for the Alice ESD         ///
   4 /// tailored for the Task framework                                 ///
   5 /// Reads in AliESDfriend to create shared hit/quality information           ///
   6 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
   7 ///                                                                          ///
   8 ////////////////////////////////////////////////////////////////////////////////
   9 #include "AliFemtoEventReaderESDChain.h"
  10
  11 #include "TFile.h"
  12 #include "TTree.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDtrack.h"
  15 #include "AliESDVertex.h"
  16
  17 #include "AliFmPhysicalHelixD.h"
  18 #include "AliFmThreeVectorF.h"
  19
  20 #include "SystemOfUnits.h"
  21
  22 #include "AliFemtoEvent.h"
  23 #include "AliFemtoModelHiddenInfo.h"
  24
  25 ClassImp(AliFemtoEventReaderESDChain)
  26
  27 #if !(ST_NO_NAMESPACES)
  28   using namespace units;
  29 #endif
  30
  31 using namespace std;
  32 //____________________________
  33 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  34   fFileName(" "),
  35   fConstrained(true),
  36   fReadInner(false),
  37   fUseTPCOnly(false),
  38   fNumberofEvent(0),
  39   fCurEvent(0),
  40   fCurFile(0),
  41   fEvent(0x0)
  42 {
  43   //constructor with 0 parameters , look at default settings
  44 //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  45 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  46 //     fClusterPerPadrow[tPad] = new list<Int_t>();
  47 //   }
  48 //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  49 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  50 //     fSharedList[tPad] = new list<Int_t>();
  51 //   }
  52 }
  53
  54 //__________________
  55 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  56   AliFemtoEventReader(aReader),
  57   fFileName(" "),
  58   fConstrained(true),
  59   fReadInner(false),
  60   fUseTPCOnly(false),
  61   fNumberofEvent(0),
  62   fCurEvent(0),
  63   fCurFile(0),
  64   fEvent(0x0)
  65 {
  66   // Copy constructor
  67   fConstrained = aReader.fConstrained;
  68   fReadInner = aReader.fReadInner;
  69   fUseTPCOnly = aReader.fUseTPCOnly;
  70   fNumberofEvent = aReader.fNumberofEvent;
  71   fCurEvent = aReader.fCurEvent;
  72   fCurFile = aReader.fCurFile;
  73   //  fEvent = new AliESD(*aReader.fEvent);
  74   fEvent = new AliESDEvent();
  75 //   fEventFriend = aReader.fEventFriend;
  76 //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  77 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  78 //     fClusterPerPadrow[tPad] = new list<Int_t>();
  79 //     list<Int_t>::iterator iter;
  80 //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
  81 //       fClusterPerPadrow[tPad]->push_back(*iter);
  82 //     }
  83 //   }
  84 //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  85 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  86 //     fSharedList[tPad] = new list<Int_t>();
  87 //     list<Int_t>::iterator iter;
  88 //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
  89 //       fSharedList[tPad]->push_back(*iter);
  90 //     }
  91 //   }
  92 }
  93 //__________________
  94 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
  95 {
  96   //Destructor
  97   delete fEvent;
  98
  99 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 100 //     fClusterPerPadrow[tPad]->clear();
 101 //     delete fClusterPerPadrow[tPad];
 102 //   }
 103 //   delete [] fClusterPerPadrow;
 104 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 105 //     fSharedList[tPad]->clear();
 106 //     delete fSharedList[tPad];
 107 //   }
 108 //   delete [] fSharedList;
 109 }
 110
 111 //__________________
 112 AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
 113 {
 114   // Assignment operator
 115   if (this == &aReader)
 116     return *this;
 117
 118   fConstrained = aReader.fConstrained;
 119   fReadInner = aReader.fReadInner;
 120   fUseTPCOnly = aReader.fUseTPCOnly;
 121   fNumberofEvent = aReader.fNumberofEvent;
 122   fCurEvent = aReader.fCurEvent;
 123   fCurFile = aReader.fCurFile;
 124   if (fEvent) delete fEvent;
 125   fEvent = new AliESDEvent();
 126
 127   //  fEventFriend = aReader.fEventFriend;
 128
 129 //   if (fClusterPerPadrow) {
 130 //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 131 //       fClusterPerPadrow[tPad]->clear();
 132 //       delete fClusterPerPadrow[tPad];
 133 //     }
 134 //     delete [] fClusterPerPadrow;
 135 //   }
 136
 137 //   if (fSharedList) {
 138 //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 139 //       fSharedList[tPad]->clear();
 140 //       delete fSharedList[tPad];
 141 //     }
 142 //     delete [] fSharedList;
 143 //   }
 144
 145 //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 146 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 147 //     fClusterPerPadrow[tPad] = new list<Int_t>();
 148 //     list<Int_t>::iterator iter;
 149 //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 150 //       fClusterPerPadrow[tPad]->push_back(*iter);
 151 //     }
 152 //   }
 153 //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 154 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 155 //     fSharedList[tPad] = new list<Int_t>();
 156 //     list<Int_t>::iterator iter;
 157 //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 158 //       fSharedList[tPad]->push_back(*iter);
 159 //     }
 160 //   }
 161
 162   return *this;
 163 }
 164 //__________________
 165 // Simple report
 166 AliFemtoString AliFemtoEventReaderESDChain::Report()
 167 {
 168   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
 169   return temp;
 170 }
 171
 172 //__________________
 173 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
 174 {
 175   // Select whether to read constrained or not constrained momentum
 176   fConstrained=constrained;
 177 }
 178 //__________________
 179 bool AliFemtoEventReaderESDChain::GetConstrained() const
 180 {
 181   // Check whether we read constrained or not constrained momentum
 182   return fConstrained;
 183 }
 184 //__________________
 185 void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
 186 {
 187   fReadInner=readinner;
 188 }
 189
 190 bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
 191 {
 192   return fReadInner;
 193 }
 194
 195 //__________________
 196 void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
 197 {
 198   fUseTPCOnly=usetpconly;
 199 }
 200
 201 bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
 202 {
 203   return fUseTPCOnly;
 204 }
 205
 206 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
 207 {
 208   // Get the event, read all the relevant information
 209   // and fill the AliFemtoEvent class
 210   // Returns a valid AliFemtoEvent
 211   AliFemtoEvent *hbtEvent = 0;
 212   string tFriendFileName;
 213
 214   // Get the friend information
 215   cout<<"starting to read event "<<fCurEvent<<endl;
 216   //  fEvent->SetESDfriend(fEventFriend);
 217
 218   hbtEvent = new AliFemtoEvent;
 219   //setting basic things
 220   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 221   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 222   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 223   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 224   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 225   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 226   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 227   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 228   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 229   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 230   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 231   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 232
 233   //Vertex
 234   double fV1[3];
 235   if (fUseTPCOnly) {
 236     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 237   }
 238   else {
 239     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 240   }
 241
 242   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 243   hbtEvent->SetPrimVertPos(vertex);
 244
 245   //starting to reading tracks
 246   int nofTracks=0;  //number of reconstructed tracks in event
 247   nofTracks=fEvent->GetNumberOfTracks();
 248   int realnofTracks=0;//number of track which we use ina analysis
 249
 250 //   // Clear the shared cluster list
 251 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 252 //     fClusterPerPadrow[tPad]->clear();
 253 //   }
 254 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 255 //     fSharedList[tPad]->clear();
 256 //   }
 257
 258
 259 //   for (int i=0;i<nofTracks;i++) {
 260 //     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 261
 262 //     list<Int_t>::iterator tClustIter;
 263
 264 //     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 265 //     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 266 //     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 267 //       if (tTrackIndices[tNcl] >= 0) {
 268 //      tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 269 //      if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 270 //        fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 271 //      }
 272 //      else {
 273 //        fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 274 //      }
 275 //       }
 276 //     }
 277
 278 //   }
 279
 280   for (int i=0;i<nofTracks;i++)
 281     {
 282       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 283
 284       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 285       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 286       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 287
 288       trackCopy->SetCharge((short)esdtrack->GetSign());
 289
 290       //in aliroot we have AliPID
 291       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 292       //we use only 5 first
 293       double esdpid[5];
 294       esdtrack->GetESDpid(esdpid);
 295       trackCopy->SetPidProbElectron(esdpid[0]);
 296       trackCopy->SetPidProbMuon(esdpid[1]);
 297       trackCopy->SetPidProbPion(esdpid[2]);
 298       trackCopy->SetPidProbKaon(esdpid[3]);
 299       trackCopy->SetPidProbProton(esdpid[4]);
 300
 301       double pxyz[3];
 302       double rxyz[3];
 303       double impact[2];
 304       double covimpact[3];
 305
 306       if (fUseTPCOnly) {
 307         if (!esdtrack->GetTPCInnerParam()) {
 308           delete trackCopy;
 309           continue;
 310         }
 311
 312
 313         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 314         param->GetXYZ(rxyz);
 315         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 316         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 317
 318         if (fReadInner == true) {
 319           AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 320           tInfo->SetPDGPid(211);
 321           tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 322           tInfo->SetMass(0.13957);
 323           //      tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 324           //      tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 325           tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 326           trackCopy->SetHiddenInfo(tInfo);
 327         }
 328
 329         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 330         if (v.mag() < 0.0001) {
 331           //    cout << "Found 0 momentum ???? " <<endl;
 332           delete trackCopy;
 333           continue;
 334         }
 335         trackCopy->SetP(v);//setting momentum
 336         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 337
 338         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 339         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 340         //setting helix I do not if it is ok
 341         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 342         trackCopy->SetHelix(helix);
 343
 344         //some stuff which could be useful
 345         trackCopy->SetImpactD(impact[0]);
 346         trackCopy->SetImpactZ(impact[1]);
 347         trackCopy->SetCdd(covimpact[0]);
 348         trackCopy->SetCdz(covimpact[1]);
 349         trackCopy->SetCzz(covimpact[2]);
 350         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 351
 352         delete param;
 353       }
 354       else {
 355         if (fReadInner == true) {
 356
 357           if (esdtrack->GetTPCInnerParam()) {
 358             AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 359             param->GetXYZ(rxyz);
 360             param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 361             param->GetPxPyPz(pxyz);//reading noconstarined momentum
 362             delete param;
 363
 364             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 365             tInfo->SetPDGPid(211);
 366             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 367             tInfo->SetMass(0.13957);
 368             //      tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 369             //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 370             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 371             trackCopy->SetHiddenInfo(tInfo);
 372           }
 373         }
 374         if (fConstrained==true)
 375           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 376         else
 377           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 378
 379         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 380         if (v.mag() < 0.0001) {
 381           //    cout << "Found 0 momentum ???? " <<endl;
 382           delete trackCopy;
 383           continue;
 384         }
 385         trackCopy->SetP(v);//setting momentum
 386         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 387         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 388         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 389         //setting helix I do not if it is ok
 390         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 391         trackCopy->SetHelix(helix);
 392
 393         //some stuff which could be useful
 394         float imp[2];
 395         float cim[3];
 396         esdtrack->GetImpactParameters(imp,cim);
 397
 398         impact[0] = imp[0];
 399         impact[1] = imp[1];
 400         covimpact[0] = cim[0];
 401         covimpact[1] = cim[1];
 402         covimpact[2] = cim[2];
 403
 404         trackCopy->SetImpactD(impact[0]);
 405         trackCopy->SetImpactZ(impact[1]);
 406         trackCopy->SetCdd(covimpact[0]);
 407         trackCopy->SetCdz(covimpact[1]);
 408         trackCopy->SetCzz(covimpact[2]);
 409         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 410       }
 411
 412       trackCopy->SetTrackId(esdtrack->GetID());
 413       trackCopy->SetFlags(esdtrack->GetStatus());
 414       trackCopy->SetLabel(esdtrack->GetLabel());
 415
 416       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 417       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 418       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 419       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 420       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 421       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 422       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 423
 424       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 425       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 426
 427       double xtpc[3];
 428       esdtrack->GetInnerXYZ(xtpc);
 429       xtpc[2] -= fV1[2];
 430       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 431
 432       esdtrack->GetOuterXYZ(xtpc);
 433       xtpc[2] -= fV1[2];
 434       trackCopy->SetNominalTPCExitPoint(xtpc);
 435
 436       int indexes[3];
 437       for (int ik=0; ik<3; ik++) {
 438         indexes[ik] = esdtrack->GetKinkIndex(ik);
 439       }
 440       trackCopy->SetKinkIndexes(indexes);
 441       //decision if we want this track
 442       //if we using diffrent labels we want that this label was use for first time
 443       //if we use hidden info we want to have match between sim data and ESD
 444       if (tGoodMomentum==true)
 445         {
 446           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 447           realnofTracks++;//real number of tracks
 448           //      delete trackCopy;
 449         }
 450       else
 451         {
 452           delete  trackCopy;
 453         }
 454
 455     }
 456
 457   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 458   fCurEvent++;
 459   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 460   return hbtEvent;
 461 }
 462 //___________________
 463 void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
 464 {
 465   // The chain loads the ESD for us
 466   // You must provide the address where it can be found
 467   fEvent = aESD;
 468 }
 469 //___________________
 470 // void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 471 // {
 472 //   // We need the ESD tree to obtain
 473 //   // information about the friend file location
 474 //   fEventFriend = aFriend;
 475 // }
 476
 477 //____________________________________________________________________
 478 Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
 479 {
 480   // Calculates the number of sigma to the vertex.
 481
 482   Float_t b[2];
 483   Float_t bRes[2];
 484   Float_t bCov[3];
 485
 486   b[0] = impact[0];
 487   b[1] = impact[1];
 488   bCov[0] = covar[0];
 489   bCov[1] = covar[1];
 490   bCov[2] = covar[2];
 491
 492   bRes[0] = TMath::Sqrt(bCov[0]);
 493   bRes[1] = TMath::Sqrt(bCov[2]);
 494
 495   // -----------------------------------
 496   // How to get to a n-sigma cut?
 497   //
 498   // The accumulated statistics from 0 to d is
 499   //
 500   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 501   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 502   //
 503   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 504   // Can this be expressed in a different way?
 505
 506   if (bRes[0] == 0 || bRes[1] ==0)
 507     return -1;
 508
 509   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 510
 511   // stupid rounding problem screws up everything:
 512   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 513   if (TMath::Exp(-d * d / 2) < 1e-10)
 514     return 1000;
 515
 516   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 517   return d;
 518 }
 519
 520
 521
 522
 523
 524
 525
 526