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   double fVCov[6];
 236   if (fUseTPCOnly) {
 237     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 238     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 239     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
 240       fVCov[4] = -1001.0;
 241   }
 242   else {
 243     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 244     fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 245     if (!fEvent->GetPrimaryVertex()->GetStatus())
 246       fVCov[4] = -1001.0;
 247   }
 248
 249   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 250   hbtEvent->SetPrimVertPos(vertex);
 251   hbtEvent->SetPrimVertCov(fVCov);
 252
 253   //starting to reading tracks
 254   int nofTracks=0;  //number of reconstructed tracks in event
 255   nofTracks=fEvent->GetNumberOfTracks();
 256   int realnofTracks=0;//number of track which we use ina analysis
 257
 258 //   // Clear the shared cluster list
 259 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 260 //     fClusterPerPadrow[tPad]->clear();
 261 //   }
 262 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 263 //     fSharedList[tPad]->clear();
 264 //   }
 265
 266
 267 //   for (int i=0;i<nofTracks;i++) {
 268 //     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 269
 270 //     list<Int_t>::iterator tClustIter;
 271
 272 //     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 273 //     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 274 //     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 275 //       if (tTrackIndices[tNcl] >= 0) {
 276 //      tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 277 //      if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 278 //        fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 279 //      }
 280 //      else {
 281 //        fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 282 //      }
 283 //       }
 284 //     }
 285
 286 //   }
 287
 288   for (int i=0;i<nofTracks;i++)
 289     {
 290       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 291
 292       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 293       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 294       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 295
 296       trackCopy->SetCharge((short)esdtrack->GetSign());
 297
 298       //in aliroot we have AliPID
 299       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 300       //we use only 5 first
 301       double esdpid[5];
 302       esdtrack->GetESDpid(esdpid);
 303       trackCopy->SetPidProbElectron(esdpid[0]);
 304       trackCopy->SetPidProbMuon(esdpid[1]);
 305       trackCopy->SetPidProbPion(esdpid[2]);
 306       trackCopy->SetPidProbKaon(esdpid[3]);
 307       trackCopy->SetPidProbProton(esdpid[4]);
 308
 309       double pxyz[3];
 310       double rxyz[3];
 311       double impact[2];
 312       double covimpact[3];
 313
 314       if (fUseTPCOnly) {
 315         if (!esdtrack->GetTPCInnerParam()) {
 316           delete trackCopy;
 317           continue;
 318         }
 319
 320
 321         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 322         param->GetXYZ(rxyz);
 323         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 324         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 325
 326         if (fReadInner == true) {
 327           AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 328           tInfo->SetPDGPid(211);
 329           tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 330           tInfo->SetMass(0.13957);
 331           //      tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 332           //      tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 333           tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 334           trackCopy->SetHiddenInfo(tInfo);
 335         }
 336
 337         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 338         if (v.mag() < 0.0001) {
 339           //    cout << "Found 0 momentum ???? " <<endl;
 340           delete trackCopy;
 341           continue;
 342         }
 343         trackCopy->SetP(v);//setting momentum
 344         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 345
 346         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 347         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 348         //setting helix I do not if it is ok
 349         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 350         trackCopy->SetHelix(helix);
 351
 352         //some stuff which could be useful
 353         trackCopy->SetImpactD(impact[0]);
 354         trackCopy->SetImpactZ(impact[1]);
 355         trackCopy->SetCdd(covimpact[0]);
 356         trackCopy->SetCdz(covimpact[1]);
 357         trackCopy->SetCzz(covimpact[2]);
 358         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 359
 360         delete param;
 361       }
 362       else {
 363         if (fReadInner == true) {
 364
 365           if (esdtrack->GetTPCInnerParam()) {
 366             AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 367             param->GetXYZ(rxyz);
 368             param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 369             param->GetPxPyPz(pxyz);//reading noconstarined momentum
 370             delete param;
 371
 372             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 373             tInfo->SetPDGPid(211);
 374             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 375             tInfo->SetMass(0.13957);
 376             //      tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 377             //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 378             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 379             trackCopy->SetHiddenInfo(tInfo);
 380           }
 381         }
 382         if (fConstrained==true)
 383           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 384         else
 385           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 386
 387         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 388         if (v.mag() < 0.0001) {
 389           //    cout << "Found 0 momentum ???? " <<endl;
 390           delete trackCopy;
 391           continue;
 392         }
 393         trackCopy->SetP(v);//setting momentum
 394         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 395         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 396         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 397         //setting helix I do not if it is ok
 398         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 399         trackCopy->SetHelix(helix);
 400
 401         //some stuff which could be useful
 402         float imp[2];
 403         float cim[3];
 404         esdtrack->GetImpactParameters(imp,cim);
 405
 406         impact[0] = imp[0];
 407         impact[1] = imp[1];
 408         covimpact[0] = cim[0];
 409         covimpact[1] = cim[1];
 410         covimpact[2] = cim[2];
 411
 412         trackCopy->SetImpactD(impact[0]);
 413         trackCopy->SetImpactZ(impact[1]);
 414         trackCopy->SetCdd(covimpact[0]);
 415         trackCopy->SetCdz(covimpact[1]);
 416         trackCopy->SetCzz(covimpact[2]);
 417         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 418       }
 419
 420       trackCopy->SetTrackId(esdtrack->GetID());
 421       trackCopy->SetFlags(esdtrack->GetStatus());
 422       trackCopy->SetLabel(esdtrack->GetLabel());
 423
 424       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 425       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 426       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 427       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 428       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 429       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 430       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 431
 432       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 433       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 434
 435       double xtpc[3];
 436       esdtrack->GetInnerXYZ(xtpc);
 437       xtpc[2] -= fV1[2];
 438       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 439
 440       esdtrack->GetOuterXYZ(xtpc);
 441       xtpc[2] -= fV1[2];
 442       trackCopy->SetNominalTPCExitPoint(xtpc);
 443
 444       int indexes[3];
 445       for (int ik=0; ik<3; ik++) {
 446         indexes[ik] = esdtrack->GetKinkIndex(ik);
 447       }
 448       trackCopy->SetKinkIndexes(indexes);
 449       //decision if we want this track
 450       //if we using diffrent labels we want that this label was use for first time
 451       //if we use hidden info we want to have match between sim data and ESD
 452       if (tGoodMomentum==true)
 453         {
 454           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 455           realnofTracks++;//real number of tracks
 456           //      delete trackCopy;
 457         }
 458       else
 459         {
 460           delete  trackCopy;
 461         }
 462
 463     }
 464
 465   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 466   fCurEvent++;
 467   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 468   return hbtEvent;
 469 }
 470 //___________________
 471 void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
 472 {
 473   // The chain loads the ESD for us
 474   // You must provide the address where it can be found
 475   fEvent = aESD;
 476 }
 477 //___________________
 478 // void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 479 // {
 480 //   // We need the ESD tree to obtain
 481 //   // information about the friend file location
 482 //   fEventFriend = aFriend;
 483 // }
 484
 485 //____________________________________________________________________
 486 Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
 487 {
 488   // Calculates the number of sigma to the vertex.
 489
 490   Float_t b[2];
 491   Float_t bRes[2];
 492   Float_t bCov[3];
 493
 494   b[0] = impact[0];
 495   b[1] = impact[1];
 496   bCov[0] = covar[0];
 497   bCov[1] = covar[1];
 498   bCov[2] = covar[2];
 499
 500   bRes[0] = TMath::Sqrt(bCov[0]);
 501   bRes[1] = TMath::Sqrt(bCov[2]);
 502
 503   // -----------------------------------
 504   // How to get to a n-sigma cut?
 505   //
 506   // The accumulated statistics from 0 to d is
 507   //
 508   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 509   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 510   //
 511   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 512   // Can this be expressed in a different way?
 513
 514   if (bRes[0] == 0 || bRes[1] ==0)
 515     return -1;
 516
 517   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 518
 519   // stupid rounding problem screws up everything:
 520   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 521   if (TMath::Exp(-d * d / 2) < 1e-10)
 522     return 1000;
 523
 524   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 525   return d;
 526 }
 527
 528
 529
 530
 531
 532
 533
 534