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 {
  43   // Default constructor
  44   // Nothing to do
  45 }
  46 //__________________
  47 AliFemtoV0TrackPairCut::~AliFemtoV0TrackPairCut(){
  48   /* no-op */
  49 }
  50 AliFemtoV0TrackPairCut& AliFemtoV0TrackPairCut::operator=(const AliFemtoV0TrackPairCut& cut)
  51 {
  52   if (this != &cut) {
  53
  54     AliFemtoPairCut::operator=(cut);
  55     fNPairsPassed = 0;
  56     fNPairsFailed =0;
  57     fV0Max = 1.0;
  58     fShareQualityMax = 1.0;
  59     fShareFractionMax = 1.0;
  60     fRemoveSameLabel = 0;
  61     fTrackTPCOnly = 0;
  62     fDataType = kAOD;
  63     fDTPCMin = 0;
  64     fDTPCExitMin = 0;
  65   }
  66
  67   return *this;
  68 }
  69 //__________________
  70 bool AliFemtoV0TrackPairCut::Pass(const AliFemtoPair* pair){
  71   // Check for pairs that are possibly shared/double reconstruction
  72
  73   bool temp = true;
  74   //Track1 - V0
  75   //Track2 - track
  76
  77   if(!(pair->Track1()->V0() && pair->Track2()->Track()))
  78     {
  79       return false;
  80     }
  81   if(fTrackTPCOnly)
  82     {
  83       if( (-(pair->Track1()->V0()->IdNeg()+1)) ==pair->Track2()->TrackId() || (-(pair->Track1()->V0()->IdPos()+1)) ==pair->Track2()->TrackId())
  84         {
  85           return false;
  86         }
  87     }
  88   else
  89     {
  90       if(pair->Track1()->V0()->IdNeg()==pair->Track2()->TrackId() || pair->Track1()->V0()->IdPos()==pair->Track2()->TrackId())
  91         {
  92           return false;
  93         }
  94     }
  95
  96   bool tempTPCEntrancePos = true;
  97   bool tempTPCEntranceNeg = true;
  98   bool tempTPCExitPos = true;
  99   bool tempTPCExitNeg = true;
 100   if(fDataType==kESD || fDataType==kAOD)
 101     {
 102       double distx = pair->Track1()->V0()->NominalTpcEntrancePointPos().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 103       double disty = pair->Track1()->V0()->NominalTpcEntrancePointPos().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 104       double distz = pair->Track1()->V0()->NominalTpcEntrancePointPos().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 105       double distPos = sqrt(distx*distx + disty*disty + distz*distz);
 106
 107       distx = pair->Track1()->V0()->NominalTpcEntrancePointNeg().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 108       disty = pair->Track1()->V0()->NominalTpcEntrancePointNeg().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 109       distz = pair->Track1()->V0()->NominalTpcEntrancePointNeg().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 110       double distNeg = sqrt(distx*distx + disty*disty + distz*distz);
 111
 112       double distExitX = pair->Track1()->V0()->NominalTpcExitPointPos().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
 113       double distExitY = pair->Track1()->V0()->NominalTpcExitPointPos().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
 114       double distExitZ = pair->Track1()->V0()->NominalTpcExitPointPos().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
 115       double distExitPos = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);
 116
 117       distExitX = pair->Track1()->V0()->NominalTpcExitPointNeg().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
 118       distExitY = pair->Track1()->V0()->NominalTpcExitPointNeg().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
 119       distExitZ = pair->Track1()->V0()->NominalTpcExitPointNeg().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
 120       double distExitNeg = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);
 121
 122       tempTPCEntrancePos = distPos > fDTPCMin;
 123       tempTPCEntranceNeg = distNeg > fDTPCMin;
 124       tempTPCExitPos = distExitPos > fDTPCExitMin;
 125       tempTPCExitNeg = distExitNeg > fDTPCExitMin;
 126     }
 127
 128   if(!tempTPCEntrancePos || !tempTPCEntranceNeg || !tempTPCExitPos || !tempTPCExitNeg) return false;
 129
 130   //reject merged trakcs in TPC
 131
 132
 133   //temp = dist > fDTPCMin;
 134   //koniec kopii
 135   //if(!temp)
 136   //return false;
 137
 138   //kopia z AliFemtoShareQualityPairCut.cxx
 139   Int_t nh = 0;
 140   Int_t an = 0;
 141   Int_t ns = 0;
 142
 143
 144
 145   if ((fShareFractionMax < 1.0) && ( fShareQualityMax < 1.0)) {
 146     for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersPos().GetNbits(); imap++) {
 147       // If both have clusters in the same row
 148       if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) &&
 149           pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 150         // Do they share it ?
 151         if (pair->Track1()->V0()->TPCsharingPos().TestBitNumber(imap) &&
 152             pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
 153           {
 154             //    cout << "A shared cluster !!!" << endl;
 155             //  cout << "imap idx1 idx2 "
 156             //       << imap << " "
 157             //       << tP1idx[imap] << " " << tP2idx[imap] << endl;
 158             an++;
 159             nh+=2;
 160             ns+=2;
 161           }
 162
 163         // Different hits on the same padrow
 164         else {
 165           an--;
 166           nh+=2;
 167         }
 168       }
 169       else if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) ||
 170                pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 171         // One track has a hit, the other does not
 172         an++;
 173         nh++;
 174       }
 175     }
 176
 177     Float_t hsmval = 0.0;
 178     Float_t hsfval = 0.0;
 179
 180     if (nh >0) {
 181       hsmval = an*1.0/nh;
 182       hsfval = ns*1.0/nh;
 183     }
 184     //  if (hsmval > -0.4) {
 185     //   cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
 186     //        << (pair->Track1()->Track()) << " "
 187     //        << (pair->Track2()->Track()) << endl;
 188     //   cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
 189     //  }
 190     //   if (hsfval > 0.0) {
 191     //     cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << "    " << hsmval << " " << an << " " << nh << endl;
 192     //   }
 193
 194     temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
 195     if(!temp) return false;
 196
 197     nh = 0;
 198     an = 0;
 199     ns = 0;
 200
 201     for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersNeg().GetNbits(); imap++) {
 202       // If both have clusters in the same row
 203       if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) &&
 204           pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 205         // Do they share it ?
 206         if (pair->Track1()->V0()->TPCsharingNeg().TestBitNumber(imap) &&
 207             pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
 208           {
 209             //    cout << "A shared cluster !!!" << endl;
 210             //  cout << "imap idx1 idx2 "
 211             //       << imap << " "
 212             //       << tP1idx[imap] << " " << tP2idx[imap] << endl;
 213             an++;
 214             nh+=2;
 215             ns+=2;
 216           }
 217
 218         // Different hits on the same padrow
 219         else {
 220           an--;
 221           nh+=2;
 222         }
 223       }
 224       else if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) ||
 225                pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 226         // One track has a hit, the other does not
 227         an++;
 228         nh++;
 229       }
 230     }
 231
 232     hsmval = 0.0;
 233     hsfval = 0.0;
 234
 235     if (nh >0) {
 236       hsmval = an*1.0/nh;
 237       hsfval = ns*1.0/nh;
 238     }
 239     //  if (hsmval > -0.4) {
 240     //   cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
 241     //        << (pair->Track1()->Track()) << " "
 242     //        << (pair->Track2()->Track()) << endl;
 243     //   cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
 244     //  }
 245     //   if (hsfval > 0.0) {
 246     //     cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << "    " << hsmval << " " << an << " " << nh << endl;
 247     //   }
 248
 249     temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
 250
 251
 252   }
 253   else
 254     temp = true;
 255   //koniec kopii
 256
 257
 258
 259
 260
 261   return temp;
 262 }
 263 //__________________
 264 AliFemtoString AliFemtoV0TrackPairCut::Report(){
 265   // Prepare the report from the execution
 266   string stemp = "AliFemtoV0 Pair Cut - remove shared and split pairs\n";  char ctemp[100];
 267   snprintf(ctemp , 100, "Number of pairs which passed:\t%ld  Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
 268   stemp += ctemp;
 269   AliFemtoString returnThis = stemp;
 270   return returnThis;}
 271 //__________________
 272
 273 void AliFemtoV0TrackPairCut::SetV0Max(Double_t aV0Max) {
 274   fV0Max = aV0Max;
 275 }
 276
 277 Double_t AliFemtoV0TrackPairCut::GetAliFemtoV0Max() const {
 278   return fV0Max;
 279 }
 280
 281
 282 TList *AliFemtoV0TrackPairCut::ListSettings()
 283 {
 284   // return a list of settings in a writable form
 285   TList *tListSetttings = new TList();
 286   char buf[200];
 287   snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharequalitymax=%f", fV0Max);
 288   snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharefractionmax=%f", fShareFractionMax);
 289   tListSetttings->AddLast(new TObjString(buf));
 290
 291   return tListSetttings;
 292 }
 293
 294 void     AliFemtoV0TrackPairCut::SetRemoveSameLabel(Bool_t aRemove)
 295 {
 296   fRemoveSameLabel = aRemove;
 297 }
 298
 299 void AliFemtoV0TrackPairCut::SetTPCOnly(Bool_t tpconly)
 300 {
 301   fTrackTPCOnly = tpconly;
 302 }
 303
 304 void AliFemtoV0TrackPairCut::SetShareQualityMax(Double_t aShareQualityMax) {
 305   fShareQualityMax = aShareQualityMax;
 306 }
 307
 308 void AliFemtoV0TrackPairCut::SetShareFractionMax(Double_t aShareFractionMax) {
 309   fShareFractionMax = aShareFractionMax;
 310 }
 311
 312 void AliFemtoV0TrackPairCut::SetDataType(AliFemtoDataType type)
 313 {
 314   fDataType = type;
 315 }
 316
 317 void AliFemtoV0TrackPairCut::SetTPCEntranceSepMinimum(double dtpc)
 318 {
 319   fDTPCMin = dtpc;
 320 }
 321
 322 void AliFemtoV0TrackPairCut::SetTPCExitSepMinimum(double dtpc)
 323 {
 324   fDTPCExitMin = dtpc;
 325 }