PWGCF/FEMTOSCOPY/AliFemto/AliFemtoV0TrackPairCut.cxx

   1 /////////////////////////////////////////////////////////////////////////////
   2 //                                                                         //
   3 // AliFemtoShareQualityPairCut - a pair cut which checks for some pair     //
   4 // qualities that attempt to identify slit/doubly reconstructed tracks     //
   5 //                                                                         //
   6 /////////////////////////////////////////////////////////////////////////////
   7 /***************************************************************************
   8  *
   9  * $Id: AliFemtoShareQualityPairCut.cxx 50722 2011-07-21 15:18:38Z akisiel $
  10  *
  11  * Author: Adam Kisiel, Ohio State, kisiel@mps.ohio-state.edu
  12  ***************************************************************************
  13  *
  14  * Description: part of STAR HBT Framework: AliFemtoMaker package
  15  *   a cut to remove "shared" and "split" pairs
  16  *
  17  ***************************************************************************
  18  *
  19  *
  20  **************************************************************************/
  21
  22 #include "AliFemtoV0TrackPairCut.h"
  23 #include <string>
  24 #include <cstdio>
  25
  26 #ifdef __ROOT__
  27 ClassImp(AliFemtoV0TrackPairCut)
  28 #endif
  29
  30 //__________________
  31 AliFemtoV0TrackPairCut::AliFemtoV0TrackPairCut():
  32   fNPairsPassed(0),
  33   fNPairsFailed(0),
  34   fV0Max(1.0),
  35   fShareQualityMax(1.0),
  36   fShareFractionMax(1.0),
  37   fRemoveSameLabel(0),
  38   fTrackTPCOnly(0),
  39   fDataType(kAOD),
  40   fDTPCMin(0),
  41   fDTPCExitMin(0),
  42   fKstarCut(0),
  43   fFirstParticleType(kLambda),
  44   fSecondParticleType(kProton),
  45   fMinAvgSepTrackPos(0),
  46   fMinAvgSepTrackNeg(0)
  47 {
  48   // Default constructor
  49   // Nothing to do
  50 }
  51 //__________________
  52 AliFemtoV0TrackPairCut::~AliFemtoV0TrackPairCut(){
  53   /* no-op */
  54 }
  55 AliFemtoV0TrackPairCut& AliFemtoV0TrackPairCut::operator=(const AliFemtoV0TrackPairCut& cut)
  56 {
  57   if (this != &cut) {
  58
  59     AliFemtoPairCut::operator=(cut);
  60     fNPairsPassed = 0;
  61     fNPairsFailed =0;
  62     fV0Max = 1.0;
  63     fShareQualityMax = 1.0;
  64     fShareFractionMax = 1.0;
  65     fRemoveSameLabel = 0;
  66     fTrackTPCOnly = 0;
  67     fDataType = kAOD;
  68     fDTPCMin = 0;
  69     fDTPCExitMin = 0;
  70     fMinAvgSepTrackPos = 0;
  71     fMinAvgSepTrackNeg = 0;
  72   }
  73
  74   return *this;
  75 }
  76 //__________________
  77 bool AliFemtoV0TrackPairCut::Pass(const AliFemtoPair* pair){
  78   // Check for pairs that are possibly shared/double reconstruction
  79
  80   bool temp = true;
  81   //Track1 - V0
  82   //Track2 - track
  83
  84   if(!(pair->Track1()->V0() && pair->Track2()->Track()))
  85     {
  86       return false;
  87     }
  88   if(fTrackTPCOnly)
  89     {
  90       if( (-(pair->Track1()->V0()->IdNeg()+1)) ==pair->Track2()->TrackId() || (-(pair->Track1()->V0()->IdPos()+1)) ==pair->Track2()->TrackId())
  91         {
  92           return false;
  93         }
  94     }
  95   else
  96     {
  97       if(pair->Track1()->V0()->IdNeg()==pair->Track2()->TrackId() || pair->Track1()->V0()->IdPos()==pair->Track2()->TrackId())
  98         {
  99           return false;
 100         }
 101     }
 102
 103   bool tempTPCEntrancePos = true;
 104   bool tempTPCEntranceNeg = true;
 105   bool tempTPCExitPos = true;
 106   bool tempTPCExitNeg = true;
 107   if(fDataType==kESD || fDataType==kAOD)
 108     {
 109       double distx = pair->Track1()->V0()->NominalTpcEntrancePointPos().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 110       double disty = pair->Track1()->V0()->NominalTpcEntrancePointPos().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 111       double distz = pair->Track1()->V0()->NominalTpcEntrancePointPos().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 112       double distPos = sqrt(distx*distx + disty*disty + distz*distz);
 113
 114       distx = pair->Track1()->V0()->NominalTpcEntrancePointNeg().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 115       disty = pair->Track1()->V0()->NominalTpcEntrancePointNeg().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 116       distz = pair->Track1()->V0()->NominalTpcEntrancePointNeg().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 117       double distNeg = sqrt(distx*distx + disty*disty + distz*distz);
 118
 119       double distExitX = pair->Track1()->V0()->NominalTpcExitPointPos().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
 120       double distExitY = pair->Track1()->V0()->NominalTpcExitPointPos().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
 121       double distExitZ = pair->Track1()->V0()->NominalTpcExitPointPos().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
 122       double distExitPos = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);
 123
 124       distExitX = pair->Track1()->V0()->NominalTpcExitPointNeg().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
 125       distExitY = pair->Track1()->V0()->NominalTpcExitPointNeg().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
 126       distExitZ = pair->Track1()->V0()->NominalTpcExitPointNeg().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
 127       double distExitNeg = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);
 128
 129       tempTPCEntrancePos = distPos > fDTPCMin;
 130       tempTPCEntranceNeg = distNeg > fDTPCMin;
 131       tempTPCExitPos = distExitPos > fDTPCExitMin;
 132       tempTPCExitNeg = distExitNeg > fDTPCExitMin;
 133     }
 134
 135   if(!tempTPCEntrancePos || !tempTPCEntranceNeg || !tempTPCExitPos || !tempTPCExitNeg) return false;
 136
 137   //reject merged trakcs in TPC
 138
 139
 140   //temp = dist > fDTPCMin;
 141   //koniec kopii
 142   //if(!temp)
 143   //return false;
 144
 145   //kopia z AliFemtoShareQualityPairCut.cxx
 146   Int_t nh = 0;
 147   Int_t an = 0;
 148   Int_t ns = 0;
 149
 150
 151
 152   if ((fShareFractionMax < 1.0) && ( fShareQualityMax < 1.0)) {
 153     for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersPos().GetNbits(); imap++) {
 154       // If both have clusters in the same row
 155       if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) &&
 156           pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 157         // Do they share it ?
 158         if (pair->Track1()->V0()->TPCsharingPos().TestBitNumber(imap) &&
 159             pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
 160           {
 161             //    cout << "A shared cluster !!!" << endl;
 162             //  cout << "imap idx1 idx2 "
 163             //       << imap << " "
 164             //       << tP1idx[imap] << " " << tP2idx[imap] << endl;
 165             an++;
 166             nh+=2;
 167             ns+=2;
 168           }
 169
 170         // Different hits on the same padrow
 171         else {
 172           an--;
 173           nh+=2;
 174         }
 175       }
 176       else if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) ||
 177                pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 178         // One track has a hit, the other does not
 179         an++;
 180         nh++;
 181       }
 182     }
 183
 184     Float_t hsmval = 0.0;
 185     Float_t hsfval = 0.0;
 186
 187     if (nh >0) {
 188       hsmval = an*1.0/nh;
 189       hsfval = ns*1.0/nh;
 190     }
 191     //  if (hsmval > -0.4) {
 192     //   cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
 193     //        << (pair->Track1()->Track()) << " "
 194     //        << (pair->Track2()->Track()) << endl;
 195     //   cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
 196     //  }
 197     //   if (hsfval > 0.0) {
 198     //     cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << "    " << hsmval << " " << an << " " << nh << endl;
 199     //   }
 200
 201     temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
 202     if(!temp) return false;
 203
 204     nh = 0;
 205     an = 0;
 206     ns = 0;
 207
 208     for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersNeg().GetNbits(); imap++) {
 209       // If both have clusters in the same row
 210       if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) &&
 211           pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 212         // Do they share it ?
 213         if (pair->Track1()->V0()->TPCsharingNeg().TestBitNumber(imap) &&
 214             pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
 215           {
 216             //    cout << "A shared cluster !!!" << endl;
 217             //  cout << "imap idx1 idx2 "
 218             //       << imap << " "
 219             //       << tP1idx[imap] << " " << tP2idx[imap] << endl;
 220             an++;
 221             nh+=2;
 222             ns+=2;
 223           }
 224
 225         // Different hits on the same padrow
 226         else {
 227           an--;
 228           nh+=2;
 229         }
 230       }
 231       else if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) ||
 232                pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 233         // One track has a hit, the other does not
 234         an++;
 235         nh++;
 236       }
 237     }
 238
 239     hsmval = 0.0;
 240     hsfval = 0.0;
 241
 242     if (nh >0) {
 243       hsmval = an*1.0/nh;
 244       hsfval = ns*1.0/nh;
 245     }
 246     //  if (hsmval > -0.4) {
 247     //   cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
 248     //        << (pair->Track1()->Track()) << " "
 249     //        << (pair->Track2()->Track()) << endl;
 250     //   cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
 251     //  }
 252     //   if (hsfval > 0.0) {
 253     //     cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << "    " << hsmval << " " << an << " " << nh << endl;
 254     //   }
 255
 256     temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
 257
 258
 259   }
 260   else
 261     temp = true;
 262   //koniec kopii
 263
 264
 265   //avg sep pair cut
 266   double avgSep=0;
 267   AliFemtoThreeVector first, second, tmp;
 268   for(int i=0; i<8 ;i++)
 269     {
 270       tmp = pair->Track1()->V0()->NominalTpcPointPos(i);
 271       //cout<<"X pos: "<<tmp.x()<<endl;
 272       first.SetX((double)(tmp.x()));
 273       first.SetY((double)tmp.y());
 274       first.SetZ((double)tmp.z());
 275
 276       tmp = pair->Track2()->Track()->NominalTpcPoint(i);
 277       second.SetX((double)tmp.x());
 278       second.SetY((double)tmp.y());
 279       second.SetZ((double)tmp.z());
 280
 281       avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())*((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())*((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
 282     }
 283   avgSep /= 8;
 284
 285   if(avgSep<fMinAvgSepTrackPos) return false;
 286
 287   avgSep = 0;
 288
 289   for(int i=0; i<8 ;i++)
 290     {
 291       tmp = pair->Track1()->V0()->NominalTpcPointNeg(i);
 292       //cout<<"X pos: "<<tmp.x()<<endl;
 293       first.SetX((double)(tmp.x()));
 294       first.SetY((double)tmp.y());
 295       first.SetZ((double)tmp.z());
 296
 297       tmp = pair->Track2()->Track()->NominalTpcPoint(i);
 298       second.SetX((double)tmp.x());
 299       second.SetY((double)tmp.y());
 300       second.SetZ((double)tmp.z());
 301
 302       avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())*((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())*((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
 303     }
 304   avgSep /= 8;
 305
 306   if(avgSep<fMinAvgSepTrackNeg) return false;
 307
 308
 309
 310
 311   //Qinv cut (we are trying to get rid of antiresidual correlation between primary protons)
 312   if(fKstarCut > 0)
 313     {
 314
 315       //2 daughters of first V0
 316       double mom1PosX = pair->Track1()->V0()->MomPosX();
 317       double mom1PosY = pair->Track1()->V0()->MomPosY();
 318       double mom1PosZ = pair->Track1()->V0()->MomPosZ();
 319       double mom1NegX = pair->Track1()->V0()->MomNegX();
 320       double mom1NegY = pair->Track1()->V0()->MomNegY();
 321       double mom1NegZ = pair->Track1()->V0()->MomNegZ();
 322
 323
 324       //double PionMass = 0.13956995;
 325       //double KaonMass = 0.493677;
 326       double ProtonMass = 0.938272;
 327       //double LambdaMass = 1.115683;
 328
 329       AliFemtoLorentzVector fFourMomentum1; // Particle momentum
 330       AliFemtoLorentzVector fFourMomentum2; // Particle momentum
 331
 332       AliFemtoThreeVector temp1;
 333       double ener1=0;
 334       if(fFirstParticleType == 0) //lambda
 335         {
 336           if(fSecondParticleType == 2) //proton
 337             {
 338               //temp1 = ::sqrt(mom1PosX*mom1PosX+mom1PosY*mom1PosY+mom1PosZ*mom1PosZ);
 339               temp1.SetX(mom1PosX);
 340               temp1.SetY(mom1PosY);
 341               temp1.SetZ(mom1PosZ);
 342               ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
 343             }
 344         }
 345       else if(fFirstParticleType == 1) //antilambda
 346         {
 347           if(fSecondParticleType == 3) //antiproton
 348             {
 349               //temp1 = ::sqrt(mom1NegX*mom1NegX+mom1NegY*mom1NegY+mom1NegZ*mom1NegZ);
 350               temp1.SetX(mom1NegX);
 351               temp1.SetY(mom1NegY);
 352               temp1.SetZ(mom1NegZ);
 353               ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
 354             }
 355         }
 356       fFourMomentum1.SetVect(temp1);
 357       fFourMomentum1.SetE(ener1);
 358
 359       //AliFemtoLorentzVector fFourMomentum2; // Particle momentum
 360       AliFemtoThreeVector temp2;
 361       double ener2=0;
 362
 363       //2 daughters of second V0
 364       temp2 = pair->Track2()->Track()->P();
 365
 366
 367       if(fSecondParticleType == 2 || fSecondParticleType == 3) //proton
 368         {
 369             ener2 = ::sqrt(temp2.Mag2()+ProtonMass*ProtonMass);
 370         }
 371
 372       fFourMomentum2.SetVect(temp2);
 373       fFourMomentum2.SetE(ener2);
 374
 375       //QInv calculation
 376       AliFemtoLorentzVector tDiff = (fFourMomentum1-fFourMomentum2);
 377
 378       double tQinv = fabs(-1.*tDiff.m());   // note - qInv() will be negative for identical pairs...
 379
 380       //cout<<"tQinv/2: "<<tQinv/2<<endl;
 381       if(tQinv/2 < fKstarCut)
 382         {
 383
 384           temp = false;
 385         }
 386     }
 387   return temp;
 388 }
 389 //__________________
 390 AliFemtoString AliFemtoV0TrackPairCut::Report(){
 391   // Prepare the report from the execution
 392   string stemp = "AliFemtoV0 Pair Cut - remove shared and split pairs\n";  char ctemp[100];
 393   snprintf(ctemp , 100, "Number of pairs which passed:\t%ld  Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
 394   stemp += ctemp;
 395   AliFemtoString returnThis = stemp;
 396   return returnThis;}
 397 //__________________
 398
 399 void AliFemtoV0TrackPairCut::SetV0Max(Double_t aV0Max) {
 400   fV0Max = aV0Max;
 401 }
 402
 403 Double_t AliFemtoV0TrackPairCut::GetAliFemtoV0Max() const {
 404   return fV0Max;
 405 }
 406
 407
 408 TList *AliFemtoV0TrackPairCut::ListSettings()
 409 {
 410   // return a list of settings in a writable form
 411   TList *tListSetttings = new TList();
 412   char buf[200];
 413   snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharequalitymax=%f", fV0Max);
 414   snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharefractionmax=%f", fShareFractionMax);
 415   tListSetttings->AddLast(new TObjString(buf));
 416
 417   return tListSetttings;
 418 }
 419
 420 void     AliFemtoV0TrackPairCut::SetRemoveSameLabel(Bool_t aRemove)
 421 {
 422   fRemoveSameLabel = aRemove;
 423 }
 424
 425 void AliFemtoV0TrackPairCut::SetTPCOnly(Bool_t tpconly)
 426 {
 427   fTrackTPCOnly = tpconly;
 428 }
 429
 430 void AliFemtoV0TrackPairCut::SetShareQualityMax(Double_t aShareQualityMax) {
 431   fShareQualityMax = aShareQualityMax;
 432 }
 433
 434 void AliFemtoV0TrackPairCut::SetShareFractionMax(Double_t aShareFractionMax) {
 435   fShareFractionMax = aShareFractionMax;
 436 }
 437
 438 void AliFemtoV0TrackPairCut::SetDataType(AliFemtoDataType type)
 439 {
 440   fDataType = type;
 441 }
 442
 443 void AliFemtoV0TrackPairCut::SetTPCEntranceSepMinimum(double dtpc)
 444 {
 445   fDTPCMin = dtpc;
 446 }
 447
 448 void AliFemtoV0TrackPairCut::SetTPCExitSepMinimum(double dtpc)
 449 {
 450   fDTPCExitMin = dtpc;
 451 }
 452
 453 void AliFemtoV0TrackPairCut::SetKstarCut(double kstar, AliFemtoParticleType firstParticle, AliFemtoParticleType secondParticle)
 454 {
 455   fKstarCut = kstar;
 456   fFirstParticleType = firstParticle;  //for kstar - first particle type (V0 type)
 457   fSecondParticleType = secondParticle;
 458 }
 459
 460 void AliFemtoV0TrackPairCut::SetMinAvgSeparation(int type, double minSep)
 461 {
 462   if(type == 0) //Track-Pos
 463     fMinAvgSepTrackPos = minSep;
 464   else if(type == 1) //Track-Neg
 465     fMinAvgSepTrackNeg = minSep;
 466 }