PWGHF/hfe/AliHFEextraCuts.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
  17 // Cuts stored here:
  18 // - ITS pixels
  19 // - TPC cluster ratio
  20 // - TRD tracklets
  21 //
  22 // Authors:
  23 //   Markus Fasel <M.Fasel@gsi.de>
  24 //
  25 #include <TBits.h>
  26 #include <TClass.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TList.h>
  30 #include <TString.h>
  31 #include <TMath.h>
  32
  33 #include "AliAODTrack.h"
  34 #include "AliAODPid.h"
  35 #include "AliAODVertex.h"
  36 #include "AliESDEvent.h"
  37 #include "AliESDVertex.h"
  38 #include "AliESDtrack.h"
  39 #include "AliLog.h"
  40 #include "AliMCParticle.h"
  41 #include "AliVEvent.h"
  42 #include "AliVTrack.h"
  43 #include "AliVParticle.h"
  44 #include "AliVertexerTracks.h"
  45 #include "AliVVertex.h"
  46
  47 #include "AliHFEextraCuts.h"
  48
  49 ClassImp(AliHFEextraCuts)
  50
  51 const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
  52
  53 //______________________________________________________
  54 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
  55   AliCFCutBase(name, title),
  56   fEvent(NULL),
  57   fCutCorrelation(0),
  58   fRequirements(0),
  59   fMinNClustersTPC(0),
  60   fMinNClustersTPCPID(0),
  61   fClusterRatioTPC(0.),
  62   fMinTrackletsTRD(0),
  63   fMinNbITScls(0),
  64   fTRDtrackletsExact(0),
  65   fPixelITS(0),
  66   fTPCclusterDef(0),
  67   fTPCclusterRatioDef(0),
  68   fFractionTPCShared(-1.0),
  69   fAbsHFEImpactParamNsigmaR(kTRUE),
  70   fCheck(kFALSE),
  71   fQAlist(0x0) ,
  72   fDebugLevel(0)
  73 {
  74   //
  75   // Default Constructor
  76   //
  77   //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
  78   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
  79   memset(fIPcutParam, 0, sizeof(Float_t) * 4);
  80 }
  81
  82 //______________________________________________________
  83 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
  84   AliCFCutBase(c),
  85   fEvent(c.fEvent),
  86   fCutCorrelation(c.fCutCorrelation),
  87   fRequirements(c.fRequirements),
  88   fMinNClustersTPC(c.fMinNClustersTPC),
  89   fMinNClustersTPCPID(c.fMinNClustersTPCPID),
  90   fClusterRatioTPC(c.fClusterRatioTPC),
  91   fMinTrackletsTRD(c.fMinTrackletsTRD),
  92   fMinNbITScls(c.fMinNbITScls),
  93   fTRDtrackletsExact(c.fTRDtrackletsExact),
  94   fPixelITS(c.fPixelITS),
  95   fTPCclusterDef(c.fTPCclusterDef),
  96   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
  97   fFractionTPCShared(c.fFractionTPCShared),
  98   fAbsHFEImpactParamNsigmaR(c.fAbsHFEImpactParamNsigmaR),
  99   fCheck(c.fCheck),
 100   fQAlist(0x0),
 101   fDebugLevel(0)
 102 {
 103   //
 104   // Copy constructor
 105   // Performs a deep copy
 106   //
 107   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 108   memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 109   if(IsQAOn()){
 110     fIsQAOn = kTRUE;
 111     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 112     if(fQAlist) fQAlist->SetOwner();
 113   }
 114 }
 115
 116 //______________________________________________________
 117 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
 118   //
 119   // Assignment operator
 120   //
 121   if(this != &c){
 122     AliCFCutBase::operator=(c);
 123     fEvent = c.fEvent;
 124     fCutCorrelation = c.fCutCorrelation;
 125     fRequirements = c.fRequirements;
 126     fClusterRatioTPC = c.fClusterRatioTPC;
 127     fMinNClustersTPC = c.fMinNClustersTPC;
 128     fMinNClustersTPCPID = c.fMinNClustersTPCPID;
 129     fMinTrackletsTRD = c.fMinTrackletsTRD;
 130     fMinNbITScls = c.fMinNbITScls;
 131     fTRDtrackletsExact = c.fTRDtrackletsExact;
 132     fPixelITS = c.fPixelITS;
 133     fTPCclusterDef = c.fTPCclusterDef;
 134     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
 135     fFractionTPCShared = c.fFractionTPCShared;
 136     fAbsHFEImpactParamNsigmaR = c.fAbsHFEImpactParamNsigmaR;
 137     fCheck = c.fCheck;
 138     fDebugLevel = c.fDebugLevel;
 139     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 140     memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 141     if(IsQAOn()){
 142       fIsQAOn = kTRUE;
 143       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 144       if(fQAlist) fQAlist->SetOwner();
 145     }else fQAlist = 0x0;
 146   }
 147   return *this;
 148 }
 149
 150 //______________________________________________________
 151 AliHFEextraCuts::~AliHFEextraCuts(){
 152   //
 153   // Destructor
 154   //
 155   if(fQAlist) delete fQAlist;
 156 }
 157
 158 //______________________________________________________
 159 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
 160   //
 161   // Set Virtual event an make a copy
 162   //
 163   if (!event) {
 164     AliError("Pointer to AliVEvent !");
 165     return;
 166   }
 167   TString className(event->ClassName());
 168   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
 169     AliError("argument must point to an AliESDEvent or AliAODEvent !");
 170     return ;
 171   }
 172   fEvent = (AliVEvent*) event;
 173
 174 }
 175
 176 //______________________________________________________
 177 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
 178   //
 179   // Steering function for the track selection
 180   //
 181   TClass *type = o->IsA();
 182   AliDebug(2, Form("Object type %s", type->GetName()));
 183   if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
 184     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
 185   }
 186   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
 187 }
 188
 189 //______________________________________________________
 190 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
 191   //
 192   // Checks cuts on reconstructed tracks
 193   // returns true if track is selected
 194   // QA histograms are filled before track selection and for
 195   // selected tracks after track selection
 196   //
 197   AliDebug(1, "Called");
 198   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
 199   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
 200   // Apply cuts
 201   Float_t impactR, impactZ;
 202   Double_t hfeimpactR, hfeimpactnsigmaR;
 203   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
 204   Double_t maximpactRcut;
 205   GetImpactParameters(track, impactR, impactZ);
 206   if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 207     // Protection for PbPb
 208     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
 209     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
 210   }
 211   Int_t  nclsITS = GetITSNbOfcls(track);
 212   UInt_t nclsTPC = GetTPCncls(track);
 213   UInt_t nclsTPCPID = track->GetTPCsignalN();
 214   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 215   Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
 216   Double_t ratioTPC = GetTPCclusterRatio(track);
 217   UChar_t trdTracklets;
 218   trdTracklets = track->GetTRDntrackletsPID();
 219   UChar_t itsPixel = track->GetITSClusterMap();
 220   Int_t status1 = GetITSstatus(track, 0);
 221   Int_t status2 = GetITSstatus(track, 1);
 222   Bool_t statusL0 = CheckITSstatus(status1);
 223   Bool_t statusL1 = CheckITSstatus(status2);
 224   if(TESTBIT(fRequirements, kTPCfractionShared)) {
 225     // cut on max fraction of shared TPC clusters
 226     if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
 227   }
 228   if(TESTBIT(fRequirements, kMinImpactParamR)){
 229     // cut on min. Impact Parameter in Radial direction
 230     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
 231   }
 232   if(TESTBIT(fRequirements, kMinImpactParamZ)){
 233     // cut on min. Impact Parameter in Z direction
 234     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
 235   }
 236   if(TESTBIT(fRequirements, kMaxImpactParamR)){
 237     // cut on max. Impact Parameter in Radial direction
 238     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
 239   }
 240   if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
 241     // HFE Impact parameter cut
 242     GetMaxImpactParameterCutR(track,maximpactRcut);
 243     if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
 244   }
 245   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
 246     // cut on max. Impact Parameter in Z direction
 247     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
 248   }
 249   if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
 250     // cut on min. HFE Impact Parameter in Radial direction
 251     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
 252   }
 253   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 254     // cut on max. HFE Impact Parameter n sigma in Radial direction
 255     if(fAbsHFEImpactParamNsigmaR) {
 256 //      if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
 257       if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
 258         SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 259 //        printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
 260       }
 261     }
 262     else {
 263       if(hfeimpactnsigmaRcut>0){
 264         if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
 265           SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 266           //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
 267         }
 268       }
 269       else{
 270         if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
 271           SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 272           //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
 273         }
 274       }
 275     }
 276   }
 277   if(TESTBIT(fRequirements, kClusterRatioTPC)){
 278     // cut on min ratio of found TPC clusters vs findable TPC clusters
 279     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
 280   }
 281   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
 282     // cut on minimum number of TRD tracklets
 283     AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
 284     if(fTRDtrackletsExact){
 285       if(trdTracklets == fMinTrackletsTRD) {
 286         SETBIT(survivedCut, kMinTrackletsTRD);
 287         AliDebug(1, "Track Selected");
 288       }
 289     }else{
 290       if(trdTracklets >= fMinTrackletsTRD){
 291         SETBIT(survivedCut, kMinTrackletsTRD);
 292         AliDebug(1, "Track Selected");
 293       }
 294       //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
 295     }
 296   }
 297
 298   if(TESTBIT(fRequirements, kMinNbITScls)){
 299     // cut on minimum number of ITS clusters
 300     //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
 301     AliDebug(1, Form("Min ITS clusters: [%d|%d]\n", fMinNbITScls, nclsITS));
 302     if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
 303   }
 304
 305   if(TESTBIT(fRequirements, kMinNClustersTPC)){
 306     // cut on minimum number of TPC tracklets
 307     //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
 308     AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
 309     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
 310   }
 311   if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
 312     AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
 313     if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
 314   }
 315   if(TESTBIT(fRequirements, kPixelITS)){
 316     // cut on ITS pixel layers
 317     AliDebug(1, "ITS cluster Map: ");
 318     //PrintBitMap(itsPixel);
 319     switch(fPixelITS){
 320       case kFirst:
 321         AliDebug(2, "First");
 322         //printf("First\n");
 323               if(TESTBIT(itsPixel, 0))
 324                 SETBIT(survivedCut, kPixelITS);
 325         else
 326                 if(fCheck && !statusL0)
 327                   SETBIT(survivedCut, kPixelITS);
 328                     break;
 329       case kSecond:
 330         //printf("Second\n");
 331         AliDebug(2, "Second");
 332               if(TESTBIT(itsPixel, 1))
 333                 SETBIT(survivedCut, kPixelITS);
 334         else
 335                 if(fCheck && !statusL1)
 336                   SETBIT(survivedCut, kPixelITS);
 337                     break;
 338       case kBoth:
 339         //printf("Both\n");
 340         AliDebug(2, "Both");
 341               if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
 342                       SETBIT(survivedCut, kPixelITS);
 343               else
 344           if(fCheck && !(statusL0 && statusL1))
 345                         SETBIT(survivedCut, kPixelITS);
 346               break;
 347       case kAny:
 348         //printf("Any\n");
 349         AliDebug(2, "Any");
 350               if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
 351                 SETBIT(survivedCut, kPixelITS);
 352         else
 353                 if(fCheck && !(statusL0 || statusL1))
 354                     SETBIT(survivedCut, kPixelITS);
 355                     break;
 356       case kExclusiveSecond:
 357         AliDebug(2, "Exlusive second");
 358         if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
 359           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
 360             SETBIT(survivedCut, kPixelITS);
 361         } else {
 362           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
 363             SETBIT(survivedCut, kPixelITS);
 364         }
 365         break;
 366       case kNone:
 367         // No cut applied, set as survived
 368         SETBIT(survivedCut, kPixelITS);
 369         break;
 370       default:
 371         // default, defined as no cut applied
 372         AliDebug(2, "None");
 373         SETBIT(survivedCut, kPixelITS);
 374         break;
 375     }
 376     AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
 377   }
 378
 379   if(TESTBIT(fRequirements, kTOFPID)){
 380     // cut on TOF pid
 381     if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
 382   }
 383
 384   if(TESTBIT(fRequirements, kTOFmismatch)){
 385     // cut on TOF mismatch
 386     if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
 387   }
 388
 389   if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
 390       // cut on TPC PID cleanup
 391       Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
 392       if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
 393   }
 394
 395   if(TESTBIT(fRequirements, kEMCALmatch)){
 396     if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
 397   }
 398
 399   if(TESTBIT(fRequirements, kRejectKinkDaughter)){
 400     //printf("test daughter\n");
 401     if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
 402     //else printf("Found kink daughter\n");
 403   }
 404
 405   if(TESTBIT(fRequirements, kRejectKinkMother)){
 406     //printf("test mother\n");
 407     if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
 408     //else printf("Found kink mother\n");
 409   }
 410
 411   if(fRequirements == survivedCut){
 412     //
 413     // Track selected
 414     //
 415     AliDebug(2, "Track Survived cuts\n");
 416     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
 417     return kTRUE;
 418   }
 419   AliDebug(2, "Track cut");
 420   if(IsQAOn()) FillCutCorrelation(survivedCut);
 421   return kFALSE;
 422 }
 423
 424 //______________________________________________________
 425 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
 426   //
 427   // Checks cuts on Monte Carlo tracks
 428   // returns true if track is selected
 429   // QA histograms are filled before track selection and for
 430   // selected tracks after track selection
 431   //
 432   return kTRUE; // not yet implemented
 433 }
 434
 435 //______________________________________________________
 436 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
 437   //
 438   // Fill the QA histograms for ESD tracks
 439   // Function can be called before cuts or after cut application (second argument)
 440   //
 441   Float_t impactR, impactZ;
 442   GetImpactParameters(track, impactR, impactZ);
 443   TH1 *htmp = NULL;
 444   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
 445   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
 446   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 447   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
 448   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
 449   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
 450   UChar_t itsPixel = track->GetITSClusterMap();
 451   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
 452   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 453   if(pixelHist){
 454     Double_t firstEntry = 0.5;
 455     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 456       pixelHist->Fill(firstEntry + 3);
 457     else{
 458       if(itsPixel & BIT(0)){
 459         pixelHist->Fill(firstEntry);
 460         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 461         else pixelHist->Fill(firstEntry + 4);
 462       }
 463       if(itsPixel & BIT(1)){
 464         pixelHist->Fill(firstEntry + 1);
 465         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 466       }
 467     }
 468   }
 469   // Fill histogram with the status bits
 470   TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
 471   if(hStatusBits) {
 472     hStatusBits->Fill(0);  // Fill first bin with all tracks
 473     if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
 474     if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
 475     if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
 476     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
 477   }
 478   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
 479 }
 480
 481 // //______________________________________________________
 482 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 483 //   //
 484 //   // Fill the QA histograms for MC tracks
 485 //   // Function can be called before cuts or after cut application (second argument)
 486 //   // Not yet implemented
 487 //   //
 488 // }
 489
 490 //______________________________________________________
 491 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 492   //
 493   // Fill cut correlation histograms for tracks that didn't pass cuts
 494   //
 495   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
 496   if(correlation){
 497     for(Int_t icut = 0; icut < kNcuts; icut++){
 498       if(!TESTBIT(fRequirements, icut)) continue;
 499       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 500         if(!TESTBIT(fRequirements, jcut)) continue;
 501         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 502                 correlation->Fill(icut, jcut);
 503       }
 504     }
 505   }
 506 }
 507
 508 //______________________________________________________
 509 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 510   //
 511   // Add QA histograms
 512   // For each cut a histogram before and after track cut is created
 513   // Histos before respectively after cut are stored in different lists
 514   // Additionally a histogram with the cut correlation is created and stored
 515   // in the top directory
 516   //
 517
 518   TH1 *histo1D = 0x0;
 519   TH2 *histo2D = 0x0;
 520   TString cutstr[2] = {"before", "after"};
 521
 522   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
 523   for(Int_t icond = 0; icond < 2; icond++){
 524     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
 525     fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
 526     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 527     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 528     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
 529     fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
 530     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 531     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 532     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
 533     fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
 534     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 535     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 536     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
 537     fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
 538     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 539     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 540     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
 541     fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
 542     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
 543     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 544     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
 545     fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
 546     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 547     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 548     Int_t first = histo1D->GetXaxis()->GetFirst();
 549     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 550     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 551       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 552     qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
 553     fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
 554     histo1D->GetXaxis()->SetTitle("Track Status Bit");
 555     histo1D->GetYaxis()->SetTitle("Number of tracks");
 556     TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
 557     for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
 558     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcdEdxClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters for dEdx calculation", 161, 0, 160)), 7 + icond * fgkNQAhistos);
 559     fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
 560     histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
 561     histo1D->GetYaxis()->SetTitle("counts");
 562   }
 563   // Add cut correlation
 564   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
 565   fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
 566   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters" "Fraction Shared TPC clusters", "TOFPID", "No TOFmismatch", "EMCALmatch", "ImpactParam"};
 567   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 568   for(Int_t icut = 0; icut < kNcuts; icut++){
 569     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 570     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 571   }
 572 }
 573
 574 //______________________________________________________
 575 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 576   for(Int_t ibit = 32; ibit--; )
 577     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 578   printf("\n");
 579 }
 580
 581 //______________________________________________________
 582 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 583   //
 584   // Check whether ITS area is dead
 585   //
 586   Bool_t status;
 587   switch(itsStatus){
 588     case 2: status = kFALSE; break;
 589     case 3: status = kFALSE; break;
 590     case 7: status = kFALSE; break;
 591     default: status = kTRUE;
 592   }
 593   return status;
 594 }
 595
 596 //______________________________________________________
 597 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
 598         //
 599         // Check ITS layer status
 600         //
 601         Int_t status = 0;
 602         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 603                 Int_t det;
 604                 Float_t xloc, zloc;
 605                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 606                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
 607         }
 608         return status;
 609 }
 610
 611
 612 //______________________________________________________
 613 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
 614   //
 615   // Checks if number of shared bits is above 1
 616   //
 617   Int_t nsharebit = 1;
 618   const TBits *shared;
 619   if(track->IsA() == AliESDtrack::Class())
 620     shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
 621   else if(track->IsA() == AliAODTrack::Class())
 622     shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
 623    else
 624     return 0;
 625
 626         if(shared->CountBits() >= 2) nsharebit=1;
 627   else nsharebit=0;
 628
 629   return nsharebit;
 630 }
 631
 632 //______________________________________________________
 633 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
 634         //
 635         // Get Number of findable clusters in the TPC
 636         //
 637         Int_t nClusters = 0; // in case no Information available consider all clusters findable
 638   TClass *type = track->IsA();
 639   if(TESTBIT(fTPCclusterDef, kFoundIter1)){
 640                 if(type == AliESDtrack::Class()){
 641       AliDebug(2, ("Using def kFoundIter1"));
 642       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 643       nClusters = esdtrack->GetTPCNclsIter1();
 644     } else {
 645       AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
 646     }
 647   } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
 648     AliDebug(2, ("Using def kCrossedRows"));
 649     nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
 650   } else{
 651     AliDebug(2, ("Using def kFound"));
 652           nClusters = track->GetTPCNcls();
 653   }
 654         return nClusters;
 655 }
 656
 657 //______________________________________________________
 658 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
 659   //
 660   // Get fraction of shared TPC clusters
 661   //
 662   Float_t fracClustersTPCShared = 0.0;
 663   Int_t nClustersTPC = track->GetTPCNcls();
 664   Int_t nClustersTPCShared = 0;
 665   TClass *type = track->IsA();
 666   if(type == AliESDtrack::Class()){
 667     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 668     nClustersTPCShared = esdtrack->GetTPCnclsS();
 669   } else if(type == AliESDtrack::Class()){
 670     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 671     const TBits &shared = aodtrack->GetTPCSharedMap();
 672     nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
 673   }
 674   if (nClustersTPC!=0) {
 675            fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
 676   }
 677   return fracClustersTPCShared;
 678 }
 679
 680 //______________________________________________________
 681 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
 682   //
 683   // Get Ratio of found / findable clusters for different definitions
 684   // Only implemented for ESD tracks
 685   //
 686   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
 687   TClass *type = track->IsA();
 688   if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
 689     AliDebug(2, "Using ratio def kCROverFindable");
 690     clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
 691   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
 692     AliDebug(2, "Using ratio def kFoundOverCR");
 693     clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
 694   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
 695     if(type == AliESDtrack::Class()){
 696       AliDebug(2, "Using ratio def kFoundOverFindableIter1");
 697       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 698       clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
 699     } else {
 700       AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
 701       clusterRatio = 1.;
 702     }
 703   } else {
 704     AliDebug(2, "Using ratio def kFoundOverFindable");
 705     clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
 706   }
 707   return clusterRatio;
 708 }
 709
 710 //______________________________________________________
 711 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
 712   //
 713   // Get impact parameter
 714   //
 715   TClass *type = track->IsA();
 716   if(type == AliESDtrack::Class()){
 717     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 718     esdtrack->GetImpactParameters(radial, z);
 719   }
 720   else if(type == AliAODTrack::Class()){
 721     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 722     if(aodtrack){
 723       Double_t xyz[3];
 724       aodtrack->XYZAtDCA(xyz);
 725       z = xyz[2];
 726       radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
 727     }
 728   }
 729 }
 730 //______________________________________________________
 731 Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
 732   //
 733   // Is kink Daughter
 734   //
 735   TClass *type = track->IsA();
 736   if(type == AliESDtrack::Class()){
 737     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 738     if(!esdtrack) return kFALSE;
 739     if(esdtrack->GetKinkIndex(0)<=0) return kFALSE;
 740     else return kTRUE;
 741
 742   }
 743   else if(type == AliAODTrack::Class()){
 744     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 745     if(aodtrack){
 746       AliAODVertex *aodvertex = aodtrack->GetProdVertex();
 747       if(!aodvertex) return kFALSE;
 748       if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
 749       else return kFALSE;
 750     }
 751     else return kFALSE;
 752   }
 753
 754   return kFALSE;
 755 }
 756 //______________________________________________________
 757 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
 758   //
 759   // Is kink Mother: only for ESD since need to loop over vertices for AOD
 760   //
 761   //
 762
 763   TClass *type = track->IsA();
 764   if(type == AliESDtrack::Class()){
 765     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 766     if(!esdtrack) return kFALSE;
 767     if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
 768     else return kFALSE;
 769   }
 770
 771   return kFALSE;
 772
 773 }
 774 //______________________________________________________
 775 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
 776   //
 777   // Get ITS nb of clusters
 778   //
 779   TClass *type = track->IsA();
 780   if(type == AliESDtrack::Class()){
 781     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 782     return esdtrack->GetITSclusters(0);
 783
 784   }
 785   else if(type == AliAODTrack::Class()){
 786     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 787     if(aodtrack){
 788       return  aodtrack->GetITSNcls();
 789     }
 790   }
 791
 792   return 0;
 793 }
 794 //______________________________________________________
 795 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
 796         //
 797         // Get HFE impact parameter (with recalculated primary vertex)
 798         //
 799         dcaxy=0;
 800         dcansigmaxy=0;
 801   if(!fEvent){
 802     AliDebug(1, "No Input event available\n");
 803     return;
 804   }
 805   const Double_t kBeampiperadius=3.;
 806   TString type = track->IsA()->GetName();
 807   Double_t dcaD[2]={-999.,-999.};
 808   Double_t covD[3]={-999.,-999.,-999.};
 809   Float_t dcaF[2]={-999.,-999.};
 810   Float_t covF[3]={-999.,-999.,-999.};
 811
 812   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 813   if(!esdevent) {
 814     AliDebug(1, "No esd event available\n");
 815     return;
 816   }
 817   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
 818   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
 819    // recalculate primary vertex for peri. and pp
 820    AliVertexerTracks vertexer(fEvent->GetMagneticField());
 821    vertexer.SetITSMode();
 822    vertexer.SetMinClusters(4);
 823          Int_t skipped[2];
 824    skipped[0] = track->GetID();
 825    vertexer.SetSkipTracks(1,skipped);
 826
 827  //----diamond constraint-----------------------------
 828    vertexer.SetConstraintOn();
 829    Float_t diamondcovxy[3];
 830    esdevent->GetDiamondCovXY(diamondcovxy);
 831    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
 832    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 833    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 834    vertexer.SetVtxStart(diamond);
 835    delete diamond; diamond=NULL;
 836  //----------------------------------------------------
 837
 838    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 839
 840    // Getting the DCA
 841    // Propagation always done on a working copy to not disturb the track params of the original track
 842    AliESDtrack *esdtrack = NULL;
 843    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 844      // Case ESD track: take copy constructor
 845      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 846      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 847    } else {
 848     // Case AOD track: take different constructor
 849     esdtrack = new AliESDtrack(track);
 850    }
 851    if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
 852     // protection
 853     dcaxy = dcaD[0];
 854     if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
 855     if(!covD[0]) dcansigmaxy = -49.;
 856    }
 857    delete esdtrack;
 858    delete vtxESDSkip;
 859   }
 860   else{
 861    AliESDtrack *esdtrack = NULL;
 862    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 863     // Case ESD track: take copy constructor
 864     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 865     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 866    } else {
 867     // Case AOD track: take different constructor
 868     esdtrack = new AliESDtrack(track);
 869    }
 870    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
 871    dcaxy = dcaF[0];
 872    if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
 873    if(!covF[0]) dcansigmaxy = -49.;
 874    delete esdtrack;
 875   }
 876 }
 877
 878 //______________________________________________________
 879 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
 880         //
 881         // Get HFE impact parameter (with recalculated primary vertex)
 882         //
 883   if(!fEvent){
 884     AliDebug(1, "No Input event available\n");
 885     return;
 886   }
 887   const Double_t kBeampiperadius=3.;
 888   TString type = track->IsA()->GetName();
 889   Float_t dcaF[2]={-999.,-999.};
 890   Float_t covF[3]={-999.,-999.,-999.};
 891
 892   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 893   if(!esdevent) {
 894     AliDebug(1, "No esd event available\n");
 895     return;
 896   }
 897   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
 898   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
 899    // recalculate primary vertex for peri. and pp
 900    AliVertexerTracks vertexer(fEvent->GetMagneticField());
 901    vertexer.SetITSMode();
 902    vertexer.SetMinClusters(4);
 903          Int_t skipped[2];
 904    skipped[0] = track->GetID();
 905    vertexer.SetSkipTracks(1,skipped);
 906
 907  //----diamond constraint-----------------------------
 908    vertexer.SetConstraintOn();
 909    Float_t diamondcovxy[3];
 910    esdevent->GetDiamondCovXY(diamondcovxy);
 911    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
 912    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 913    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 914    vertexer.SetVtxStart(diamond);
 915    delete diamond; diamond=NULL;
 916  //----------------------------------------------------
 917
 918    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 919
 920    // Getting the DCA
 921    // Propagation always done on a working copy to not disturb the track params of the original track
 922    AliESDtrack *esdtrack = NULL;
 923    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 924      // Case ESD track: take copy constructor
 925      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 926      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 927    } else {
 928     // Case AOD track: take different constructor
 929     esdtrack = new AliESDtrack(track);
 930    }
 931    if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
 932    delete esdtrack;
 933    delete vtxESDSkip;
 934   }
 935   else{
 936    AliESDtrack *esdtrack = NULL;
 937    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 938     // Case ESD track: take copy constructor
 939     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 940     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 941    } else {
 942     // Case AOD track: take different constructor
 943     esdtrack = new AliESDtrack(track);
 944    }
 945    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
 946    dcaD[0] = dcaF[0];
 947    dcaD[1] = dcaF[1];
 948    covD[0] = covF[0];
 949    covD[1] = covF[1];
 950    covD[2] = covF[2];
 951    delete esdtrack;
 952   }
 953 }
 954
 955 //______________________________________________________
 956 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
 957         //
 958         // Get HFE impact parameter cut(pt dependent)
 959         //
 960
 961   TString type = track->IsA()->GetName();
 962   if(!type.CompareTo("AliESDtrack")){
 963     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 964     if(!esdtrack) return;
 965     Double_t pt = esdtrack->Pt();
 966     hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
 967     hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
 968   }
 969 }
 970 //______________________________________________________
 971 void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
 972         //
 973         // Get max impact parameter cut r (pt dependent)
 974         //
 975
 976   TString type = track->IsA()->GetName();
 977   if(!type.CompareTo("AliESDtrack")){
 978     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 979     if(!esdtrack) return;
 980     Double_t pt = esdtrack->Pt();
 981     if(pt > 0.15) {
 982       maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
 983     }
 984     else maximpactRcut = 9999999999.0;
 985   }
 986 }